BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
QCMCFilter.cxx
Go to the documentation of this file.
1//QCMCFilter-00-00-02 Jan. 2013 Dan Ambrose
2//Based on QCMCReweightProc program by Werner Sun
3//This version has a corrected Y input from original version
4#include "QCMCFilterAlg/QCMCFilter.h"
5#include "QCMCFilterAlg/Dalitz.h"
6#include "QCMCFilterAlg/D0To2pip2pim.h"
7#include "QCMCFilterAlg/D0Topippim2pi0.h"
8#include "QCMCFilterAlg/D0ToKSpipipi0.h"
9#include "QCMCFilterAlg/D0ToKSpipi.h"
10#include "QCMCFilterAlg/D0ToKLpipi.h"
11#include "QCMCFilterAlg/D0ToKLpipi2018.h"
12#include "QCMCFilterAlg/D0ToKSpipi2018.h"
13#include "QCMCFilterAlg/D0ToKKpipi.h"
14#include "QCMCFilterAlg/D0TopiKpipi.h"
15#include "QCMCFilterAlg/D0ToKpipipiLHCb.h"
16//#include "models/dcs_LHCb_BESIIIIMP.h"
17//#include "models/cf_LHCb_BESIIIIMP.h"
18
19#include "GaudiKernel/MsgStream.h"
20#include "GaudiKernel/AlgFactory.h"
21#include "GaudiKernel/ISvcLocator.h"
22#include "GaudiKernel/SmartDataPtr.h"
23#include "GaudiKernel/IDataProviderSvc.h"
24#include "GaudiKernel/PropertyMgr.h"
25#include "GaudiKernel/Bootstrap.h"
26#include "GaudiKernel/RegistryEntry.h"
27
28#include "TMath.h"
29
30#include <cmath>
31#include "HepPDT/ParticleDataTable.hh"
32#include "HepPDT/ParticleData.hh"
33//#include "PartPropSvc/PartPropSvc.h"
34#include "GaudiKernel/IPartPropSvc.h"
35
36// Monte-Carlo data
37#include "McTruth/McParticle.h"
38#include "McTruth/MdcMcHit.h"
39
40///////////////////
41#include "EventModel/EventModel.h"
42#include "EventModel/Event.h"
43#include "EventModel/EventHeader.h"
44#include "EvtGenBase/EvtComplex.hh"
45#include "EvtGenBase/EvtVector4R.hh"
46
47#include "EvtRecEvent/EvtRecEvent.h"
48
49#include "CLHEP/Random/RandFlat.h"
50#include "CLHEP/Matrix/Vector.h"
51#include "CLHEP/Matrix/Matrix.h"
52#include "CLHEP/Matrix/SymMatrix.h"
53#include "CLHEP/Vector/ThreeVector.h"
54#include "CLHEP/Vector/LorentzVector.h"
55#include "CLHEP/Vector/TwoVector.h"
56using CLHEP::HepVector;
57using CLHEP::Hep3Vector;
58using CLHEP::Hep2Vector;
59using CLHEP::HepLorentzVector;
60
61#include <vector>
62#include <complex>
63
64//
65// constants, enums and typedefs
66//
67// Partilcle masses
68const double xmpi0 = 0.134976; // BOSS 6.4.1 pdt.table value
69const double xmeta = 0.54784; // PDG08 = 0.547853, BOSS 6.4.1 pdt.table = 0.54784
70const double xmkaon = 0.49368;
71const double xmpion = 0.13957;
72const double xmk0 = 0.49761;
73const double xmrho = 0.77549;
74const double xmd0 = 1.86484; // PDG08
75
76const double PI = 3.1415926; //pi
77
78// Antiparticles have negative ID.
79
80static const int kPsi3770ID = 30443 ;
81static const int kD0ID = 421 ;
82static const int kD0BarID = -421 ;
83static const int kDpID = 411 ;
84static const int kDmID = -411 ;
85static const int kGammaID = 22 ;
86static const int kGammaFSRID = -22 ;
87static const int kPiPlusID = 211 ;
88static const int kPiMinusID = -211 ;
89static const int kPi0ID = 111 ;
90static const int kEtaID = 221 ;
91static const int kEtaPrimeID = 331 ;
92static const int kF0ID = 9010221;
93static const int kFPrime0ID = 10221 ;
94static const int kF0m1500ID = 50221;
95static const int kF2ID = 225;
96static const int kA00ID = 10111 ;
97static const int kA0PlusID = 10211 ;
98static const int kA0MinusID = -10211 ;
99static const int kRhoPlusID = 213 ;
100static const int kRhoMinusID = -213 ;
101static const int kRho0ID = 113 ;
102static const int kRho2SPlusID = 30213 ;
103static const int kRho2SMinusID = -30213 ;
104static const int kRho2S0ID = 30113 ;
105static const int kA1PlusID = 20213 ;
106static const int kA1MinusID = -20213 ;
107static const int kA10ID = 20113 ;
108static const int kOmegaID = 223 ;
109static const int kPhiID = 333 ;
110static const int kKPlusID = 321 ;
111static const int kKMinusID = -321 ;
112static const int kK0SID = 310 ;
113static const int kK0LID = 130 ;
114static const int kK0ID = 311 ;
115static const int kK0BarID = -311 ;
116static const int kKStarPlusID = 323 ;
117static const int kKStarMinusID = -323 ;
118static const int kKStar0ID = 313 ;
119static const int kKStar0BarID = -313 ;
120static const int kKDPStar0ID = 30313 ;
121static const int kKDPStar0BarID = -30313 ;
122static const int kK0Star0ID = 10311;
123static const int kK0Star0BarID = -10311;
124static const int kK0StarPlusID = 10321;
125static const int kK0StarMinusID = -10321;
126static const int kK1PlusID = 10323 ;
127static const int kK1MinusID = -10323 ;
128static const int kK10ID = 10313 ;
129static const int kK10BarID = -10313 ;
130static const int kK1PrimePlusID = 20323 ;
131static const int kK1PrimeMinusID = -20323 ;
132static const int kK1Prime0ID = 20313 ;
133static const int kK1Prime0BarID = -20313 ;
134static const int kK2StarPlusID = 325 ;
135static const int kK2StarMinusID = -325 ;
136static const int kK2Star0ID = 315 ;
137static const int kK2Star0BarID = -315 ;
138static const int kEMinusID = 11 ;
139static const int kEPlusID = -11 ;
140static const int kMuMinusID = 13 ;
141static const int kMuPlusID = -13 ;
142static const int kNuEID = 12 ;
143static const int kNuEBarID = -12 ;
144static const int kNuMuID = 14 ;
145static const int kNuMuBarID = -14 ;
146
147//D0 decay modes
148static const int kFlavoredCF_0 = 0;
149static const int kFlavoredCFBar_0 = 1;
150static const int kFlavoredCF_1 = 2;
151static const int kFlavoredCFBar_1 = 3;
152static const int kFlavoredCF_3 = 4;
153static const int kFlavoredCFBar_3 = 5;
154static const int kFlavoredCS = 6;
155static const int kFlavoredCSBar = 7;
156static const int kSLPlus = 8;
157static const int kSLMinus = 9;
158static const int kCPPlus = 10;
159static const int kCPMinus = 11;
160static const int kDalitz = 12;
161static const int k2Pip2Pim = 13;
162static const int kPipPim2Pi0 = 14;
163static const int kK0PiPiPi0 = 15;
164static const int kKKPiPi = 16;
165static const int kNDecayTypes = 17;
166
167
168
169//Varibales to keep track of the Events
170int m_nUnknownEvents = 0;
171int m_nUnknownDecays= 0;
172int m_nD0D0barEvents= 0;
173int m_nD0bar = 0;
174int m_nDpDmEvents = 0;
175int m_nD0D0barDiscarded = 0;
176int m_nDpDmDiscarded = 0;
177int m_nCPPlus = 0;
178int m_nCPMinus = 0;
179int m_nFlavoredCFD0_0 = 0;
180int m_nFlavoredCFD0_1 = 0;
181int m_nFlavoredCFD0_3 = 0;
182int m_nFlavoredCSD0 = 0;
183int m_nFlavoredDCSD0_0 = 0;
184int m_nFlavoredDCSD0_1 = 0;
185int m_nFlavoredDCSD0_3 = 0;
186int m_nSL = 0;
187int m_nDalitz = 0;
188int m_n2Pip2Pim = 0;
189int m_nPipPim2Pi0 = 0;
190int m_nK0PiPiPi0 = 0;
191int m_nKKPiPi = 0;
192double m_dalitzNumer1 = 0;
193double m_dalitzNumer2 = 0;
194double m_dalitzDenom = 0;
195double dalitzNumer1_fil = 0;
196double dalitzNumer2_fil = 0;
197double dalitzDenom_fil = 0;
198double m_2Pip2PimNumer1 = 0;
199double m_2Pip2PimNumer2 = 0;
200double m_2Pip2PimDenom = 0;
201double Pip2Pim2Numer1_fil = 0;
202double Pip2Pim2Numer2_fil = 0;
203double Pip2Pim2Denom_fil = 0;
204double m_PipPim2Pi0Numer1 = 0;
205double m_PipPim2Pi0Numer2 = 0;
206double m_PipPim2Pi0Denom = 0;
207double PipPim2Pi0Numer1_fil = 0;
208double PipPim2Pi0Numer2_fil = 0;
209double PipPim2Pi0Denom_fil = 0;
210double m_KSPiPiPi0Numer1 = 0;
211double m_KSPiPiPi0Numer2 = 0;
212double m_KSPiPiPi0Denom = 0;
213double KSPiPiPi0Numer1_fil = 0;
214double KSPiPiPi0Numer2_fil = 0;
215double KSPiPiPi0Denom_fil = 0;
216HepSymMatrix m_weights(18,0);
217double m_rwsCF=0.;
218double m_rwsCF_0=0.;
219double m_rwsCF_1=0.;
220double m_rwsCF_3=0.;
221double m_rwsCS=0.;
222double m_deltaCF=0.;
223double m_deltaCF_0=0.;
224double m_deltaCF_1=0.;
225double m_deltaCF_3=0.;
226double m_deltaCS=0.;
227HepMatrix m_modeCounter(19,19,0);
228HepMatrix m_keptModeCounter(19,19,0);
229
230int m_rho_flag;
231
232/////////////////////////////////////////////////////////////////
233DECLARE_COMPONENT( QCMCFilter )
234
235 QCMCFilter::QCMCFilter(const std::string& name, ISvcLocator* pSvcLocator) :
236 Algorithm(name, pSvcLocator){
237 //Declare the properties
238 declareProperty("Debug", m_debug= true );
239 declareProperty("x", m_x= 0.00407);
240 declareProperty("y", m_y= 0.00647); //For values outside of (-.11 , .06) measured y value may vary from input y
241
242 declareProperty("rForCFModes", m_rCF= 0.0621); //should be kept near this. based on MC for d0->kpi and d0bar->kpi 0.0621 = sqrt( .00015/.0389)
243 declareProperty("zForCFModes", m_zCF= 2.0);
244 declareProperty("wSignForCFModes", m_wCFSign= true);
245 declareProperty("rForCFMode_0", m_rCF_0= 0.0587); //should be kept near this. based on MC for d0->kpi and d0bar->kpi 0.0621 = sqrt( .00015/.0389)
246 declareProperty("RForCFMode_0", m_RCF_0= 1.); //should be kept near this.
247 declareProperty("zForCFMode_0", m_zCF_0= 1.9585);
248 declareProperty("dForCFMode_0", m_dCF_0= 191.7*3.14159/180.);
249 declareProperty("wSignForCFMode_0", m_wCFSign_0 = true);
250 declareProperty("rForCFModes_1", m_rCF_1= 0.044); //should be kept near this. based on MC for d0->kpi and d0bar->kpi 0.0621 = sqrt( .00015/.0389)
251 declareProperty("RForCFMode_1", m_RCF_1= 0.79); //should be kept near this.
252 declareProperty("zForCFModes_1", m_zCF_1= 1.9225);
253 declareProperty("dForCFMode_1", m_dCF_1= 196.*3.14159/180.);
254 declareProperty("wSignForCFModes_1", m_wCFSign_1= true);
255 declareProperty("rForCFModes_3", m_rCF_3= 0.0549); //should be kept near this. based on MC for d0->kpi and d0bar->kpi 0.0621 = sqrt( .00015/.0389)
256 declareProperty("RForCFMode_3", m_RCF_3= 0.43); //should be kept near this.
257 declareProperty("zForCFModes_3", m_zCF_3= 1.2313);
258 declareProperty("dForCFMode_3", m_dCF_3= 161.*3.14159/180.);
259 declareProperty("wSignForCFModes_3", m_wCFSign_3= true);
260
261 declareProperty("rForCSModes", m_rCS= 1.0);
262 declareProperty("zForCSModes", m_zCS= -2.0);
263 declareProperty("wSignForCSModes", m_wCSSign= true);
264 declareProperty("DplusDminusWeight", m_dplusDminusWeight= 0.5);
265 declareProperty("dalitzModel", m_dalitzModel= 0); // 0 CLEO-c, 1 Babar, 2 Belle
266 declareProperty("UseNewWeights", m_useNewWeights= true);
267 declareProperty("InvertSelection", m_invertSelection= false);
268 }
269
270//***************************************************************
271
272StatusCode QCMCFilter::initialize() {
273 MsgStream log(msgSvc(), name());
274
275 log << MSG::INFO << "in initialize()" << endmsg;
276
277 // Get the Particle Properties Service
278 IPartPropSvc* p_PartPropSvc;
279 static const bool CREATEIFNOTTHERE(true);
280 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
281 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc)
282 {
283 log << MSG::ERROR << " Could not initialize Particle Properties Service" << endmsg;
284 return PartPropStatus;
285 }
286 m_particleTable = p_PartPropSvc->PDT();
287
288
289 //Calculates parameters here based on Parameter Input by User
290 // (e.g. run expensive algorithms that are based on parameters
291 // specified by user at run-time)
292
293 //m_y = (m_y - 0.008)/0.292;
294 //corrects the y input so that we get the desired measurable y from CP:Semi-leptonic doubletag method
295 // This equation is a result of solving equations for y dependence.
296 //If parent MC Branching ratio's this equation will need resolved. January 2013 -DA
297 double x = m_x ;
298 double y = m_y ;
299 double rCF_0 = m_rCF_0 ;
300 double dCF_0 = m_dCF_0 ;
301 double RCF_0 = m_RCF_0 ;
302 //double zCF_0 = m_zCF_0 ;
303 double zCF_0 = 2*cos(dCF_0) ;
304 double rCF2_0 = rCF_0 * rCF_0 ;
305 double zCF2_0 = zCF_0 * zCF_0 ;
306 double rCF_1 = m_rCF_1 ;
307 double dCF_1 = m_dCF_1 ;
308 double RCF_1 = m_RCF_1 ;
309 //double zCF_1 = m_zCF_1 ;
310 double zCF_1 = 2*cos(dCF_1) ;
311 double rCF2_1 = rCF_1 * rCF_1 ;
312 double zCF2_1 = zCF_1 * zCF_1 ;
313 double rCF_3 = m_rCF_3 ;
314 double dCF_3 = m_dCF_3 ;
315 double RCF_3 = m_RCF_3 ;
316 //double zCF_3 = m_zCF_3 ;
317 double zCF_3 = 2*cos(dCF_3) ;
318 double rCF2_3 = rCF_3 * rCF_3 ;
319 double zCF2_3 = zCF_3 * zCF_3 ;
320 double rCS = m_rCS ;
321 double zCS = m_zCS ;
322 double rCS2 = rCS * rCS ;
323 double zCS2 = zCS * zCS ;
324
325 double rmix = ( x * x + y * y ) / 2. ;
326 double wCF_0 = 2*cos(dCF_0) ;
327 double wCF_1 = 2*cos(dCF_1) ;
328 double wCF_3 = 2*cos(dCF_3) ;
329 double wCS = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 ) ;
330 double vCF_0CSPlus = ( zCF_0 * zCS + wCF_0 * wCS ) / 2. ;
331 double vCF_1CSPlus = ( zCF_1 * zCS + wCF_1 * wCS ) / 2. ;
332 double vCF_3CSPlus = ( zCF_3 * zCS + wCF_3 * wCS ) / 2. ;
333 double vCF_0CSMinus = ( zCF_0 * zCS - wCF_0 * wCS ) / 2. ;
334 double vCF_1CSMinus = ( zCF_1 * zCS - wCF_1 * wCS ) / 2. ;
335 double vCF_3CSMinus = ( zCF_3 * zCS - wCF_3 * wCS ) / 2. ;
336
337 m_deltaCF_0 = acos( zCF_0 / 2. ) * ( m_wCFSign_0 ? 1. : -1. ) ;
338 m_deltaCF_1 = acos( zCF_1 / 2. ) * ( m_wCFSign_1 ? 1. : -1. ) ;
339 m_deltaCF_3 = acos( zCF_3 / 2. ) * ( m_wCFSign_3 ? 1. : -1. ) ;
340 m_deltaCS = acos( zCS / 2. ) * ( m_wCSSign ? 1. : -1. ) ;
341
342 double rCF = m_rCF ;
343 double zCF = m_zCF ;
344 double rCF2 = rCF * rCF ;
345 double zCF2 = zCF * zCF ;
346 double wCF = ( m_wCFSign ? 1. : -1. ) * sqrt( 4. - zCF2 ) ;
347 double vCFCSPlus = ( zCF * zCS + wCF * wCS ) / 2. ;
348 double vCFCSMinus = ( zCF * zCS - wCF * wCS ) / 2. ;
349 m_deltaCF = acos( zCF / 2. ) * ( m_wCFSign ? 1. : -1. ) ;
350
351 double bCF_0 = 1. + 0.5 * rCF_0 * ( zCF_0 * y + wCF_0 * x ) + 0.5 * ( y * y - x * x );//second order y term added -DA
352 double bCF_1 = 1. + 0.5 * rCF_1 * ( zCF_1 * y + wCF_1 * x ) + 0.5 * ( y * y - x * x );//second order y term added -DA
353 double bCF_3 = 1. + 0.5 * rCF_3 * ( zCF_3 * y + wCF_3 * x ) + 0.5 * ( y * y - x * x );//second order y term added -DA
354 double bCS = 1. + 0.5 * rCS * ( zCS * y + wCS * x ) + 0.5 * ( y * y - x * x );
355 double bCPPlus = 1. - y ;
356 double bCPMinus = 1. + y ;
357
358 if( !m_useNewWeights )
359 {
360 bCF_0 = 1. ;
361 bCF_1 = 1. ;
362 bCF_3 = 1. ;
363 m_rwsCF_0 = rCF2_0 ;
364 m_rwsCF_1 = rCF2_1 ;
365 m_rwsCF_3 = rCF2_3 ;
366 bCS = 1. ;
367 m_rwsCS = rCS2 ;
368 bCPPlus = 1. ;
369 bCPMinus = 1. ;
370 }
371
372 m_rwsCF_0 = ( rCF2_0 + 0.5 * rCF_0 * ( zCF_0 * y - wCF_0 * x ) + rmix ) / bCF_0;//though division by b not mentioned in memo it is correct
373 m_rwsCF_1 = ( rCF2_1 + 0.5 * rCF_1 * ( zCF_1 * y - wCF_1 * x ) + rmix ) / bCF_1;//though division by b not mentioned in memo it is correct
374 m_rwsCF_3 = ( rCF2_3 + 0.5 * rCF_3 * ( zCF_3 * y - wCF_3 * x ) + rmix ) / bCF_3;//though division by b not mentioned in memo it is correct
375 m_rwsCS = ( rCS2 + 0.5 * rCS * ( zCS * y - wCS * x ) + rmix ) / bCS;
376
377 double bfCF_0 = ( 1. - RCF_0 * rCF_0 * ( 0.5*zCF_0*y + 0.5*wCF_0*x ) + 0.5 * ( -1*x*x + y*y) );//B(D0->F)=AF*AF*brCF_0
378 double bfCF_1 = ( 1. - RCF_1 * rCF_1 * ( 0.5*zCF_1*y + 0.5*wCF_1*x ) + 0.5 * ( -1*x*x + y*y) );
379 double bfCF_3 = ( 1. - RCF_3 * rCF_3 * ( 0.5*zCF_3*y + 0.5*wCF_3*x ) + 0.5 * ( -1*x*x + y*y) );
380 double bfCFbar_0 = ( rCF2_0 - RCF_0 * rCF_0 * ( 0.5*zCF_0*y - 0.5*wCF_0*x ) + 0.5 * ( x*x + y*y) );//B(D0->Fbar)=AF*AF*brCFbar_0
381 double bfCFbar_1 = ( rCF2_1 - RCF_1 * rCF_1 * ( 0.5*zCF_1*y - 0.5*wCF_1*x ) + 0.5 * ( x*x + y*y) );
382 double bfCFbar_3 = ( rCF2_3 - RCF_3 * rCF_3 * ( 0.5*zCF_3*y - 0.5*wCF_3*x ) + 0.5 * ( x*x + y*y) );
383
384
385 // Calculate reweighting factors, normalized to the largest factor.
386 // Note:A fake Dalitz weight is added in to just to make initial predictions
387 // These Dalitz weights are not used in actual code
388
389 // FlavoredCF
390
391 // First calculate rate factors.
392 // m_weights[ kFlavoredCF_0 ][ kFlavoredCF_0 ] = ( ( 1 - RCF_0 * RCF_0 ) / ( 1 - RCF_0 * (1./rCF_0) * (0.5*zCF_0*y - 0.5*wCF_0*x) + rmix/(2*rCF2_0) ) ) ;
393 // rmix * m_weights[ kFlavoredCF ][ kFlavoredCFBar ] ;
394 // m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_0 ] = ( 1. + rCF2_0 * rCF2_0 - 2. * RCF_0 * RCF_0 * rCF2_0 ) / ( brCF_0 * brCF_0 );
395 // 1. + rCF2 * ( 2. - zCF2 ) + rCF2 * rCF2 ;
396
397 // Nov 2007: correct for r2 -> RWS and BR != A^2
398 // m_weights[ kFlavoredCF ][ kFlavoredCF ] /= m_rwsCF * bCF * bCF ;
399 // m_weights[ kFlavoredCF ][ kFlavoredCFBar ] /=
400 // ( 1. + m_rwsCF * m_rwsCF ) * bCF * bCF ;
401 // m_weights[ kFlavoredCF ][ kFlavoredCS ] /=
402 // ( m_rwsCF + m_rwsCS ) * bCF * bCS ;
403 // m_weights[ kFlavoredCF ][ kFlavoredCSBar ] /=
404 // ( 1. + m_rwsCF * m_rwsCS ) * bCF * bCS ;
405 m_weights[ kFlavoredCF_0 ][ kFlavoredCF_0 ] = ( ( 1 - RCF_0 * RCF_0 ) / ( 1 - RCF_0 * (1./rCF_0) * (0.5*zCF_0*y - 0.5*wCF_0*x) + rmix/(2*rCF2_0) ) ) ;
406 m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_0 ] = ( 1. + rCF2_0 * rCF2_0 - 2. * RCF_0 * RCF_0 * rCF2_0 ) / ( bfCF_0 * bfCF_0 );
407 m_weights[ kFlavoredCF_0 ][ kFlavoredCF_1 ] = ( 1. + ( rCF2_0 / rCF2_1 ) * ( rCF2_0 / rCF2_1 ) - 2. * ( rCF2_0 / rCF2_1 ) * RCF_0 * RCF_1 * cos( dCF_0 - dCF_1 ) ) / ( 1. + ( rCF2_0 / rCF2_1 ) * ( rCF2_0 / rCF2_1 ) - RCF_1 * (1./rCF_1) * ( 0.5*zCF_1*y - 0.5*wCF_1*x ) - (RCF_0*rCF_1/(rCF_0*rCF_0)) * ( 0.5*zCF_0*y - 0.5*wCF_0*x ) + rmix/(rCF_0*rCF_0) );
408 cout<<( 1. + ( rCF2_0 / rCF2_1 ) * ( rCF2_0 / rCF2_1 ) - RCF_1 * (1./rCF_1) * ( 0.5*zCF_1*y - 0.5*wCF_1*x ) - (RCF_0*rCF_1/(rCF_0*rCF_0)) * ( 0.5*zCF_0*y - 0.5*wCF_0*x ) + rmix/(rCF_0*rCF_0) )<<endl;
409 m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_1 ] = ( 1. + rCF2_0 * rCF2_1 - 2. * RCF_0 * RCF_1 * rCF_0 * rCF_1 ) / ( bfCF_0 * bfCF_1 );
410 m_weights[ kFlavoredCF_0 ][ kFlavoredCF_3 ] = ( 1. + ( rCF2_0 / rCF2_3 ) * ( rCF2_0 / rCF2_3 ) - 2. * ( rCF2_0 / rCF2_3 ) * RCF_0 * RCF_3 * cos( dCF_0 - dCF_3 ) ) / ( 1. + ( rCF2_0 / rCF2_3 ) * ( rCF2_0 / rCF2_3 ) - RCF_3 * (1./rCF_3) * ( 0.5*zCF_3*y - 0.5*wCF_3*x ) - (RCF_0/(rCF_0*rCF_0)) * ( 0.5*zCF_0*y - 0.5*wCF_0*x ) + rmix/(rCF_0*rCF_0) );
411 m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_3 ] = ( 1. + rCF2_0 * rCF2_3 - 2. * RCF_0 * RCF_3 * rCF_0 * rCF_3 ) / ( bfCF_0 * bfCF_3 );
412 m_weights[ kFlavoredCF_0 ][ kFlavoredCSBar ] =
413 1. + rCF2_0 * rCS2 - rCF_0 * rCS * vCFCSMinus ;
414 m_weights[ kFlavoredCF_0 ][ kFlavoredCS ] =
415 rCF2_0 + rCS2 - rCF_0 * rCS * vCF_0CSPlus +
416 rmix * m_weights[ kFlavoredCF_0 ][ kFlavoredCSBar ] ;
417
418 m_weights[ kFlavoredCF_0 ][ kSLPlus ] = 1. / bfCF_0 ; //inexist
419 m_weights[ kFlavoredCF_0 ][ kSLMinus ] = 1. / bfCF_0 ; //
420
421 //m_weights[ kFlavoredCF_0 ][ kCPPlus ] =
422 // ( 1. + rCF2_0 + rCF_0 * zCF_0 ) / ( 1. + m_rwsCF_0 ) / bCF_0 / bCPPlus ;
423 //m_weights[ kFlavoredCF_0 ][ kCPMinus ] =
424 // ( 1. + rCF2_0 - rCF_0 * zCF_0 ) / ( 1. + m_rwsCF_0 ) / bCF_0 / bCPMinus ;
425 m_weights[ kFlavoredCF_0 ][ kCPPlus ] =
426 ( 1. + rCF2_0 - 2. * rCF_0 * RCF_0 * cos(dCF_0) ) / ( ( 1. + rCF2_0 - 2. * rCF_0 * RCF_0 * cos(dCF_0) * y + y * y ) * bCPPlus );
427 m_weights[ kFlavoredCF_0 ][ kCPMinus ] =
428 ( 1. + rCF2_0 + 2. * rCF_0 * RCF_0 * cos(dCF_0) ) / ( ( 1. + rCF2_0 - 2. * rCF_0 * RCF_0 * cos(dCF_0) * y + y * y ) * bCPMinus );
429 m_weights[ kFlavoredCF_0 ][ kDalitz ] = 1. / bfCF_0;
430 m_weights[ kFlavoredCF_0 ][ k2Pip2Pim ] = 1. / bfCF_0;
431 m_weights[ kFlavoredCF_0 ][ kPipPim2Pi0 ] = 1. / bfCF_0;
432 m_weights[ kFlavoredCF_0 ][ kK0PiPiPi0 ] = 1. / bfCF_0;
433 m_weights[ kFlavoredCF_0 ][ kKKPiPi ] = 1. / bfCF_0;
434
435 // FlavoredCFBar
436 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCFBar_0 ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCF_0 ] ;
437 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCF_1 ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_1 ] ;
438 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCFBar_1 ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCF_1 ] ;
439 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCF_3 ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCFBar_3 ] ;
440 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCFBar_3 ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCF_3 ] ;
441 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCS ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCSBar ] ;
442 m_weights[ kFlavoredCFBar_0 ][ kFlavoredCSBar ] = m_weights[ kFlavoredCF_0 ][ kFlavoredCS ] ;
443 m_weights[ kFlavoredCFBar_0 ][ kSLPlus ] = 1. / bCF_0 ;
444 m_weights[ kFlavoredCFBar_0 ][ kSLMinus ] = 1. / bCF_0 ;
445 m_weights[ kFlavoredCFBar_0 ][ kCPPlus ] = m_weights[ kFlavoredCF_0 ][ kCPPlus ] ;
446 m_weights[ kFlavoredCFBar_0 ][ kCPMinus ] = m_weights[ kFlavoredCF_0 ][ kCPMinus ] ;
447 m_weights[ kFlavoredCFBar_0 ][ kDalitz ] = 1. / bCF_0;
448 m_weights[ kFlavoredCFBar_0 ][ k2Pip2Pim ] = 1. / bCF_0;
449 m_weights[ kFlavoredCFBar_0 ][ kPipPim2Pi0 ] = 1. / bCF_0;
450 m_weights[ kFlavoredCFBar_0 ][ kK0PiPiPi0 ] = 1. / bCF_0;
451 m_weights[ kFlavoredCFBar_0 ][ kKKPiPi ] = 1. / bCF_0;
452
453 m_weights[ kFlavoredCF_1 ][ kFlavoredCF_1 ] = ( ( 1 - RCF_1 * RCF_1 ) / ( 1 - RCF_1 * (1./rCF_1) * (0.5*zCF_1*y - 0.5*wCF_1*x) + rmix/(2*rCF2_1) ) ) ;
454 m_weights[ kFlavoredCF_1 ][ kFlavoredCFBar_1 ] = ( 1. + rCF2_1 * rCF2_1 - 2. * RCF_1 * RCF_1 * rCF2_1 ) / ( bfCF_1 * bfCF_1 );
455 m_weights[ kFlavoredCF_1 ][ kFlavoredCF_3 ] = ( 1. + ( rCF2_1 / rCF2_3 ) * ( rCF2_1 / rCF2_3 ) - 2. * ( rCF2_1 / rCF2_3 ) * RCF_1 * RCF_3 * cos( dCF_1 - dCF_3 ) ) / ( 1. + ( rCF2_1 / rCF2_3 ) * ( rCF2_1 / rCF2_3 ) - RCF_3 * (1./rCF_3) * ( 0.5*zCF_3*y - 0.5*wCF_3*x ) - (RCF_1/(rCF_1*rCF_1)) * ( 0.5*zCF_1*y - 0.5*wCF_1*x ) + rmix/(rCF_1*rCF_1) );;
456 m_weights[ kFlavoredCF_1 ][ kFlavoredCFBar_3 ] = ( 1. + rCF2_1 * rCF2_3 - 2. * RCF_1 * RCF_3 * rCF_1 * rCF_3 ) / ( bfCF_1 * bfCF_3 );;
457 m_weights[ kFlavoredCF_1 ][ kFlavoredCSBar ] =
458 1. + rCF2_1 * rCS2 - rCF_1 * rCS * vCFCSMinus ;
459 m_weights[ kFlavoredCF_1 ][ kFlavoredCS ] =
460 rCF2_1 + rCS2 - rCF_1 * rCS * vCF_1CSPlus +
461 rmix * m_weights[ kFlavoredCF_1 ][ kFlavoredCSBar ] ;
462
463 m_weights[ kFlavoredCF_1 ][ kSLPlus ] = 1. / bfCF_1 ;
464 m_weights[ kFlavoredCF_1 ][ kSLMinus ] = 1. / bfCF_1 ;
465
466 // m_weights[ kFlavoredCF_1 ][ kCPPlus ] =
467 // ( 1. + rCF2_1 + rCF_1 * zCF_1 ) / ( 1. + m_rwsCF_1 ) / bCF_1 / bCPPlus ;
468 // m_weights[ kFlavoredCF_1 ][ kCPMinus ] =
469 // ( 1. + rCF2_1 - rCF_1 * zCF_1 ) / ( 1. + m_rwsCF_1 ) / bCF_1 / bCPMinus ;
470 m_weights[ kFlavoredCF_1 ][ kCPPlus ] =
471 ( 1. + rCF2_1 - 2. * rCF_1 * RCF_1 * cos(dCF_1) ) / ( ( 1. + rCF2_1 - 2. * rCF_1 * RCF_1 * cos(dCF_1) *y + y * y ) * bCPPlus );
472 m_weights[ kFlavoredCF_1 ][ kCPMinus ] =
473 ( 1. + rCF2_1 + 2. * rCF_1 * RCF_1 * cos(dCF_1) ) / ( ( 1. + rCF2_1 - 2. * rCF_1 * RCF_1 * cos(dCF_1) *y + y * y ) * bCPMinus );
474 m_weights[ kFlavoredCF_1 ][ kDalitz ] = 1. / bCF_1;
475 m_weights[ kFlavoredCF_1 ][ k2Pip2Pim ] = 1. / bCF_1;
476 m_weights[ kFlavoredCF_1 ][ kPipPim2Pi0 ] = 1. / bCF_1;
477 m_weights[ kFlavoredCF_1 ][ kK0PiPiPi0 ] = 1. / bCF_1;
478 m_weights[ kFlavoredCF_1 ][ kKKPiPi ] = 1. / bCF_1;
479
480 // FlavoredCFBar
481 m_weights[ kFlavoredCFBar_1 ][ kFlavoredCFBar_1 ] = m_weights[ kFlavoredCF_1 ][ kFlavoredCF_1 ] ;
482 m_weights[ kFlavoredCFBar_1 ][ kFlavoredCF_3 ] = m_weights[ kFlavoredCF_1 ][ kFlavoredCFBar_3 ] ;
483 m_weights[ kFlavoredCFBar_1 ][ kFlavoredCFBar_3 ] = m_weights[ kFlavoredCF_1 ][ kFlavoredCF_3 ] ;
484 m_weights[ kFlavoredCFBar_1 ][ kFlavoredCS ] = m_weights[ kFlavoredCF_1 ][ kFlavoredCSBar ] ;
485 m_weights[ kFlavoredCFBar_1 ][ kFlavoredCSBar ] = m_weights[ kFlavoredCF_1 ][ kFlavoredCS ] ;
486 m_weights[ kFlavoredCFBar_1 ][ kSLPlus ] = 1. / bCF_1 ;
487 m_weights[ kFlavoredCFBar_1 ][ kSLMinus ] = 1. / bCF_1 ;
488 m_weights[ kFlavoredCFBar_1 ][ kCPPlus ] = m_weights[ kFlavoredCF_1 ][ kCPPlus ] ;
489 m_weights[ kFlavoredCFBar_1 ][ kCPMinus ] = m_weights[ kFlavoredCF_1 ][ kCPMinus ] ;
490 m_weights[ kFlavoredCFBar_1 ][ kDalitz ] = 1. / bCF_1;
491 m_weights[ kFlavoredCFBar_1 ][ k2Pip2Pim ] = 1. / bCF_1;
492 m_weights[ kFlavoredCFBar_1 ][ kPipPim2Pi0 ] = 1. / bCF_1;
493 m_weights[ kFlavoredCFBar_1 ][ kK0PiPiPi0 ] = 1. / bCF_1;
494 m_weights[ kFlavoredCFBar_1 ][ kKKPiPi ] = 1. / bCF_1;
495
496 m_weights[ kFlavoredCF_3 ][ kFlavoredCF_3 ] = ( ( 1 - RCF_3 * RCF_3 ) / ( 1 - RCF_3 * (1./rCF_3) * (0.5*zCF_3*y - 0.5*wCF_3*x) + rmix/(2*rCF2_3) ) ) ;
497 m_weights[ kFlavoredCF_3 ][ kFlavoredCFBar_3 ] = ( 1. + rCF2_3 * rCF2_3 - 2. * RCF_3 * RCF_3 * rCF2_3 ) / ( bfCF_3 * bfCF_3 );
498 m_weights[ kFlavoredCF_3 ][ kFlavoredCSBar ] =
499 1. + rCF2_3 * rCS2 - rCF_3 * rCS * vCFCSMinus ;
500 m_weights[ kFlavoredCF_3 ][ kFlavoredCS ] =
501 rCF2_3 + rCS2 - rCF_3 * rCS * vCF_3CSPlus +
502 rmix * m_weights[ kFlavoredCF_3 ][ kFlavoredCSBar ] ;
503
504 m_weights[ kFlavoredCF_3 ][ kSLPlus ] = 1. / bfCF_3 ;
505 m_weights[ kFlavoredCF_3 ][ kSLMinus ] = 1. / bfCF_3 ;
506
507 //m_weights[ kFlavoredCF_3 ][ kCPPlus ] =
508 // ( 1. + rCF2_3 + rCF_3 * zCF_3 ) / ( 1. + m_rwsCF_3 ) / bCF_3 / bCPPlus ;
509 //m_weights[ kFlavoredCF_3 ][ kCPMinus ] =
510 // ( 1. + rCF2_3 - rCF_3 * zCF_3 ) / ( 1. + m_rwsCF_3 ) / bCF_3 / bCPMinus ;
511 m_weights[ kFlavoredCF_3 ][ kCPPlus ] =
512 ( 1. + rCF2_3 - 2. * rCF_3 * RCF_3 * cos(dCF_3) ) / ( ( 1. + rCF2_3 - 2. * rCF_3 * RCF_3 * cos(dCF_3) * y + y * y ) * bCPPlus );
513 m_weights[ kFlavoredCF_3 ][ kCPMinus ] =
514 ( 1. + rCF2_3 + 2. * rCF_3 * RCF_3 * cos(dCF_3) ) / ( ( 1. + rCF2_3 - 2. * rCF_3 * RCF_3 * cos(dCF_3) * y + y * y ) * bCPMinus );
515 m_weights[ kFlavoredCF_3 ][ kDalitz ] = 1. / bCF_3;
516 m_weights[ kFlavoredCF_3 ][ k2Pip2Pim ] = 1. / bCF_3;
517 m_weights[ kFlavoredCF_3 ][ kPipPim2Pi0 ] = 1. / bCF_3;
518 m_weights[ kFlavoredCF_3 ][ kK0PiPiPi0 ] = 1. / bCF_3;
519 m_weights[ kFlavoredCF_3 ][ kKKPiPi ] = 1. / bCF_3;
520
521 // FlavoredCFBar
522 m_weights[ kFlavoredCFBar_3 ][ kFlavoredCFBar_3 ] = m_weights[ kFlavoredCF_3 ][ kFlavoredCF_3 ] ;
523 m_weights[ kFlavoredCFBar_3 ][ kFlavoredCS ] = m_weights[ kFlavoredCF_3 ][ kFlavoredCSBar ] ;
524 m_weights[ kFlavoredCFBar_3 ][ kFlavoredCSBar ] = m_weights[ kFlavoredCF_3 ][ kFlavoredCS ] ;
525 m_weights[ kFlavoredCFBar_3 ][ kSLPlus ] = 1. / bCF_3 ;
526 m_weights[ kFlavoredCFBar_3 ][ kSLMinus ] = 1. / bCF_3 ;
527 m_weights[ kFlavoredCFBar_3 ][ kCPPlus ] = m_weights[ kFlavoredCF_3 ][ kCPPlus ] ;
528 m_weights[ kFlavoredCFBar_3 ][ kCPMinus ] = m_weights[ kFlavoredCF_3 ][ kCPMinus ] ;
529 m_weights[ kFlavoredCFBar_3 ][ kDalitz ] = 1. / bCF_3;
530 m_weights[ kFlavoredCFBar_3 ][ k2Pip2Pim ] = 1. / bCF_3;
531 m_weights[ kFlavoredCFBar_3 ][ kPipPim2Pi0 ] = 1. / bCF_3;
532 m_weights[ kFlavoredCFBar_3 ][ kK0PiPiPi0 ] = 1. / bCF_3;
533 m_weights[ kFlavoredCFBar_3 ][ kKKPiPi ] = 1. / bCF_3;
534
535 // FlavoredCS
536
537 // First calculate rate factors.
538 m_weights[ kFlavoredCS ][ kFlavoredCSBar ] = 1. + rCS2 * ( 2. - zCS2 ) + rCS2 * rCS2 ;
539 m_weights[ kFlavoredCS ][ kFlavoredCS ] = rmix * m_weights[ kFlavoredCS ][ kFlavoredCSBar ] ;
540
541 // Nov 2007: correct for r2 -> RWS and BR != A^2
542 m_weights[ kFlavoredCS ][ kFlavoredCS ] /= m_rwsCS * bCS * bCS ;
543 m_weights[ kFlavoredCS ][ kFlavoredCSBar ] /= ( 1. + m_rwsCS * m_rwsCS ) * bCS * bCS ;
544
545 m_weights[ kFlavoredCS ][ kSLPlus ] = 1. / bCS ;
546 m_weights[ kFlavoredCS ][ kSLMinus ] = 1. / bCS ;
547
548 m_weights[ kFlavoredCS ][ kCPPlus ] = ( 1. + rCS2 + rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPPlus ;
549 m_weights[ kFlavoredCS ][ kCPMinus ] = ( 1. + rCS2 - rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPMinus ;
550
551 m_weights[ kFlavoredCS ][ kDalitz ] = 1. / bCS;
552 m_weights[ kFlavoredCS ][ k2Pip2Pim ] = 1. / bCS;
553 m_weights[ kFlavoredCS ][ kPipPim2Pi0 ] = 1. / bCS;
554 m_weights[ kFlavoredCS ][ kK0PiPiPi0 ] = 1. / bCS;
555 m_weights[ kFlavoredCS ][ kKKPiPi ] = 1. / bCS;
556
557 // FlavoredCSBar
558
559 m_weights[ kFlavoredCSBar ][ kFlavoredCSBar ] = m_weights[ kFlavoredCS ][ kFlavoredCS ] ;
560 m_weights[ kFlavoredCSBar ][ kSLPlus ] = 1. / bCS ;
561 m_weights[ kFlavoredCSBar ][ kSLMinus ] = 1. / bCS ;
562 m_weights[ kFlavoredCSBar ][ kCPPlus ] = m_weights[ kFlavoredCS ][ kCPPlus ] ;
563 m_weights[ kFlavoredCSBar ][ kCPMinus ] = m_weights[ kFlavoredCS ][ kCPMinus ] ;
564 m_weights[ kFlavoredCSBar ][ kDalitz ] = 1. / bCS;
565 m_weights[ kFlavoredCSBar ][ k2Pip2Pim ] = 1. / bCS;
566 m_weights[ kFlavoredCSBar ][ kPipPim2Pi0 ] = 1. / bCS;
567 m_weights[ kFlavoredCSBar ][ kK0PiPiPi0 ] = 1. / bCS;
568 m_weights[ kFlavoredCSBar ][ kKKPiPi ] = 1. / bCS;
569
570 // SLPlus
571
572 // SLPlus/SLPlus should have rate factor = Rmix, but none of these are
573 // generated in uncorrelated MC.
574 m_weights[ kSLPlus ][ kSLPlus ] = 0. ;
575 m_weights[ kSLPlus ][ kSLMinus ] = 1. ;
576 m_weights[ kSLPlus ][ kCPPlus ] = 1. / bCPPlus ;
577 m_weights[ kSLPlus ][ kCPMinus ] = 1. / bCPMinus ;
578 m_weights[ kSLPlus ][ kDalitz ] = 1. ;
579 m_weights[ kSLPlus ][ k2Pip2Pim ] = 1. ;
580 m_weights[ kSLPlus ][ kPipPim2Pi0 ] = 1. ;
581 m_weights[ kSLPlus ][ kK0PiPiPi0 ] = 1. ;
582 m_weights[ kSLPlus ][ kKKPiPi ] = 1. ;
583
584 // SLMinus
585
586 m_weights[ kSLMinus ][ kSLMinus ] = 0. ;
587 m_weights[ kSLMinus ][ kCPPlus ] = 1. / bCPPlus ;
588 m_weights[ kSLMinus ][ kCPMinus ] = 1. / bCPMinus ;
589 m_weights[ kSLMinus ][ kDalitz ] = 1. ;
590 m_weights[ kSLMinus ][ k2Pip2Pim ] = 1. ;
591 m_weights[ kSLMinus ][ kPipPim2Pi0 ] = 1. ;
592 m_weights[ kSLMinus ][ kK0PiPiPi0 ] = 1. ;
593 m_weights[ kSLMinus ][ kKKPiPi ] = 1. ;
594
595 // CPPlus
596
597 m_weights[ kCPPlus ][ kCPPlus ] = 0. ;
598 m_weights[ kCPPlus ][ kCPMinus ] = 2. / bCPPlus / bCPMinus ;
599 m_weights[ kCPPlus ][ kDalitz ] = 1. / bCPPlus;
600 m_weights[ kCPPlus ][ k2Pip2Pim ] = 1. / bCPPlus;
601 m_weights[ kCPPlus ][ kPipPim2Pi0 ] = 1. / bCPPlus;
602 m_weights[ kCPPlus ][ kK0PiPiPi0 ] = 1. / bCPPlus;
603 m_weights[ kCPPlus ][ kKKPiPi ] = 1. / bCPPlus;
604
605 // CPMinus
606
607 m_weights[ kCPMinus ][ kCPMinus ] = 0. ;
608 m_weights[ kCPMinus ][ kDalitz ] = 1. / bCPMinus;
609 m_weights[ kCPMinus ][ k2Pip2Pim ] = 1. / bCPMinus;
610 m_weights[ kCPMinus ][ kPipPim2Pi0 ] = 1. / bCPMinus;
611 m_weights[ kCPMinus ][ kK0PiPiPi0 ] = 1. / bCPMinus;
612 m_weights[ kCPMinus ][ kKKPiPi ] = 1. / bCPMinus;
613
614 //Dalitz estimate
615 m_weights[ kDalitz ][ kDalitz ] = 1;
616 m_weights[ kDalitz ][ k2Pip2Pim ] = 1;
617 m_weights[ kDalitz ][ kPipPim2Pi0 ] = 1;
618 m_weights[ kDalitz ][ kK0PiPiPi0 ] = 1;
619 m_weights[ kDalitz ][ kKKPiPi ] = 1;
620
621 m_weights[ k2Pip2Pim ][ k2Pip2Pim ] = 1;
622 m_weights[ k2Pip2Pim ][ kPipPim2Pi0 ] = 1;
623 m_weights[ k2Pip2Pim ][ kK0PiPiPi0 ] = 1;
624 m_weights[ k2Pip2Pim ][ kKKPiPi ] = 1;
625
626 m_weights[ kPipPim2Pi0 ][ kPipPim2Pi0 ] = 1;
627 m_weights[ kPipPim2Pi0 ][ kK0PiPiPi0 ] = 1;
628 m_weights[ kPipPim2Pi0 ][ kKKPiPi ] = 1;
629
630 m_weights[ kK0PiPiPi0 ][ kK0PiPiPi0 ] = 1;
631 m_weights[ kK0PiPiPi0 ][ kKKPiPi ] = 1;
632
633 m_weights[ kKKPiPi ][ kKKPiPi ] = 1;
634 log << MSG::INFO << "Weight matrix" << m_weights << endmsg ;
635
636 //Boss 6.6.2 MC Branching fractions
637 HepVector fractionsD0( kNDecayTypes+1, 0 ) ;
638 fractionsD0[ kFlavoredCF_0 ] = 0.305 ;
639 fractionsD0[ kFlavoredCFBar_0 ] = 0.0076 ;
640 fractionsD0[ kFlavoredCF_1 ] = 0.144 ;
641 fractionsD0[ kFlavoredCFBar_1 ] = 0.00031 ;
642 fractionsD0[ kFlavoredCF_3 ] = 0.0822 ;
643 fractionsD0[ kFlavoredCFBar_3 ] = 0.00027 ;
644 fractionsD0[ kFlavoredCS ] = 0.031 ;
645 fractionsD0[ kFlavoredCSBar ] = 0.031 ;
646 fractionsD0[ kSLPlus ] = 0.125 ;
647 fractionsD0[ kSLMinus ] = 0. ;
648 fractionsD0[ kCPPlus ] = 0.113 ;
649 fractionsD0[ kCPMinus ] = 0.102 ;
650 fractionsD0[ kDalitz ] = 0.056 ;
651 fractionsD0[ k2Pip2Pim ] = 0.00720 ;
652 fractionsD0[ kPipPim2Pi0 ] = 0.0102 ;
653 fractionsD0[ kK0PiPiPi0 ] = 0.104 ;
654 fractionsD0[ kKKPiPi ] = 0.00273 ;
655 HepVector scales = m_weights * fractionsD0 ;
656 log << MSG::INFO << "Integrated scale factors " <<scales << endmsg ;
657
658
659 HepVector fractionsD0bar( kNDecayTypes+1, 0 ) ;
660 fractionsD0bar[ kFlavoredCF_0 ] = 0.0076 ;
661 fractionsD0bar[ kFlavoredCFBar_0 ] = 0.305 ;
662 fractionsD0bar[ kFlavoredCF_1 ] = 0.00031 ;
663 fractionsD0bar[ kFlavoredCFBar_1 ] = 0.144 ;
664 fractionsD0bar[ kFlavoredCF_3 ] = 0.00027 ;
665 fractionsD0bar[ kFlavoredCFBar_3 ] = 0.0822 ;
666 fractionsD0bar[ kFlavoredCS ] = 0.031 ;
667 fractionsD0bar[ kFlavoredCSBar ] = 0.031 ;
668 fractionsD0bar[ kSLPlus ] = 0. ;
669 fractionsD0bar[ kSLMinus ] = 0.125 ;
670 fractionsD0bar[ kCPPlus ] = 0.113 ;
671 fractionsD0bar[ kCPMinus ] = 0.102 ;
672 fractionsD0bar[ kDalitz ] = 0.056 ;
673 fractionsD0bar[ k2Pip2Pim ] = 0.00720 ;
674 fractionsD0bar[ kPipPim2Pi0 ] = 0.0102 ;
675 fractionsD0bar[ kK0PiPiPi0 ] = 0.104 ;
676 fractionsD0bar[ kKKPiPi ] = 0.00273 ;
677 HepVector weight_norm = fractionsD0bar.T() * scales;
678 log << MSG::INFO << "Overall scale factor " << fractionsD0bar.T() * scales << endmsg ;
679
680
681
682 //Original MC
683 //Branching fraction predictions after weight
684 double brCF_0 = ( fractionsD0bar[ kFlavoredCFBar_0 ] * scales[ kFlavoredCFBar_0 ] ) / weight_norm[0] ;
685 double brCF_1 = ( fractionsD0bar[ kFlavoredCFBar_1 ] * scales[ kFlavoredCFBar_1 ] ) / weight_norm[0] ;
686 double brCF_3 = ( fractionsD0bar[ kFlavoredCFBar_3 ] * scales[ kFlavoredCFBar_3 ] ) / weight_norm[0] ;
687 double brCS = ( fractionsD0bar[ kFlavoredCS ] * scales[ kFlavoredCS ] + fractionsD0bar[ kFlavoredCSBar ] * scales[ kFlavoredCSBar ] ) / weight_norm[0] ;
688 double brDCS_0 = ( fractionsD0bar[ kFlavoredCF_0 ] * scales[ kFlavoredCF_0 ] ) / weight_norm[0] ;
689 double brDCS_1 = ( fractionsD0bar[ kFlavoredCF_1 ] * scales[ kFlavoredCF_1 ] ) / weight_norm[0] ;
690 double brDCS_3 = ( fractionsD0bar[ kFlavoredCF_3 ] * scales[ kFlavoredCF_3 ] ) / weight_norm[0] ;
691 double brCPPlus = ( fractionsD0bar[ kCPPlus ] * scales[ kCPPlus ] ) / weight_norm[0] ;
692 double brCPMinus = ( fractionsD0bar[ kCPMinus ] * scales[ kCPMinus ] ) / weight_norm[0] ;
693
694 //y=0.0 dalitz example used for prediction
695 double dalitz_y_num = -0.000177536;
696 double dalitz_y_den = -0.0150846;
697
698 double numFactCF_0 =
699 -rCF_0 * zCF_0 * ( 1. - 0.5 * rCF_0 * wCF_0 * m_x ) ;
700 double numFactCF_1 =
701 -rCF_1 * zCF_1 * ( 1. - 0.5 * rCF_1 * wCF_1 * m_x ) ;
702 double numFactCF_3 =
703 -rCF_3 * zCF_3 * ( 1. - 0.5 * rCF_3 * wCF_3 * m_x ) ;
704 double numFactCS =
705 -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x ) ;
706 double denFactCF_0 = 0.5 * rCF2_0 * zCF2_0 ;
707 double denFactCF_1 = 0.5 * rCF2_1 * zCF2_1 ;
708 double denFactCF_3 = 0.5 * rCF2_3 * zCF2_3 ;
709 double denFactCS = 0.5 * rCS2 * zCS2 ;
710
711 double num_pre =
712 brCPPlus - brCPMinus + brCF_0 * numFactCF_0 + brCF_1 * numFactCF_1 + brCF_3 * numFactCF_3 + 0.5 * brCS * numFactCS + dalitz_y_num;
713 double den_pre =
714 1. - brCPPlus - brCPMinus - brCF_0 * denFactCF_0 - brCF_1 * denFactCF_1 - brCF_3 * denFactCF_3 - 0.5 * brCS * denFactCS + dalitz_y_den ;
715 double y_pre = num_pre / den_pre ;
716 double num =
717 fractionsD0bar[ kCPPlus ] - fractionsD0bar[ kCPMinus ] + fractionsD0bar[ kFlavoredCFBar_0 ] * numFactCF_0 + fractionsD0bar[ kFlavoredCFBar_1 ] * numFactCF_1 + fractionsD0bar[ kFlavoredCFBar_3 ] * numFactCF_3 + fractionsD0bar[ kFlavoredCS ] * numFactCS + dalitz_y_num;
718 double den =
719 1. - fractionsD0bar[ kCPPlus ] - fractionsD0bar[ kCPMinus ] - fractionsD0bar[ kFlavoredCFBar_0 ] * denFactCF_0 - fractionsD0bar[ kFlavoredCFBar_1 ] * denFactCF_1 - fractionsD0bar[ kFlavoredCFBar_3 ] * denFactCF_3 - fractionsD0bar[ kFlavoredCS ] * denFactCS + dalitz_y_den;
720 double y_prematrix = num / den ;
721 double y_test1 = 0.25 * ( ( ( scales[ kCPMinus ] * m_weights[ kCPPlus ][ kSLPlus ] ) / ( scales[ kCPPlus ] * m_weights[ kCPMinus ][ kSLPlus ] ) ) -
722 ( ( scales[ kCPPlus ] * m_weights[ kCPMinus ][ kSLPlus ] ) / ( scales[ kCPMinus ] * m_weights[ kCPPlus ][ kSLPlus ] ) ) ) ;
723
724 log << MSG::INFO << "An Input Y of " << m_y << endmsg ;
725 log << MSG::INFO << "Parent MC has an intrinsic value of y as " << y_prematrix << endmsg ;
726 log << MSG::INFO << "The BR method predics a y of " << y_pre << endmsg ;
727 log << MSG::INFO << "The CP-SL double tag method predicts y of " <<y_test1 << endmsg ;
728
729 // Renormalize
730 m_largestWeight = 0. ;
731 for( int i = 0 ; i < kNDecayTypes ; ++i )
732 {
733 for( int j = 0 ; j < kNDecayTypes ; ++j )
734 {
735 if( m_weights[ i ][ j ] < 0 ) m_weights[ i ][ j ] = 0; //gets rid of negative weights in early calculations
736 if( m_weights[ i ][ j ] > m_largestWeight )
737 {
738 m_largestWeight = m_weights[ i ][ j ] ;
739 }
740 }
741 }
742 m_weights /= m_largestWeight ;
743
744 log << MSG::INFO <<"Renormalized weight matrix " << m_weights << endmsg ;
745
746
747 //Set up weight parameters
748 Dalitz t(8);
749 double D0Sum[ 10 ], D0barSum[ 10 ] ; // last index is for sum over bins
750 TComplex D0D0barSum[ 10 ] ;
751 int points[ 10 ] ;
752
753 for( int i = 0 ; i < 10 ; ++i )
754 {
755 D0Sum[ i ] = 0. ;
756 D0barSum[ i ] = 0. ;
757 D0D0barSum[ i ] = TComplex( 0., 0. ) ;
758 points[ i ] = 0 ;
759 }
760
761 double m2min = 0. ;
762 double m2max = 3. ;
763 int nsteps = 1000 ;
764 double stepsize = ( m2max - m2min ) / ( double ) nsteps ;
765
766 for( int i = 0 ; i < nsteps ; ++i )
767 {
768 double m2i = m2min + ( ( double ) i + 0.5 ) * stepsize ;
769
770 for( int j = 0 ; j < nsteps ; ++j )
771 {
772 double m2j = m2min + ( ( double ) j + 0.5 ) * stepsize ;
773
774 if( t.Point_on_DP( m2i, m2j ) )
775 {
776 double m2k = 1.86450*1.86450 + 0.497648*0.497648 +
777 0.139570*0.139570 + 0.139570*0.139570 - m2i - m2j ;
778
779 TComplex D0, D0bar;
780 if (m_dalitzModel == 0) { //Cleo model
781 D0 = t.CLEO_Amplitude( m2i, m2j, m2k ) ;
782 D0bar = t.CLEO_Amplitude( m2j, m2i, m2k ) ;}
783 if (m_dalitzModel == 1) { //Babar model
784 D0 = t.Babar_Amplitude( m2i, m2j, m2k ) ;
785 D0bar = t.Babar_Amplitude( m2j, m2i, m2k ) ;}
786 if (m_dalitzModel == 2) { //Belle model
787 D0 = t.Amplitude( m2i, m2j, m2k ) ;
788 D0bar = t.Amplitude( m2j, m2i, m2k ) ;}
789
790 if ( j <= i ) {// lower half only
791 int bin = t.getBin( m2i, m2j, m2k ) ;
792 if( bin == -1 ) bin = 8 ;
793
794 ++points[ bin ] ;
795 D0Sum[ bin ] += norm( D0 ) ;
796 D0barSum[ bin ] += norm( D0bar ) ;
797 D0D0barSum[ bin ] += D0 * conj( D0bar ) ;
798
799 ++points[ 9 ] ;
800 D0Sum[ 9 ] += norm( D0 ) ;
801 D0barSum[ 9 ] += norm( D0bar ) ;
802 D0D0barSum[ 9 ] += D0 * conj( D0bar ) ;
803 }
804
805 }
806 }
807 }
808
809 for( int i = 0 ; i < 10 ; ++i )
810 {
811 double r2 = D0barSum[ i ] / D0Sum[ i ] ;
812 double c = real( D0D0barSum[ i ] ) / sqrt( D0barSum[ i ] ) / sqrt( D0Sum[ i ] ) ;
813 double s = imag( D0D0barSum[ i ] ) / sqrt( D0barSum[ i ] ) / sqrt( D0Sum[ i ] ) ;
814
815 cout << "BIN " << i << ": "
816 << points[ i ] << " pts "
817 << r2 << " " << c << " " << s << " " << c*c+s*s
818 << endmsg ;
819 }
820
821 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
822 return StatusCode::SUCCESS;
823}
824
825//********************************************************************
826StatusCode QCMCFilter::execute() {
827
828 //interface to event data service
829 ISvcLocator* svcLocator = Gaudi::svcLocator();
830 StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
831 if (sc.isFailure())
832 std::cout<<"Could not access EventDataSvc!"<<std::endl;
833
834 //setFilterPassed(false);
835
836 MsgStream log(msgSvc(), name());
837 log<< MSG::INFO << "in execute()" << endmsg;
838
839 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
840
841 if(!eventHeader)
842 {
843 cout<<" eventHeader "<<endl;
844 return StatusCode::FAILURE;
845 }
846 long runNo = eventHeader->runNumber();
847 long event = eventHeader->eventNumber();
848 log << MSG::DEBUG << "run, evtnum = "
849 << runNo << " , "
850 << event << endmsg;
851
852 unsigned int QC_status = 0;
853 //eventHeader->setEventTag(0);
854
855 // Get McParticle collection
856 SmartDataPtr<Event::McParticleCol> mcParticles(eventSvc(), EventModel::MC::McParticleCol);
857
858 //Check if event is Psi(3770)
859 int psipp_daughter = 0;
860 int d0_count = 0;
861 int chd_count = 0;
862 for( Event::McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end(); it++)
863 {
864 int pdg_code = (*it)->particleProperty();
865 if ((*it)->primaryParticle()) continue;
866 Event::McParticle it2 = (*it)->mother();
867 int mother_pdg = 0;
868 //finds if the particle is daughter of Psi(3770)
869 mother_pdg = it2.particleProperty();
870 if (mother_pdg == kPsi3770ID)
871 {
872 psipp_daughter++; //Found daughter
873 if (abs(pdg_code) == kD0ID) d0_count++;
874 if (abs(pdg_code) == kDpID) chd_count++;
875 if (pdg_code == kD0BarID) m_nD0bar++;
876 }
877 }
878
879 if( psipp_daughter == 0 ) //didn't find Psi(3770)
880 {
881 ++m_nUnknownEvents ;
882
883 if( m_invertSelection )
884 {
885 //eventHeader->setEventTag(0);
886 QC_status = 0;
887 eventHeader->setEventTag(QC_status);
888
889 log << MSG::INFO << "QC_status " << QC_status << endmsg;
890 log<< MSG::DEBUG << "Discarding event -- did not find psi(3770). " << endmsg ;
891 return StatusCode::SUCCESS; //discard event
892 }
893 else
894 {
895 // -------- Write to root file
896 //setFilterPassed(true);
897 //eventHeader->setEventTag(1);
898 QC_status = 1;
899 eventHeader->setEventTag(QC_status);
900
901 log << MSG::INFO << "QC_status " << QC_status << endmsg;
902
903 return StatusCode::SUCCESS; //kept event
904 }
905 }
906
907 log<< MSG::DEBUG << "Found psi(3770) -->" << endmsg ;
908
909 // Check that top particle has only two daughters + possible FSR photon
910 if ( psipp_daughter > 3 )
911 {
912 ++m_nUnknownEvents ;
913
914 if( m_invertSelection )
915 {
916 //eventHeader->setEventTag(0);
917 QC_status = 0;
918 eventHeader->setEventTag(QC_status);
919
920 log << MSG::INFO << "QC_status " << QC_status << endmsg;
921 log<< MSG::DEBUG << "Discarding event -- psi(3770) has >3 daughters." << endmsg ;
922 return StatusCode::SUCCESS; //discard event
923 }
924 else
925 {
926 // -------- Write to root file
927 //eventHeader->setEventTag(1);
928 QC_status = 1;
929 eventHeader->setEventTag(QC_status);
930
931 log << MSG::INFO << "QC_status " << QC_status << endmsg;
932 //setFilterPassed(true);
933
934
935 return StatusCode::SUCCESS; //kept event
936 }
937 }
938
939 //Check if D+D- state
940 if (chd_count == 2) { // D+D- event
941 ++m_nDpDmEvents ;
942
943 // D+D- must be rewieghted otherwise too many D+D- events relative to D0D0bar.
944 double random = RandFlat::shoot() ;
945
946 if( ( random < m_dplusDminusWeight && !m_invertSelection ) ||
947 ( random > m_dplusDminusWeight && m_invertSelection ) )
948 {
949 // -------- Write to root file
950 //setFilterPassed(true);
951 //eventHeader->setEventTag(1);
952 QC_status = 1;
953 eventHeader->setEventTag(QC_status);
954
955 log << MSG::INFO << "QC_status " << QC_status << endmsg;
956
957
958 return StatusCode::SUCCESS; //event kept
959 }
960 else
961 {
962 //eventHeader->setEventTag(0);
963 QC_status = 0;
964 eventHeader->setEventTag(QC_status);
965
966 log << MSG::INFO << "QC_status " << QC_status << endmsg;
967 log<< MSG::DEBUG << "Discarding event -- D+D- event failed reweighting." << endmsg ;
968 ++m_nDpDmDiscarded ;
969 return StatusCode::SUCCESS; //event discarded
970 }
971 }
972
973 //Check if D0D0Bar
974 if (d0_count != 2) {
975 if( !m_invertSelection )
976 {
977 //eventHeader->setEventTag(0);
978 QC_status = 0;
979 eventHeader->setEventTag(QC_status);
980
981 log << MSG::INFO << "QC_status " << QC_status << endmsg;
982 log<< MSG::DEBUG << "Discarding event -- non-D0-D0bar event." << endmsg ;
983 return StatusCode::SUCCESS; //discard event
984 }
985 else
986 {
987 // -------- Write to root file
988 //eventHeader->setEventTag(1);
989 QC_status = 1;
990 eventHeader->setEventTag(QC_status);
991
992 log << MSG::INFO << "QC_status " << QC_status << endmsg;
993 //setFilterPassed(true);
994
995
996 return StatusCode::SUCCESS; //kept event
997 }
998 }
999
1000 log<< MSG::INFO << "D0-D0bar event." << endmsg ;
1001 m_nD0D0barEvents++ ;
1002
1003 m_rho_flag = 0;
1004
1005 //Get D0-D0bar modes and info
1006 vector<double> d0_decay = findD0Decay(1);
1007 vector<double> d0bar_decay = findD0Decay(-1);
1008
1009 int d0mode = d0_decay[0];
1010 int d0barmode = d0bar_decay[0];
1011 if(m_debug)cout<<"D0 "<<d0mode<<endl;
1012 if(m_debug)cout<<"D0bar "<<d0barmode<<endl;
1013 m_modeCounter[d0mode][d0barmode]++;
1014
1015
1016
1017 //if any event unidentified remove
1018 if (d0_decay[0] == 18 || d0bar_decay[0] == 18 ) {
1019 if( !m_invertSelection )
1020 {
1021 //eventHeader->setEventTag(0);
1022 QC_status = 0;
1023 eventHeader->setEventTag(QC_status);
1024
1025 log << MSG::INFO << "QC_status " << QC_status << endmsg;
1026 log<< MSG::DEBUG << "Discarding event -- unknown D0-D0bar event." << endmsg ;
1027 return StatusCode::SUCCESS; //discard event
1028 }
1029 else
1030 {
1031 // -------- Write to root file
1032 //eventHeader->setEventTag(1);
1033 QC_status = 1;
1034 eventHeader->setEventTag(QC_status);
1035
1036 log << MSG::INFO << "QC_status " << QC_status << endmsg;
1037 //setFilterPassed(true);
1038
1039
1040
1041 return StatusCode::SUCCESS; //kept event
1042 }
1043 }
1044
1045 // Loop over D's
1046 double r2D0 = d0_decay[1] ;
1047 double deltaD0 = d0_decay[2] ;
1048 double RD0 = d0_decay[3] ;
1049 double FCP = d0_decay[4] ;
1050 double r2D0bar = d0bar_decay[1] ;
1051 double deltaD0bar = d0bar_decay[2] ;
1052 double RD0bar = d0bar_decay[3] ;
1053 double FCPbar = d0bar_decay[4] ;
1054
1055
1056 //Weight based on DDbar combination
1057 double weight ;
1058 if( d0_decay[0] == kDalitz || d0bar_decay[0] == kDalitz || d0_decay[0] == k2Pip2Pim || d0bar_decay[0] == k2Pip2Pim || d0_decay[0] == kPipPim2Pi0 || d0bar_decay[0] == kPipPim2Pi0 || d0_decay[0] == kK0PiPiPi0 || d0bar_decay[0] == kK0PiPiPi0
1059 || d0_decay[0] == kKKPiPi || d0bar_decay[0] == kKKPiPi ||d0bar_decay[0] == kFlavoredCF_3 || d0bar_decay[0] == kFlavoredCFBar_3 )
1060 {
1061 double r2prod = r2D0 * r2D0bar ;
1062 double x = m_x ;
1063 double y = m_y ;
1064
1065 //double bD0 = 1. + sqrt( r2D0 ) * ( cos( deltaD0 ) * y + sin( deltaD0 ) * x ) ;
1066 //double rwsD0 = ( r2D0 + sqrt( r2D0 ) * ( cos( deltaD0 ) * y - sin( deltaD0 ) * x ) ) / bD0 ;
1067 //double bD0bar = 1. + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y + sin( deltaD0bar ) * x ) ;
1068 //double rwsD0bar = ( r2D0bar + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y - sin( deltaD0bar ) * x ) ) / bD0bar ;
1069
1070 //weight = 1. + r2prod - 2. * sqrt( r2prod ) * cos( deltaD0 + deltaD0bar ) ;
1071 //weight /= ( 1. + rwsD0 * rwsD0bar ) * bD0 * bD0bar ;
1072 double bD0 = 1. - (2*FCP - 1.) * RD0 * sqrt( r2D0 ) * ( cos( deltaD0 ) * y + sin( deltaD0 ) * x ) ;
1073 double rwsD0 = ( r2D0 - (2*FCP - 1.) * RD0 * sqrt( r2D0 ) * ( cos( deltaD0 ) * y - sin( deltaD0 ) * x ) ) / bD0 ;
1074 double bD0bar = 1. - (2*FCPbar - 1.) * RD0bar * sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y + sin( deltaD0bar ) * x ) ;
1075 double rwsD0bar = ( r2D0bar - (2*FCPbar - 1.) * RD0bar * sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y - sin( deltaD0bar ) * x ) ) / bD0bar ;
1076
1077 log << MSG::INFO << "bD0 == " << bD0<< endmsg;
1078 log << MSG::INFO << "bD0bar == " << bD0bar<< endmsg;
1079 log << MSG::INFO << "rwsD0 == " << rwsD0<< endmsg;
1080 log << MSG::INFO << "rwsD0bar == " << rwsD0bar<< endmsg;
1081 weight = 1. + r2prod - 2. * RD0 * (2*FCP - 1.) * sqrt( r2prod ) * cos( deltaD0 + deltaD0bar ) ;
1082 log << MSG::INFO << "weight == " << weight<< endmsg;
1083 weight /= ( 1. + rwsD0 * rwsD0bar ) * bD0 * bD0bar ;
1084 log << MSG::INFO << "weight == " << weight<< endmsg;
1085 weight /= m_largestWeight ;
1086 log << MSG::INFO << "weight == " << weight<< endmsg;
1087
1088 //get dalitz graph information before cut
1089 int k = -10;
1090 double m2i= 0;
1091 double m2j= 0;
1092 double m2k= 0;
1093 if( d0_decay[0] == kDalitz) {
1094 k = d0bar_decay[0];
1095 m2i= d0_decay[4];
1096 m2j= d0_decay[5];
1097 m2k= d0_decay[7];}
1098 if( d0bar_decay[0] == kDalitz) {
1099 k = d0_decay[0];
1100 m2i= d0bar_decay[4];
1101 m2j= d0bar_decay[5];
1102 m2k= d0bar_decay[7];}
1103
1104 }
1105 else
1106 {
1107 weight = m_weights[ d0_decay[0] ][ d0bar_decay[0] ] ;
1108 }
1109
1110 double random = RandFlat::shoot() ;
1111 log<< MSG::DEBUG << "type: " << d0_decay[0] << " " << d0bar_decay[0] << ", weight " <<
1112 weight << " <? random " << random << endmsg ;
1113
1114 if( ( random < weight && !m_invertSelection ) || ( random > weight && m_invertSelection ) )
1115 {
1116 // -------- Write to root file
1117 QC_status = 1;
1118 eventHeader->setEventTag(QC_status);
1119
1120 log << MSG::INFO << "QC_status " << QC_status << endmsg;
1121 //eventHeader->setEventTag(1);
1122 //setFilterPassed(true);
1123
1124 //get dalitz graph information after filter
1125 /*int k = -10;
1126 double m2i= 0;
1127 double m2j= 0;
1128 double m2k= 0;
1129 double cosd = cos( deltaD0 ) ;
1130 double sind = sin( deltaD0 ) ;
1131 if( d0_decay[0] == kDalitz) {
1132 k = d0bar_decay[0];
1133 m2i= d0_decay[4];
1134 m2j= d0_decay[5];
1135 m2k= d0_decay[7];
1136 if ( m2j - m2i > 0. ) { //avoids doublecounting the branching fraction
1137 dalitzNumer1_fil += -2. * sqrt( r2D0 ) * cosd ;
1138 dalitzNumer2_fil += 2. * r2D0 * cosd * sind * m_x ;
1139 dalitzDenom_fil += -2. * r2D0 * cosd * cosd ; }
1140 }
1141 if( d0bar_decay[0] == kDalitz) {
1142 k = d0_decay[0];
1143 m2i= d0bar_decay[4];
1144 m2j= d0bar_decay[5];
1145 m2k= d0bar_decay[7];
1146 cosd = cos( deltaD0bar ) ;
1147 sind = sin( deltaD0bar ) ;
1148 if( m2j - m2i < 0. ) { //avoids doublecounting the branching fraction
1149 dalitzNumer1_fil += -2. * sqrt( r2D0bar ) * cosd ;
1150 dalitzNumer2_fil += 2. * r2D0bar * cosd * sind * m_x ;
1151 dalitzDenom_fil += -2. * r2D0bar * cosd * cosd ; }
1152 }
1153
1154 m_keptModeCounter[d0mode][d0barmode]++;*/
1155 return StatusCode::SUCCESS; //event kept
1156 }
1157 else
1158 {
1159 QC_status = 0;
1160 eventHeader->setEventTag(QC_status);
1161
1162 log << MSG::INFO << "QC_status " << QC_status << endmsg;
1163 //eventHeader->setEventTag(0);
1164 log << MSG::DEBUG << "Discarding event -- failed reweighting." << endmsg ;
1165
1166 ++m_nD0D0barDiscarded ;
1167 return StatusCode::SUCCESS; //event discarded
1168 }
1169
1170}
1171
1172
1173//*********************************************************************
1174StatusCode QCMCFilter::finalize() {
1175 MsgStream log(msgSvc(), name());
1176 log << MSG::INFO << "in finalize()" << endmsg;
1177 log << MSG::INFO << "Number of unknown events: " << m_nUnknownEvents << endmsg ;
1178 log << MSG::INFO << "Number of D+D- events seen: " << m_nDpDmEvents << endmsg ;
1179 log << MSG::INFO << "Number of D+D- events discarded: " << m_nDpDmDiscarded << endmsg ;
1180 if( m_nDpDmEvents > 0 )
1181 {
1182 log << MSG::INFO <<"Fraction discarded: " << double( m_nDpDmDiscarded ) / double( m_nDpDmEvents ) << endmsg ;
1183 }
1184
1185 log << MSG::INFO << "Number of D0D0bar events seen: " << m_nD0D0barEvents << " " << m_nD0bar << endmsg ;
1186 log << MSG::INFO << "Number of D0D0bar events discarded: " << m_nD0D0barDiscarded << endmsg ;
1187 /*if( m_nD0D0barEvents > 0 && m_nCPPlus > 0 && m_nCPMinus > 0)
1188 {
1189 log << MSG::INFO << "Fraction discarded: " << double( m_nD0D0barDiscarded ) / double( m_nD0D0barEvents ) << endl <<
1190 "Fraction unknown decays: " << 0.5 * double( m_nUnknownDecays ) / double( m_nD0D0barEvents ) <<endmsg ;
1191
1192 double nD0 = 2. * m_nD0D0barEvents ;
1193 double brSL = m_nSL / nD0 ;
1194 double dBrSL = sqrt( brSL * ( 1. - brSL ) / nD0 ) ;
1195 double brCF_0 = m_nFlavoredCFD0_0 / nD0 ;
1196 double dBrCF_0 = sqrt( brCF_0 * ( 1. - brCF_0 ) / nD0 ) ;
1197 double brCF_1 = m_nFlavoredCFD0_1 / nD0 ;
1198 double dBrCF_1 = sqrt( brCF_1 * ( 1. - brCF_1 ) / nD0 ) ;
1199 double brCF_3 = m_nFlavoredCFD0_3 / nD0 ;
1200 double dBrCF_3 = sqrt( brCF_3 * ( 1. - brCF_3 ) / nD0 ) ;
1201 double brCS = m_nFlavoredCSD0 / nD0 ;
1202 double dBrCS = sqrt( brCS * ( 1. - brCS ) / nD0 ) ;
1203 double brDCS_0 = m_nFlavoredDCSD0_0 / nD0 ;
1204 double dBrDCS_0 = sqrt( brDCS_0 * ( 1. - brDCS_0 ) / nD0 ) ;
1205 double brDCS_1 = m_nFlavoredDCSD0_1 / nD0 ;
1206 double dBrDCS_1 = sqrt( brDCS_1 * ( 1. - brDCS_1 ) / nD0 ) ;
1207 double brDCS_3 = m_nFlavoredDCSD0_3 / nD0 ;
1208 double dBrDCS_3 = sqrt( brDCS_3 * ( 1. - brDCS_3 ) / nD0 ) ;
1209 double brCPPlus = m_nCPPlus / nD0 ;
1210 double dBrCPPlus = sqrt( brCPPlus * ( 1. - brCPPlus ) / nD0 ) ;
1211 double brCPMinus = m_nCPMinus / nD0 ;
1212 double dBrCPMinus = sqrt( brCPMinus * ( 1. - brCPMinus ) / nD0 ) ;
1213 double brDalitz = m_nDalitz / nD0 ;
1214 double dBrDalitz = sqrt( brDalitz * ( 1. - brDalitz ) / nD0 ) ;
1215 double fdBrDalitz = dBrDalitz / brDalitz ;
1216 double dalitzNumer1Norm = m_dalitzNumer1 / nD0 ;
1217 double dDalitzNumer1Norm = dalitzNumer1Norm * fdBrDalitz ;
1218 double dalitzNumer2Norm = m_dalitzNumer2 / nD0 ;
1219 double dDalitzNumer2Norm = dalitzNumer2Norm * fdBrDalitz ;
1220 double dalitzDenomNorm = m_dalitzDenom / nD0 ;
1221 double dDalitzDenomNorm = dalitzDenomNorm * fdBrDalitz ;
1222
1223 double br2Pip2Pim = m_n2Pip2Pim / nD0 ;
1224 double dBr2Pip2Pim = sqrt( br2Pip2Pim * ( 1. - br2Pip2Pim ) / nD0 ) ;
1225 double fdBr2Pip2Pim= dBr2Pip2Pim / br2Pip2Pim ;
1226 double Pip2Pim2Numer1Norm = m_2Pip2PimNumer1 / nD0 ;
1227 double d2Pip2PimNumer1Norm = Pip2Pim2Numer1Norm * fdBr2Pip2Pim ;
1228 double Pip2Pim2Numer2Norm = m_2Pip2PimNumer2 / nD0 ;
1229 double d2Pip2PimNumer2Norm = Pip2Pim2Numer2Norm * fdBr2Pip2Pim ;
1230 double Pip2Pim2DenomNorm = m_2Pip2PimDenom / nD0 ;
1231 double d2Pip2PimDenomNorm = Pip2Pim2DenomNorm * fdBr2Pip2Pim ;
1232
1233 double brPipPim2Pi0 = m_nPipPim2Pi0 / nD0 ;
1234 double dBrPipPim2Pi0 = sqrt( brPipPim2Pi0 * ( 1. - brPipPim2Pi0 ) / nD0 ) ;
1235 double fdBrPipPim2Pi0= dBrPipPim2Pi0 / brPipPim2Pi0 ;
1236 double Pip2Pim2Numer1Norm = m_PipPim2Pi0Numer1 / nD0 ;
1237 double dPipPim2Pi0Numer1Norm = Pip2Pim2Numer1Norm * fdBrPipPim2Pi0 ;
1238 double Pip2Pim2Numer2Norm = m_PipPim2Pi0Numer2 / nD0 ;
1239 double dPipPim2Pi0Numer2Norm = Pip2Pim2Numer2Norm * fdBrPipPim2Pi0 ;
1240 double Pip2Pim2DenomNorm = m_PipPim2Pi0Denom / nD0 ;
1241 double dPipPim2Pi0DenomNorm = Pip2Pim2DenomNorm * fdBrPipPim2Pi0 ;
1242
1243 double brKSPiPiPi0 = m_nKSPiPiPi0 / nD0 ;
1244 double dBrKSPiPiPi0 = sqrt( brKSPiPiPi0 * ( 1. - brKSPiPiPi0 ) / nD0 ) ;
1245 double fdBrKSPiPiPi0= dBrKSPiPiPi0 / brKSPiPiPi0 ;
1246 double Pip2Pim2Numer1Norm = m_KSPiPiPi0Numer1 / nD0 ;
1247 double dKSPiPiPi0Numer1Norm = Pip2Pim2Numer1Norm * fdBrKSPiPiPi0 ;
1248 double Pip2Pim2Numer2Norm = m_KSPiPiPi0Numer2 / nD0 ;
1249 double dKSPiPiPi0Numer2Norm = Pip2Pim2Numer2Norm * fdBrKSPiPiPi0 ;
1250 double Pip2Pim2DenomNorm = m_KSPiPiPi0Denom / nD0 ;
1251 double dKSPiPiPi0DenomNorm = Pip2Pim2DenomNorm * fdBrKSPiPiPi0 ;
1252
1253 //after the filter was applied
1254 double fil_SL = 0, fil_FlavoredCFD0_0 = 0, fil_FlavoredCFD0_1 = 0, fil_FlavoredCFD0_3 = 0, fil_FlavoredCSD0 = 0, fil_FlavoredDCSD0_0 =0, fil_FlavoredDCSD0_1 =0, fil_FlavoredDCSD0_3 =0, fil_CPPlus = 0, fil_CPMinus = 0, fil_Dalitz = 0, fil_2Pip2Pim = 0;
1255 for( int i = 0 ; i < 14 ; ++i )
1256 {
1257 fil_SL += m_keptModeCounter[i][9] + m_keptModeCounter[9][i] + m_keptModeCounter[i][8] + m_keptModeCounter[8][i] ;
1258 fil_FlavoredCFD0_0 += m_keptModeCounter[i][1] + m_keptModeCounter[0][i];
1259 fil_FlavoredCFD0_1 += m_keptModeCounter[i][3] + m_keptModeCounter[2][i];
1260 fil_FlavoredCFD0_3 += m_keptModeCounter[i][5] + m_keptModeCounter[4][i];
1261 fil_FlavoredCSD0 += m_keptModeCounter[i][7] + m_keptModeCounter[7][i] + m_keptModeCounter[i][6] + m_keptModeCounter[6][i];
1262 fil_FlavoredDCSD0_0 += m_keptModeCounter[i][0] + m_keptModeCounter[1][i];
1263 fil_FlavoredDCSD0_1 += m_keptModeCounter[i][2] + m_keptModeCounter[3][i];
1264 fil_FlavoredDCSD0_3 += m_keptModeCounter[i][4] + m_keptModeCounter[5][i];
1265 fil_CPPlus += m_keptModeCounter[i][10] + m_keptModeCounter[10][i];
1266 fil_CPMinus += m_keptModeCounter[i][11] + m_keptModeCounter[11][i];
1267 fil_Dalitz += m_keptModeCounter[i][12] + m_keptModeCounter[12][i];
1268 fil_2Pip2Pim += m_keptModeCounter[i][13] + m_keptModeCounter[13][i];
1269 fil_PipPim2Pi0 += m_keptModeCounter[i][14] + m_keptModeCounter[14][i];
1270 fil_KSPiPiPi0 += m_keptModeCounter[i][15] + m_keptModeCounter[15][i];
1271}
1272double nD0_fil = 2. * ( m_nD0D0barEvents - m_nD0D0barDiscarded ) ;
1273double brSL_fil = fil_SL / nD0_fil ;
1274double dBrSL_fil = sqrt( brSL_fil * ( 1. - brSL_fil ) / nD0_fil ) ;
1275double brCF_0_fil = fil_FlavoredCFD0_0 / nD0_fil ;
1276double dBrCF_0_fil = sqrt( brCF_0_fil * ( 1. - brCF_0_fil ) / nD0_fil ) ;
1277double brCF_1_fil = fil_FlavoredCFD0_1 / nD0_fil ;
1278double dBrCF_1_fil = sqrt( brCF_1_fil * ( 1. - brCF_1_fil ) / nD0_fil ) ;
1279double brCF_3_fil = fil_FlavoredCFD0_3 / nD0_fil ;
1280double dBrCF_3_fil = sqrt( brCF_3_fil * ( 1. - brCF_3_fil ) / nD0_fil ) ;
1281double brCS_fil = fil_FlavoredCSD0 / nD0_fil ;
1282double dBrCS_fil = sqrt( brCS_fil * ( 1. - brCS_fil ) / nD0_fil ) ;
1283double brDCS_0_fil = fil_FlavoredDCSD0_0 / nD0_fil ;
1284double brDCS_1_fil = fil_FlavoredDCSD0_1 / nD0_fil ;
1285double brDCS_3_fil = fil_FlavoredDCSD0_3 / nD0_fil ;
1286double dBrDCS_0_fil = sqrt( brDCS_0_fil * ( 1. - brDCS_0_fil ) / nD0_fil ) ;
1287double dBrDCS_1_fil = sqrt( brDCS_1_fil * ( 1. - brDCS_1_fil ) / nD0_fil ) ;
1288double dBrDCS_3_fil = sqrt( brDCS_3_fil * ( 1. - brDCS_3_fil ) / nD0_fil ) ;
1289double brCPPlus_fil = fil_CPPlus / nD0_fil ;
1290double dBrCPPlus_fil = sqrt( brCPPlus_fil * ( 1. - brCPPlus_fil ) / nD0_fil ) ;
1291double brCPMinus_fil = fil_CPMinus / nD0_fil ;
1292double dBrCPMinus_fil = sqrt( brCPMinus_fil * ( 1. - brCPMinus_fil ) / nD0_fil ) ;
1293double brDalitz_fil = fil_Dalitz / nD0_fil ;
1294double dBrDalitz_fil = sqrt( brDalitz_fil * ( 1. - brDalitz_fil ) / nD0_fil ) ;
1295double fdBrDalitz_fil = dBrDalitz_fil / brDalitz_fil ;
1296double dalitzNumer1Norm_fil = dalitzNumer1_fil / nD0_fil ;
1297double dDalitzNumer1Norm_fil = dalitzNumer1Norm_fil * fdBrDalitz_fil ;
1298double dalitzNumer2Norm_fil = dalitzNumer2_fil / nD0_fil ;
1299double dDalitzNumer2Norm_fil = dalitzNumer2Norm_fil * fdBrDalitz_fil ;
1300double dalitzDenomNorm_fil = dalitzDenom_fil / nD0_fil ;
1301double dDalitzDenomNorm_fil = dalitzDenomNorm_fil * fdBrDalitz_fil ;
1302
1303double br2Pip2Pim_fil = fil_2Pip2Pim / nD0_fil ;
1304double dBr2Pip2Pim_fil = sqrt( br2Pip2Pim_fil * ( 1. - br2Pip2Pim_fil ) / nD0_fil ) ;
1305double fdBr2Pip2Pim_fil = dBr2Pip2Pim_fil / br2Pip2Pim_fil ;
1306double Pip2Pim2Numer1Norm_fil = Pip2Pim2Numer1_fil / nD0_fil ;
1307double d2Pip2PimNumer1Norm_fil = Pip2Pim2Numer1Norm_fil * fdBr2Pip2Pim_fil ;
1308double Pip2Pim2Numer2Norm_fil = Pip2Pim2Numer2_fil / nD0_fil ;
1309double d2Pip2PimNumer2Norm_fil = Pip2Pim2Numer2Norm_fil * fdBr2Pip2Pim_fil ;
1310double Pip2Pim2DenomNorm_fil = Pip2Pim2Denom_fil / nD0_fil ;
1311double d2Pip2PimDenomNorm_fil = Pip2Pim2DenomNorm_fil * fdBr2Pip2Pim_fil ;
1312
1313log << MSG::INFO <<
1314"Original number of SL decays: " << m_nSL <<" ==> Br(SL) = " << brSL << " +- " << dBrSL << endl <<
1315"Filtered number of SL decays: " << fil_SL <<" ==> Br(SL) = " << brSL_fil << " +- " << dBrSL_fil << endl <<
1316"Original number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: " <<m_nFlavoredCFD0_0 << "/" << m_nFlavoredCSD0 << "/" << m_nFlavoredDCSD0_0 << endl <<
1317" ==> Br(CF) = " << brCF_0 << " +- " << dBrCF_0 << endl <<
1318"Original number of D0->CF1/DCS1 or D0bar->CFbar3/DCSbar3: " <<m_nFlavoredCFD0_1 << "/" << m_nFlavoredDCSD0_1 << endl <<
1319" ==> Br(CF) = " << brCF_1 << " +- " << dBrCF_1 << endl <<
1320"Original number of D0->CF1/DCS1 or D0bar->CFbar3/DCSbar3: " <<m_nFlavoredCFD0_3 << "/" << m_nFlavoredDCSD0_3 << endl <<
1321" ==> Br(CF) = " << brCF_3 << " +- " << dBrCF_3 << endl <<
1322" ==> Br(CS) = " << brCS << " +- " << dBrCS << endl <<
1323" ==> Br(DCS) = " << brDCS_0 << " +- " << dBrDCS_0 << endl <<
1324"Filtered number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: " <<fil_FlavoredCFD0_0 << "/" << fil_FlavoredCSD0 << "/" << fil_FlavoredDCSD0_0 << endl <<
1325" ==> Br(DCS) = " << brDCS_1 << " +- " << dBrDCS_1 << endl <<
1326"Filtered number of D0->CF1/DCS1 or D0bar->CFbar3/DCSbar3: " <<fil_FlavoredCFD0_1 << "/" << fil_FlavoredDCSD0_1 << endl <<
1327" ==> Br(DCS) = " << brDCS_3 << " +- " << dBrDCS_3 << endl <<
1328"Filtered number of D0->CF1/DCS1 or D0bar->CFbar3/DCSbar3: " <<fil_FlavoredCFD0_3 << "/" << fil_FlavoredDCSD0_3 << endl <<
1329" ==> Br(CF) = " << brCF_0_fil << " +- " << dBrCF_0_fil << endl <<
1330" ==> Br(CF) = " << brCF_1_fil << " +- " << dBrCF_1_fil << endl <<
1331" ==> Br(CF) = " << brCF_3_fil << " +- " << dBrCF_3_fil << endl <<
1332" ==> Br(CS) = " << brCS_fil << " +- " << dBrCS_fil << endl <<
1333" ==> Br(DCS) = " << brDCS_0_fil << " +- " << dBrDCS_0_fil << endl <<
1334" ==> Br(DCS) = " << brDCS_1_fil << " +- " << dBrDCS_1_fil << endl <<
1335" ==> Br(DCS) = " << brDCS_3_fil << " +- " << dBrDCS_3_fil << endl <<
1336
1337"Original number of CP+/-: " << m_nCPPlus << "/" << m_nCPMinus <<endl <<
1338" ==> Br(CP+) = " << brCPPlus << " +- " << dBrCPPlus <<endl <<
1339" ==> Br(CP-) = " << brCPMinus << " +- " << dBrCPMinus << endl <<
1340"Filtered number of CP+/-: " << fil_CPPlus << "/" << fil_CPMinus <<endl <<
1341" ==> Br(CP+) = " << brCPPlus_fil << " +- " << dBrCPPlus_fil <<endl <<
1342" ==> Br(CP-) = " << brCPMinus_fil << " +- " << dBrCPMinus_fil << endl <<
1343
1344"Original number of Dalitz: " << m_nDalitz << endl <<
1345" ==> Br(Dalitz) = " << brDalitz << " +- " << dBrDalitz <<endl <<
1346" y contrib. numer1 " << dalitzNumer1Norm << " +- " << dDalitzNumer1Norm << endl <<
1347" numer2 " << dalitzNumer2Norm << " +- " << dDalitzNumer2Norm << endl <<
1348" denom " << dalitzDenomNorm << " +- " << dDalitzDenomNorm << endmsg <<
1349"Filtered number of Dalitz: " << fil_Dalitz << endl <<
1350" ==> Br(Dalitz) = " << brDalitz_fil << " +- " << dBrDalitz_fil <<endl <<
1351" y contrib. numer1 " << dalitzNumer1Norm_fil << " +- " << dDalitzNumer1Norm_fil << endl <<
1352" numer2 " << dalitzNumer2Norm_fil << " +- " << dDalitzNumer2Norm_fil << endl <<
1353" denom " << dalitzDenomNorm_fil << " +- " << dDalitzDenomNorm_fil << endmsg <<
1354
1355"Original number of 2Pip2Pim: " << m_n2Pip2Pim << endl <<
1356" ==> Br(2Pip2Pim) = " << br2Pip2Pim << " +- " << dBr2Pip2Pim <<endl <<
1357" y contrib. numer1 " << Pip2Pim2Numer1Norm << " +- " << d2Pip2PimNumer1Norm << endl <<
1358" numer2 " << Pip2Pim2Numer2Norm << " +- " << d2Pip2PimNumer2Norm << endl <<
1359" denom " << Pip2Pim2DenomNorm << " +- " << d2Pip2PimDenomNorm << endmsg <<
1360"Filtered number of 2Pip2Pim: " << fil_2Pip2Pim << endl <<
1361" ==> Br(2Pip2Pim) = " << br2Pip2Pim_fil << " +- " << dBr2Pip2Pim_fil <<endl <<
1362" y contrib. numer1 " << Pip2Pim2Numer1Norm_fil << " +- " << d2Pip2PimNumer1Norm_fil << endl <<
1363" numer2 " << Pip2Pim2Numer2Norm_fil << " +- " << d2Pip2PimNumer2Norm_fil << endl <<
1364" denom " << Pip2Pim2DenomNorm_fil << " +- " << d2Pip2PimDenomNorm_fil << endmsg ;
1365
1366
1367HepMatrix differencetoWeight(17,17,0);
1368for( int i = 0 ; i < 16 ; ++i ) {
1369 for( int j = 0 ; j < 16 ; ++j ) {
1370 if (m_modeCounter[i][j] != 0) {
1371 differencetoWeight[i][j] = m_keptModeCounter[i][j] / m_modeCounter[i][j] - m_weights[i][j]; } } }
1372
1373 log << MSG::INFO << "Weight matrix" << m_weights << endmsg ;
1374 log << MSG::INFO << "D0 Modes before filter" << m_modeCounter << endmsg ;
1375 log << MSG::INFO << "D0 Modes after filter" << m_keptModeCounter << endmsg ;
1376 log << MSG::INFO << "Compare percentage difference to weights " << differencetoWeight << endmsg ;
1377
1378 //Calculating and comparing Y before/after the filter
1379 // To calculate y, we want half the CS BF
1380 brCS /= 2. ;
1381 dBrCS /= 2. ;
1382 brCS_fil /= 2 ;
1383 dBrCS_fil /= 2 ;
1384
1385 double y, dy, y_fil, dy_fil ;
1386 double y_cpsl, dy_cpsl, y_fil_cpsl, dy_fil_cpsl ;
1387
1388 //Original MC
1389 double rCF_0 = m_rCF_0 ;
1390 double zCF_0 = m_zCF_0 ;
1391 double rCF2_0 = rCF_0 * rCF_0 ;
1392 double zCF2_0 = zCF_0 * zCF_0 ;
1393 double rCF_1 = m_rCF_1 ;
1394 double zCF_1 = m_zCF_1 ;
1395 double rCF2_1 = rCF_1 * rCF_1 ;
1396 double zCF2_1 = zCF_1 * zCF_1 ;
1397 double rCF_3 = m_rCF_3 ;
1398 double zCF_3 = m_zCF_3 ;
1399 double rCF2_3 = rCF_3 * rCF_3 ;
1400 double zCF2_3 = zCF_3 * zCF_3 ;
1401 double rCS = m_rCS ;
1402 double zCS = m_zCS ;
1403 double rCS2 = rCS * rCS ;
1404 double zCS2 = zCS * zCS ;
1405
1406 double wCF_0 = ( m_wCFSign_0 ? 1. : -1. ) * sqrt( 4. - zCF2_0 ) ;
1407 double wCF_1 = ( m_wCFSign_1 ? 1. : -1. ) * sqrt( 4. - zCF2_1 ) ;
1408 double wCF_3 = ( m_wCFSign_3 ? 1. : -1. ) * sqrt( 4. - zCF2_3 ) ;
1409 double wCS = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 ) ;
1410
1411 double numFactCF_0 =
1412 -rCF_0 * zCF_0 * ( 1. - 0.5 * rCF_0 * wCF_0 * m_x ) ;
1413 double numFactCF_1 =
1414 -rCF_1 * zCF_1 * ( 1. - 0.5 * rCF_1 * wCF_1 * m_x ) ;
1415 double numFactCF_3 =
1416 -rCF_3 * zCF_3 * ( 1. - 0.5 * rCF_3 * wCF_3 * m_x ) ;
1417 double numFactCS =
1418 -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x ) ;
1419 double denFactCF_0 = 0.5 * rCF2_0 * zCF2_0 ;
1420 double denFactCF_1 = 0.5 * rCF2_1 * zCF2_1 ;
1421 double denFactCF_3 = 0.5 * rCF2_3 * zCF2_3 ;
1422 double denFactCS = 0.5 * rCS2 * zCS2 ;
1423 double dalitzNumerNorm = dalitzNumer1Norm + dalitzNumer2Norm ;
1424 double Pip2Pim2NumerNorm = Pip2Pim2Numer1Norm + Pip2Pim2Numer2Norm ;
1425
1426 double num =
1427 brCPPlus - brCPMinus + brCF_0 * numFactCF_0 + brCF_1 * numFactCF_1 + brCF_3 * numFactCF_3 + brCS * numFactCS +
1428 dalitzNumerNorm + Pip2Pim2NumerNorm;
1429 double den =
1430 1. - brCPPlus - brCPMinus - brCF_0 * denFactCF_0 - brCF_1 * denFactCF_1 - brCF_3 * denFactCF_3 - brCS * denFactCS +
1431 dalitzDenomNorm + Pip2Pim2DenomNorm;
1432 y = num / den ;
1433 dy = sqrt(
1434 ( num + den ) * ( num + den ) * dBrCPPlus * dBrCPPlus +
1435 ( num - den ) * ( num - den ) * dBrCPMinus * dBrCPMinus +
1436 ( numFactCF_0 * den + denFactCF_0 * num ) * dBrCF_0 *
1437 ( numFactCF_0 * den + denFactCF_0 * num ) * dBrCF_0 +
1438 ( numFactCF_1 * den + denFactCF_1 * num ) * dBrCF_1 *
1439 ( numFactCF_1 * den + denFactCF_1 * num ) * dBrCF_1 +
1440 ( numFactCF_3 * den + denFactCF_3 * num ) * dBrCF_3 *
1441 ( numFactCF_3 * den + denFactCF_3 * num ) * dBrCF_3 +
1442 ( numFactCS * den + denFactCS * num ) * dBrCS *
1443 ( numFactCS * den + denFactCS * num ) * dBrCS +
1444 ( dalitzNumerNorm * den + dalitzDenomNorm * num ) * fdBrDalitz *
1445 ( dalitzNumerNorm * den + dalitzDenomNorm * num ) * fdBrDalitz +
1446 ( Pip2Pim2NumerNorm * den + Pip2Pim2DenomNorm * num ) * fdBrDalitz *
1447 ( Pip2Pim2NumerNorm * den + Pip2Pim2DenomNorm * num ) * fdBrDalitz
1448 ) / ( den * den ) ;
1449
1450 double n_cpplussl = m_modeCounter[6][4] + m_modeCounter[6][5] + m_modeCounter[4][6] + m_modeCounter[5][6] ;
1451 double n_cpminussl = m_modeCounter[7][4] + m_modeCounter[7][5] + m_modeCounter[4][7] + m_modeCounter[5][7] ;
1452 double a_cpsl = ( n_cpplussl * m_nCPMinus ) / ( n_cpminussl * m_nCPPlus ) ;
1453 double b_cpsl = ( n_cpminussl * m_nCPPlus ) / ( n_cpplussl * m_nCPMinus ) ;
1454 y_cpsl = 0.25 * ( a_cpsl - b_cpsl ) ;
1455 dy_cpsl = 0.25 * ( a_cpsl + b_cpsl ) * sqrt( ( 1/n_cpplussl ) + ( 1/n_cpminussl ) + ( 1/m_nCPPlus ) + ( 1/m_nCPMinus ) ) ;
1456
1457
1458 //Filtered MC
1459 double dalitzNumerNorm_fil = dalitzNumer1Norm_fil + dalitzNumer2Norm_fil ;
1460 double Pip2Pim2NumerNorm_fil = Pip2Pim2Numer1Norm_fil + Pip2Pim2Numer2Norm_fil ;
1461 double num_fil =
1462 brCPPlus_fil - brCPMinus_fil + brCF_0_fil * numFactCF_0 + brCF_1_fil * numFactCF_1 + brCF_3_fil * numFactCF_3 + brCS_fil * numFactCS +
1463 dalitzNumerNorm_fil + Pip2Pim2NumerNorm_fil;
1464 double den_fil =
1465 1. - brCPPlus_fil - brCPMinus_fil - brCF_0_fil * denFactCF_0 - brCF_1_fil * denFactCF_1 - brCF_3_fil * denFactCF_3 - brCS_fil * denFactCS +
1466 dalitzDenomNorm_fil + Pip2Pim2DenomNorm_fil;
1467 y_fil = num_fil / den_fil ;
1468 dy_fil = sqrt(
1469 ( num_fil + den_fil ) * ( num_fil + den_fil ) * dBrCPPlus_fil * dBrCPPlus_fil +
1470 ( num_fil - den_fil ) * ( num_fil - den_fil ) * dBrCPMinus_fil * dBrCPMinus_fil +
1471 ( numFactCF_0 * den_fil + denFactCF_0 * num_fil ) * dBrCF_0_fil *
1472 ( numFactCF_0 * den_fil + denFactCF_0 * num_fil ) * dBrCF_0_fil +
1473 ( numFactCF_1 * den_fil + denFactCF_1 * num_fil ) * dBrCF_1_fil *
1474 ( numFactCF_1 * den_fil + denFactCF_1 * num_fil ) * dBrCF_1_fil +
1475 ( numFactCF_3 * den_fil + denFactCF_3 * num_fil ) * dBrCF_3_fil *
1476 ( numFactCF_3 * den_fil + denFactCF_3 * num_fil ) * dBrCF_3_fil +
1477 ( numFactCS * den_fil + denFactCS * num_fil ) * dBrCS_fil *
1478 ( numFactCS * den_fil + denFactCS * num_fil ) * dBrCS_fil +
1479 ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil * num_fil ) * fdBrDalitz_fil *
1480 ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil * num_fil ) * fdBrDalitz_fil +
1481 ( Pip2Pim2NumerNorm_fil * den_fil + Pip2Pim2DenomNorm_fil * num_fil ) * fdBr2Pip2Pim_fil *
1482 ( Pip2Pim2NumerNorm_fil * den_fil + Pip2Pim2DenomNorm_fil * num_fil ) * fdBr2Pip2Pim_fil
1483 ) / ( den_fil * den_fil ) ;
1484
1485 double fil_cpplussl = m_keptModeCounter[6][4] + m_keptModeCounter[6][5] + m_keptModeCounter[4][6] + m_keptModeCounter[5][6] ;
1486 double fil_cpminussl = m_keptModeCounter[7][4] + m_keptModeCounter[7][5] + m_keptModeCounter[4][7] + m_keptModeCounter[5][7] ;
1487 a_cpsl = ( fil_cpplussl * fil_CPMinus ) / ( fil_cpminussl * fil_CPPlus ) ;
1488 b_cpsl = ( fil_cpminussl * fil_CPPlus ) / ( fil_cpplussl * fil_CPMinus ) ;
1489 y_fil_cpsl = 0.25 * ( a_cpsl - b_cpsl ) ;
1490 dy_fil_cpsl = 0.25 * ( a_cpsl + b_cpsl ) * sqrt( ( 1/n_cpplussl ) + ( 1/n_cpminussl ) + ( 1/fil_CPPlus ) + ( 1/fil_CPMinus ) ) ;
1491
1492
1493
1494 log << MSG::INFO <<"BR Method : Parent MC ==> y = " << y << " +- " << dy << endmsg ;
1495 log << MSG::INFO <<"BR Method : Filtered MC ==> y = " << y_fil << " +- " << dy_fil << endmsg ;
1496 log << MSG::INFO <<"CP-SL doubletag : Parent MC ==> y = " << y_cpsl << " +- " << dy_cpsl <<endl<<"Previous line should be equivalent with 0 as parent MC not correalated"<< endmsg ;
1497 log << MSG::INFO <<"CP-SL doubletag : Filtered MC ==> y = " << y_fil_cpsl << " +- " << dy_fil_cpsl << endmsg ;
1498
1499 }*/
1500
1501return StatusCode::SUCCESS;
1502}
1503
1504
1505///////////////////////////////////////////////////////////////////////////////////////////
1506//////////////////////////////////////////////////////////////////////////////////////////
1507vector<double> QCMCFilter::findD0Decay(int charm)
1508{
1509 MsgStream log(msgSvc(), name());
1510 log << MSG::INFO <<" In findD0Decay " << endmsg ;
1511 vector<double> d0_info(9);
1512 for(int i=0; i< 8; i++) d0_info[i]=0;
1513 double decay_mode = 18; //d0_info[0] = decay mode;
1514 double r2D0 = -1; //d0_info[1] = r2D0;
1515 double deltaD0 = -1; //d0_info[2] = deltaD0;
1516 double RD0 = 1; //d0_info[3] = RD0;
1517 double FCP = 1; //d0_info[4] = FCP;
1518
1519 double num[26];
1520 for(int i=0; i< 26; i++) num[i]=0;
1521 // Number of each partile type
1522 int kPiPlus = 0;// num[0] = Pi+, a0+
1523 int kPiMinus = 1;// num[1] = Pi-, a0-
1524 int kPi0 = 2;// num[2] = Pi0
1525 int kEta = 3;// num[3] = Eta
1526 int kEtaPrime = 4;// num[4] = Eta Prime
1527 int kNeutralScalar = 5;// num[5] = Neutral Scalar : f0, f'0, f0(1500), a0, (f_2 can be included here because it is only used in f_2 pi0 decay)
1528 int kUFVPlus = 6;// num[6] = unflavored positive (axial) vectors: rho+, a1+, rho(2S)+
1529 int kUFVMinus = 7;// num[7] = unflavored negative (axial) vectors: rho-, a1-, rho(2S)-
1530 int kRho0 = 8;// num[8] = Rho0, Rho(2S)0
1531 int kOmega = 9;// num[9] = Omega
1532 int kPhi = 10;// num[10] = Phi
1533 int kKPlus = 11;// num[11] = K+
1534 int kKMinus = 12;// num[12] = K-
1535 int kK0S = 13;// num[13] = K short
1536 int kK0L = 14;// num[14] = K long
1537 int kFVPlus = 15;// num[15] = flavored positive (axial) vectors: K*+, K_1+, K_2*+
1538 int kFVMinus = 16;// num[16] = flavored negative (axial) vectors: K*-, K_1-, K_2*-
1539 int kKStar0 = 17;// num[17] = K*0
1540 int kKStar0Bar = 18;// num[18] = anti-K*0
1541 int kK10 = 19;// num[19] = K1_0
1542 int kK10Bar = 20;// num[20] = anti-K1_0
1543 int kLeptonPlus = 21;// num[21] = LeptonPlus: e+, mu+
1544 int kLeptonMinus = 22;// num[22] = LeptonMinus: e-, mu-
1545 int kNeutrino = 23;// num[23] = Neutrino: nu, nubar
1546 int kGamma = 24;// num[24] = gamma
1547 int kUnknown = 25;// num[25] = Not identified particle
1548 //int kGammaFSR GammaFSR is ignored as we are interested in the PreFSR quantum #
1549
1550 // Get McParticle collection
1551 SmartDataPtr<Event::McParticleCol> mcParticles(eventSvc(), EventModel::MC::McParticleCol);
1552 int d0_pdg = charm == 1 ? kD0ID : kD0BarID;
1553
1554 HepLorentzVector piPlus, piMinus, k0;
1555 vector<HepLorentzVector> piPlus_4pi, piMinus_4pi, pi0_p4;
1556 vector<HepLorentzVector> kPlus_kkpipi, kMinus_kkpipi;
1557 piPlus_4pi.clear();piMinus_4pi.clear(), pi0_p4.clear();
1558 kPlus_kkpipi.clear();kMinus_kkpipi.clear();
1559
1560 for( Event::McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end(); it++)
1561 {
1562 int pdg_code = (*it)->particleProperty();
1563 if ((*it)->primaryParticle()) continue;
1564 Event::McParticle it2 = (*it)->mother();
1565 int mother_pdg = 0;
1566
1567 //finds if the particle is related to the d0 we need
1568 mother_pdg = it2.particleProperty();
1569
1570 //special cases
1571 if (pdg_code == kKStar0ID || pdg_code == kKStar0BarID || pdg_code == kK0ID || pdg_code == kK0BarID || pdg_code == kKDPStar0ID || pdg_code == kKDPStar0BarID
1572 || pdg_code == kK10ID || pdg_code == kK10BarID
1573 || pdg_code == kK0Star0ID || pdg_code == kK0Star0BarID || pdg_code == kK0StarPlusID || pdg_code == kK0BarID ) continue; //don't count special intermediate states
1574
1575 if (mother_pdg == kK0ID || mother_pdg == kK0BarID) {//Looks another step up if mother was k0 or k0bar
1576 it2 = it2.mother();
1577 mother_pdg = it2.particleProperty(); }
1578
1579 if (mother_pdg == kKStar0ID || mother_pdg == kKStar0BarID) {//Looks another step up if mother was k*0 or k*0bar
1580 //if code is found to be K*0 -> K+Pi- we save it as a KStar0 and K*0Bar -> K-Pi+ as KStar0Bar
1581 //if is it a neutral decay K*0(Bar)->K0(bar) Pi0 or Gamma. We save them as K0(bar) and Pi0 or Gamma.
1582 if (pdg_code == kPiPlusID || pdg_code == kPiMinusID) continue;
1583 if (mother_pdg == kKStar0ID && pdg_code == kKPlusID) pdg_code = kKStar0ID;
1584 if (mother_pdg == kKStar0BarID && pdg_code == kKMinusID) pdg_code = kKStar0BarID;
1585 it2 = it2.mother();
1586 mother_pdg = it2.particleProperty();
1587 }
1588
1589 if ( mother_pdg == kK0Star0ID || mother_pdg == kK0Star0BarID ) { // Looks another step up if mother was K_0
1590 it2 = it2.mother();
1591 mother_pdg = it2.particleProperty();
1592 }
1593
1594 if (mother_pdg == kK10ID || mother_pdg == kK10BarID) {//Looks another step up if mother was k_10 or k_10bar
1595 it2 = it2.mother();
1596 mother_pdg = it2.particleProperty();
1597 }
1598
1599
1600
1601 //Identifies what particles are the daughters
1602 if (mother_pdg == d0_pdg) //Found daughter
1603 {
1604 if (pdg_code == kPiPlusID || pdg_code == kA0PlusID ) num[0]++;
1605 else if (pdg_code == kPiMinusID || pdg_code == kA0MinusID ) num[1]++;
1606 else if (pdg_code == kPi0ID ) num[2]++;
1607 else if (pdg_code == kEtaID ) num[3]++;
1608 else if (pdg_code == kEtaPrimeID ) num[4]++;
1609 else if (pdg_code == kF0ID || pdg_code == kA00ID
1610 || pdg_code == kFPrime0ID|| pdg_code == kF0m1500ID
1611 || pdg_code == kF2ID) num[5]++;
1612 else if (pdg_code == kRhoPlusID || pdg_code == kRho2SPlusID
1613 || pdg_code == kA1PlusID ) num[6]++;
1614 else if (pdg_code == kRhoMinusID || pdg_code == kRho2SMinusID
1615 || pdg_code == kA1MinusID) num[7]++;
1616 else if (pdg_code == kRho0ID || pdg_code == kRho2S0ID) num[8]++;
1617 else if (pdg_code == kOmegaID ) num[9]++;
1618 else if (pdg_code == kPhiID ) num[10]++;
1619 else if (pdg_code == kKPlusID ) num[11]++;
1620 else if (pdg_code == kKMinusID ) num[12]++;
1621 else if (pdg_code == kK0SID ) num[13]++;
1622 else if (pdg_code == kK0LID ) num[14]++;
1623 else if (pdg_code == kKStarPlusID || pdg_code == kK1PlusID
1624 || pdg_code == kK2StarPlusID || pdg_code == kK1PrimePlusID
1625 || pdg_code == kK0StarPlusID) num[15]++;
1626 else if (pdg_code == kKStarMinusID || pdg_code == kK1MinusID
1627 || pdg_code == kK2StarMinusID || pdg_code == kK1PrimeMinusID
1628 || pdg_code == kK0StarMinusID ) num[16]++;
1629 else if (pdg_code == kKStar0ID ) num[17]++;
1630 else if (pdg_code == kKStar0BarID ) num[18]++;
1631 else if (pdg_code == kK10ID || pdg_code == kK1Prime0ID) num[19]++;
1632 else if (pdg_code == kK10BarID || pdg_code == kK1Prime0BarID) num[20]++;
1633 else if (pdg_code == kEPlusID || pdg_code == kMuPlusID ) num[21]++;
1634 else if (pdg_code == kEMinusID || pdg_code == kMuMinusID ) num[22]++;
1635 else if (pdg_code == kNuEID || pdg_code == kNuEBarID
1636 || pdg_code == kNuMuID || pdg_code == kNuMuBarID ) num[23]++;
1637 else if (pdg_code == kGammaID ) num[24]++;
1638 else if (pdg_code == kGammaFSRID ) continue;
1639 else {
1640 num[25]++;
1641 cout <<"Unknown particle: "<< pdg_code << endl;
1642 }
1643
1644 //if (pdg_code == kA10ID ) num[0]++; // not present
1645
1646 //Enter Momentums for Dalitz canidates
1647 //It would be prefered if we could get PreFSR Momentum
1648 //m_D0daughter_P4.push_back( (*it)->initialFourMomentum().vect() );
1649 //m_D0daughter_PID.push_back( pdg_code );
1650
1651 if (pdg_code == kPiPlusID) piPlus.setVectM((*it)->initialFourMomentum().vect(), xmpion);
1652 if (pdg_code == kPiMinusID) piMinus.setVectM((*it)->initialFourMomentum().vect(), xmpion);
1653 if (pdg_code == kK0SID) k0.setVectM((*it)->initialFourMomentum().vect(), xmk0);
1654 if (pdg_code == kK0LID) k0.setVectM((*it)->initialFourMomentum().vect(), xmk0);
1655 HepLorentzVector daughterP4;
1656 if(piPlus_4pi.size()>2||piMinus_4pi.size()>2)continue;
1657 if(kPlus_kkpipi.size()>1||kMinus_kkpipi.size()>1)continue;
1658 //log << MSG::INFO << "in D0->pipipipi;" << endmsg;
1659 if (pdg_code == kPiPlusID) {daughterP4.setVectM((*it)->initialFourMomentum().vect(), xmpion);piPlus_4pi.push_back(daughterP4);}
1660 if (pdg_code == kPiMinusID) {daughterP4.setVectM((*it)->initialFourMomentum().vect(), xmpion);piMinus_4pi.push_back(daughterP4);}
1661 if (pdg_code == kKPlusID) {daughterP4.setVectM((*it)->initialFourMomentum().vect(), xmkaon);kPlus_kkpipi.push_back(daughterP4);}
1662 if (pdg_code == kKMinusID) {daughterP4.setVectM((*it)->initialFourMomentum().vect(), xmkaon);kMinus_kkpipi.push_back(daughterP4);}
1663 if (pdg_code == kPi0ID) {daughterP4.setVectM((*it)->initialFourMomentum().vect(), xmpi0);pi0_p4.push_back(daughterP4);}
1664 }
1665 }
1666
1667 // D decay daughters have been tabulated. Now, classify decay.
1668 int nDaughters = 0 ;
1669 for( int i = 0 ; i < 26 ; i++ )
1670 {
1671 nDaughters += num[ i ] ;
1672 }
1673
1674 int nUFP0 = num[ kPi0 ] + num[ kEta ] + num[ kEtaPrime ] ;
1675 int nUFV0 = num[ kRho0 ] + num[ kPhi ] + num[ kOmega ] + num[ kGamma ] ;
1676 int nFV0 = num[ kKStar0 ] + num[ kK10 ] ;
1677 int nFV0Bar = num[ kKStar0Bar ] + num[ kK10Bar ] ;
1678 int nStrange = num[ kKMinus ] + num[ kFVMinus ] + nFV0Bar ;
1679 int nAntiStrange = num[ kKPlus ] + num[ kFVPlus ] + nFV0 ;
1680 int nCPPlusEig = num[ kNeutralScalar ] + num[ kRho0 ] + num[ kOmega ] + num[ kPhi ] + num[ kK0S ] + num[ kGamma ] ;
1681 int nCPMinusEig = num[ kPi0 ] + num[ kEta ] + num[ kEtaPrime ] + num[ kK0L ] ;
1682 int nCPEig = nCPPlusEig + nCPMinusEig ;
1683 int nChargedPiK = num[ kPiPlus ] + num[ kPiMinus ] + num[ kKPlus ] + num[ kKMinus ] ;
1684 int nK0 = num[ kK0S ] + num[ kK0L ] ;
1685
1686
1687 //check Dalitz modes: KsPi+Pi- or KlPi+Pi-
1688 double mnm_gen = 0. ;
1689 double mpn_gen = 0. ;
1690 double mpm_gen = 0. ;
1691 bool isKsPiPi = false ;
1692 bool isKsPiPiPi0 = false ;
1693
1694 //check K3Pi modes:
1695
1696
1697 if( nK0 == 1 && num[ kPiPlus ] == 1 && num[ kPiMinus ] && nDaughters == 3 )
1698 {
1699 decay_mode = kDalitz ;
1700 //k0.boost(-0.011, 0, 0);
1701 //piMinus.boost(-0.011, 0, 0);
1702 //piPlus.boost(-0.011, 0, 0);
1703 //mnm_gen = (k0 + piMinus).m2();
1704 //mpn_gen = (k0 + piPlus).m2();
1705 //mpm_gen = (piPlus + piMinus).m2();
1706 if( num[ kK0S ] == 1) isKsPiPi = true ;
1707 }
1708 if( num[ kPiPlus ] == 2&&num[ kPiMinus ]==2 && nDaughters == 4 )
1709 {
1710 decay_mode = k2Pip2Pim;
1711 }
1712 if( num[ kPiPlus ] == 1&&num[ kPiMinus ]==1 && num[ kPi0 ]==2 && nDaughters == 4 )
1713 {
1714 decay_mode = kPipPim2Pi0;
1715 }
1716 if( nK0 == 1 && num[ kPiPlus ] == 1&&num[ kPiMinus ]==1 && num[ kPi0 ]==1 && nDaughters == 4 )
1717 {
1718 decay_mode = kK0PiPiPi0;
1719 if( num[ kK0S ] == 1) isKsPiPiPi0 = true ;
1720 }
1721 if( num[ kPiPlus ] == 1&&num[ kPiMinus ]==1 &&num[ kKPlus ] == 1&&num[ kKMinus ]==1 && nDaughters == 4 )
1722 {
1723 decay_mode = kKKPiPi;
1724 }
1725
1726 // The order of if-statements below matters.
1727 if( decay_mode == kDalitz )
1728 {
1729 ++m_nDalitz ;
1730
1731 // Dalitz t(8) ;
1732 // TComplex D0, D0bar;
1733 // if (m_dalitzModel == 0) { //Cleo model
1734 // D0 = t.CLEO_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1735 // D0bar = t.CLEO_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1736 // if (m_dalitzModel == 1) { //Babar model
1737 // D0 = t.Babar_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1738 // D0bar = t.Babar_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1739 // if (m_dalitzModel == 2) { //Belle model
1740 // D0 = t.Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1741 // D0bar = t.Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1742
1743 complex<double> D0, D0bar;
1744 vector<double> k0l;vector<double> pip;vector<double> pim;
1745 k0l.clear();pip.clear();pim.clear();
1746 k0l.push_back(k0[0]);k0l.push_back(k0[1]);k0l.push_back(k0[2]);k0l.push_back(k0[3]);
1747 pip.push_back(piPlus[0]);pip.push_back(piPlus[1]);pip.push_back(piPlus[2]);pip.push_back(piPlus[3]);
1748 pim.push_back(piMinus[0]);pim.push_back(piMinus[1]);pim.push_back(piMinus[2]);pim.push_back(piMinus[3]);
1749
1750 if(isKsPiPi){
1751 D0ToKSpipi2018 tKSpipi;tKSpipi.init();
1752 D0 = tKSpipi.Amp_PFT(k0l,pip,pim);
1753
1754 vector<double> cpk0l;vector<double> cppip;vector<double> cppim;
1755 cpk0l.clear();pip.clear();pim.clear();
1756 cpk0l.push_back(-k0l[0]);cpk0l.push_back(-k0l[1]);cpk0l.push_back(-k0l[2]);cpk0l.push_back(k0l[3]);
1757 cppim.push_back(-piPlus[0]);cppim.push_back(-piPlus[1]);cppim.push_back(-piPlus[2]);cppim.push_back(piPlus[3]);
1758 cppip.push_back(-piMinus[0]);cppip.push_back(-piMinus[1]);cppip.push_back(-piMinus[2]);cppip.push_back(piMinus[3]);
1759 D0bar = tKSpipi.Amp_PFT(cpk0l, cppip, cppim);
1760 }else{
1761 D0ToKLpipi2018 tKLpipi;tKLpipi.init();//be aware of the convention in D0ToKLpipi.cxx is not the same as used in other analyses
1762 D0 = tKLpipi.Amp_PFT(k0l,pip,pim);
1763 vector<double> cpk0l;vector<double> cppip;vector<double> cppim;
1764 cpk0l.clear();pip.clear();pim.clear();
1765 cpk0l.push_back(-k0l[0]);cpk0l.push_back(-k0l[1]);cpk0l.push_back(-k0l[2]);cpk0l.push_back(k0l[3]);
1766 cppim.push_back(-piPlus[0]);cppim.push_back(-piPlus[1]);cppim.push_back(-piPlus[2]);cppim.push_back(piPlus[3]);
1767 cppip.push_back(-piMinus[0]);cppip.push_back(-piMinus[1]);cppip.push_back(-piMinus[2]);cppip.push_back(piMinus[3]);
1768 D0bar = tKLpipi.Amp_PFT(cpk0l, cppip, cppim);
1769 }
1770
1771 if( charm == 1){
1772
1773 r2D0 = std::norm(D0bar)/std::norm(D0);
1774 double temp = std::arg(D0)-std::arg(D0bar);
1775 log << MSG::INFO << "D0 calculation " << sqrt(r2D0) << " "<< temp << endmsg;
1776 while (temp < -TMath::Pi())
1777 {
1778 temp += 2.0 * TMath::Pi();
1779 }
1780 while (temp > TMath::Pi())
1781 {
1782 temp -= 2.0 * TMath::Pi();
1783 }
1784 deltaD0 = temp;
1785 deltaD0 = isKsPiPi ? deltaD0 : (deltaD0+PI);
1786
1787 double cosd = cos(deltaD0);
1788 double sind = sin(deltaD0);
1789 if( mpn_gen - mnm_gen > 0. )
1790 {
1791 m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd ;
1792 m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x ;
1793 m_dalitzDenom += -2. * r2D0 * cosd * cosd ;
1794 }
1795 }
1796 else {
1797 r2D0 = std::norm(D0)/std::norm(D0bar);
1798 double temp = std::arg(D0bar)-std::arg(D0);
1799 while (temp < -TMath::Pi())
1800 {
1801 temp += 2.0 * TMath::Pi();
1802 }
1803 while (temp > TMath::Pi())
1804 {
1805 temp -= 2.0 * TMath::Pi();
1806 }
1807 deltaD0 = temp;
1808 deltaD0 = isKsPiPi ? deltaD0 : (deltaD0+PI);
1809
1810 double cosd = cos( deltaD0 ) ;
1811 double sind = sin( deltaD0 ) ;
1812 if( mpn_gen - mnm_gen < 0. )
1813 {
1814 m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd ;
1815 m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x ;
1816 m_dalitzDenom += -2. * r2D0 * cosd * cosd ;
1817 }
1818 }
1819 }else if( decay_mode == k2Pip2Pim ){
1820
1821 ++m_n2Pip2Pim ;
1822
1823 HepLorentzVector pip1_4pi = (piPlus_4pi.at(0));
1824 HepLorentzVector pim1_4pi = (piMinus_4pi.at(0));
1825 HepLorentzVector pip2_4pi = (piPlus_4pi.at(1));
1826 HepLorentzVector pim2_4pi = (piMinus_4pi.at(1));
1827 pip1_4pi.boost(-0.011, 0, 0);
1828 pip2_4pi.boost(-0.011, 0, 0);
1829 pim1_4pi.boost(-0.011, 0, 0);
1830 pim2_4pi.boost(-0.011, 0, 0);
1831 std::complex<double> D0, D0bar;
1832 //EvtVector4R pip1( pip1_4pi[3],pip1_4pi[0],pip1_4pi[1],pip1_4pi[2]);
1833 //EvtVector4R pip2( pip2_4pi[3],pip2_4pi[0],pip2_4pi[1],pip2_4pi[2]);
1834 //EvtVector4R pim1( pim1_4pi[3],pim1_4pi[0],pim1_4pi[1],pim1_4pi[2]);
1835 //EvtVector4R pim2( pim2_4pi[3],pim2_4pi[0],pim2_4pi[1],pim2_4pi[2]);
1836 std::vector<double> pip1; pip1.clear(); pip1.push_back(pip1_4pi[0]);pip1.push_back(pip1_4pi[1]);pip1.push_back(pip1_4pi[2]);pip1.push_back(pip1_4pi[3]);
1837 std::vector<double> pim1; pim1.clear(); pim1.push_back(pim1_4pi[0]);pim1.push_back(pim1_4pi[1]);pim1.push_back(pim1_4pi[2]);pim1.push_back(pim1_4pi[3]);
1838 std::vector<double> pip2; pip2.clear(); pip2.push_back(pip2_4pi[0]);pip2.push_back(pip2_4pi[1]);pip2.push_back(pip2_4pi[2]);pip2.push_back(pip2_4pi[3]);
1839 std::vector<double> pim2; pim2.clear(); pim2.push_back(pim2_4pi[0]);pim2.push_back(pim2_4pi[1]);pim2.push_back(pim2_4pi[2]);pim2.push_back(pim2_4pi[3]);
1840
1841 D0To2pip2pim t2pip2pim;t2pip2pim.init();
1842 D0 = t2pip2pim.Amp(pip1,pim1,pip2,pim2);
1843 //EvtVector4R cppim1(pip1_4pi[3], -1*pip1_4pi[0], -1*pip1_4pi[1], -1*pip1_4pi[2]);
1844 //EvtVector4R cppip1(pim1_4pi[3], -1*pim1_4pi[0], -1*pim1_4pi[1], -1*pim1_4pi[2]);
1845 //EvtVector4R cppim2(pip2_4pi[3], -1*pip2_4pi[0], -1*pip2_4pi[1], -1*pip2_4pi[2]);
1846 //EvtVector4R cppip2(pim2_4pi[3], -1*pim2_4pi[0], -1*pim2_4pi[1], -1*pim2_4pi[2]);
1847 std::vector<double> cppip1; cppip1.clear(); cppip1.push_back(-pim1_4pi[0]);cppip1.push_back(-pim1_4pi[1]);cppip1.push_back(-pim1_4pi[2]);cppip1.push_back(pim1_4pi[3]);
1848 std::vector<double> cppim1; cppim1.clear(); cppim1.push_back(-pip1_4pi[0]);cppim1.push_back(-pip1_4pi[1]);cppim1.push_back(-pip1_4pi[2]);cppim1.push_back(pip1_4pi[3]);
1849 std::vector<double> cppip2; cppip2.clear(); cppip2.push_back(-pim2_4pi[0]);cppip2.push_back(-pim2_4pi[1]);cppip2.push_back(-pim2_4pi[2]);cppip2.push_back(pim2_4pi[3]);
1850 std::vector<double> cppim2; cppim2.clear(); cppim2.push_back(-pip2_4pi[0]);cppim2.push_back(-pip2_4pi[1]);cppim2.push_back(-pip2_4pi[2]);cppim2.push_back(pip2_4pi[3]);
1851
1852 D0bar = t2pip2pim.Amp(cppip1, cppim1, cppip2, cppim2);
1853
1854 if( charm == 1){
1855
1856 r2D0 = std::norm(D0bar)/std::norm(D0);
1857 double temp = std::arg(D0)-std::arg(D0bar);
1858 while (temp < -TMath::Pi())
1859 {
1860 temp += 2.0 * TMath::Pi();
1861 }
1862 while (temp > TMath::Pi())
1863 {
1864 temp -= 2.0 * TMath::Pi();
1865 }
1866 deltaD0 = temp;
1867
1868 double cosd = cos(deltaD0);
1869 double sind = sin(deltaD0);
1870 //if( mpn_gen - mnm_gen > 0. )
1871 //{
1872 m_2Pip2PimNumer1 += -2. * sqrt( r2D0 ) * cosd ;
1873 m_2Pip2PimNumer2 += 2. * r2D0 * cosd * sind * m_x ;
1874 m_2Pip2PimDenom += -2. * r2D0 * cosd * cosd ;
1875 //}
1876 }
1877 else {
1878 r2D0 = std::norm(D0)/std::norm(D0bar);
1879 double temp = std::arg(D0bar)-std::arg(D0);
1880 while (temp < -TMath::Pi())
1881 {
1882 temp += 2.0 * TMath::Pi();
1883 }
1884 while (temp > TMath::Pi())
1885 {
1886 temp -= 2.0 * TMath::Pi();
1887 }
1888 deltaD0 = temp;
1889
1890 double cosd = cos( deltaD0 ) ;
1891 double sind = sin( deltaD0 ) ;
1892 //if( mpn_gen - mnm_gen < 0. )
1893 //{
1894 m_2Pip2PimNumer1 += -2. * sqrt( r2D0 ) * cosd ;
1895 m_2Pip2PimNumer2 += 2. * r2D0 * cosd * sind * m_x ;
1896 m_2Pip2PimDenom += -2. * r2D0 * cosd * cosd ;
1897 //}
1898 }
1899 }else if( decay_mode == kPipPim2Pi0 ){
1900
1901 ++m_nPipPim2Pi0 ;
1902
1903
1904 HepLorentzVector pip1_4pi = (piPlus_4pi.at(0));
1905 HepLorentzVector pim1_4pi = (piMinus_4pi.at(0));
1906 HepLorentzVector pi01_p4 = (pi0_p4.at(0));
1907 HepLorentzVector pi02_p4 = (pi0_p4.at(1));
1908 pip1_4pi.boost(-0.011, 0, 0);
1909 pim1_4pi.boost(-0.011, 0, 0);
1910 pi01_p4.boost(-0.011, 0, 0);
1911 pi02_p4.boost(-0.011, 0, 0);
1912 std::complex<double> D0, D0bar;
1913 //EvtVector4R pip1( pip1_4pi[3],pip1_4pi[0],pip1_4pi[1],pip1_4pi[2]);
1914 //EvtVector4R pip2( pip2_4pi[3],pip2_4pi[0],pip2_4pi[1],pip2_4pi[2]);
1915 //EvtVector4R pim1( pim1_4pi[3],pim1_4pi[0],pim1_4pi[1],pim1_4pi[2]);
1916 //EvtVector4R pim2( pim2_4pi[3],pim2_4pi[0],pim2_4pi[1],pim2_4pi[2]);
1917 std::vector<double> pip1; pip1.clear(); pip1.push_back(pip1_4pi[0]);pip1.push_back(pip1_4pi[1]);pip1.push_back(pip1_4pi[2]);pip1.push_back(pip1_4pi[3]);
1918 std::vector<double> pim1; pim1.clear(); pim1.push_back(pim1_4pi[0]);pim1.push_back(pim1_4pi[1]);pim1.push_back(pim1_4pi[2]);pim1.push_back(pim1_4pi[3]);
1919 std::vector<double> pi01; pi01.clear(); pi01.push_back(pi01_p4[0]);pi01.push_back(pi01_p4[1]);pi01.push_back(pi01_p4[2]);pi01.push_back(pi01_p4[3]);
1920 std::vector<double> pi02; pi02.clear(); pi02.push_back(pi02_p4[0]);pi02.push_back(pi02_p4[1]);pi02.push_back(pi02_p4[2]);pi02.push_back(pi02_p4[3]);
1921
1922 D0Topippim2pi0 tpippim2pi0;tpippim2pi0.init();
1923 D0 = tpippim2pi0.Amp(pip1,pim1,pi01,pi02);
1924 //EvtVector4R cppim1(pip1_4pi[3], -1*pip1_4pi[0], -1*pip1_4pi[1], -1*pip1_4pi[2]);
1925 //EvtVector4R cppip1(pim1_4pi[3], -1*pim1_4pi[0], -1*pim1_4pi[1], -1*pim1_4pi[2]);
1926 //EvtVector4R cppim2(pip2_4pi[3], -1*pip2_4pi[0], -1*pip2_4pi[1], -1*pip2_4pi[2]);
1927 //EvtVector4R cppip2(pim2_4pi[3], -1*pim2_4pi[0], -1*pim2_4pi[1], -1*pim2_4pi[2]);
1928 std::vector<double> cppip1; cppip1.clear(); cppip1.push_back(-pim1_4pi[0]);cppip1.push_back(-pim1_4pi[1]);cppip1.push_back(-pim1_4pi[2]);cppip1.push_back(pim1_4pi[3]);
1929 std::vector<double> cppim1; cppim1.clear(); cppim1.push_back(-pip1_4pi[0]);cppim1.push_back(-pip1_4pi[1]);cppim1.push_back(-pip1_4pi[2]);cppim1.push_back(pip1_4pi[3]);
1930 std::vector<double> cppi01; cppi01.clear(); cppi01.push_back(-pi01_p4[0]);cppi01.push_back(-pi01_p4[1]);cppi01.push_back(-pi01_p4[2]);cppi01.push_back(pi01_p4[3]);
1931 std::vector<double> cppi02; cppi02.clear(); cppi02.push_back(-pi02_p4[0]);cppi02.push_back(-pi02_p4[1]);cppi02.push_back(-pi02_p4[2]);cppi02.push_back(pi02_p4[3]);
1932
1933 D0bar = tpippim2pi0.Amp(cppip1, cppim1, cppi01, cppi02);
1934
1935 if( charm == 1){
1936
1937 r2D0 = std::norm(D0bar)/std::norm(D0);
1938 double temp = std::arg(D0)-std::arg(D0bar);
1939 while (temp < -TMath::Pi())
1940 {
1941 temp += 2.0 * TMath::Pi();
1942 }
1943 while (temp > TMath::Pi())
1944 {
1945 temp -= 2.0 * TMath::Pi();
1946 }
1947 deltaD0 = temp;
1948
1949 double cosd = cos(deltaD0);
1950 double sind = sin(deltaD0);
1951 //if( mpn_gen - mnm_gen > 0. )
1952 //{
1953 m_PipPim2Pi0Numer1 += -2. * sqrt( r2D0 ) * cosd ;
1954 m_PipPim2Pi0Numer2 += 2. * r2D0 * cosd * sind * m_x ;
1955 m_PipPim2Pi0Denom += -2. * r2D0 * cosd * cosd ;
1956 //}
1957 }
1958 else {
1959 r2D0 = std::norm(D0)/std::norm(D0bar);
1960 double temp = std::arg(D0bar)-std::arg(D0);
1961 while (temp < -TMath::Pi())
1962 {
1963 temp += 2.0 * TMath::Pi();
1964 }
1965 while (temp > TMath::Pi())
1966 {
1967 temp -= 2.0 * TMath::Pi();
1968 }
1969 deltaD0 = temp;
1970
1971 double cosd = cos( deltaD0 ) ;
1972 double sind = sin( deltaD0 ) ;
1973
1974 m_PipPim2Pi0Numer1 += -2. * sqrt( r2D0 ) * cosd ;
1975 m_PipPim2Pi0Numer2 += 2. * r2D0 * cosd * sind * m_x ;
1976 m_PipPim2Pi0Denom += -2. * r2D0 * cosd * cosd ;
1977 }
1978 }else if( decay_mode == kK0PiPiPi0 )
1979 {
1980 ++m_nK0PiPiPi0 ;
1981
1982 // Dalitz t(8) ;
1983 // TComplex D0, D0bar;
1984 // if (m_dalitzModel == 0) { //Cleo model
1985 // D0 = t.CLEO_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1986 // D0bar = t.CLEO_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1987 // if (m_dalitzModel == 1) { //Babar model
1988 // D0 = t.Babar_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1989 // D0bar = t.Babar_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1990 // if (m_dalitzModel == 2) { //Belle model
1991 // D0 = t.Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1992 // D0bar = t.Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1993 HepLorentzVector pip_kspipipi0 = (piPlus_4pi.at(0));
1994 HepLorentzVector pim_kspipipi0 = (piMinus_4pi.at(0));
1995 HepLorentzVector k0_kspipipi0 = k0;
1996 HepLorentzVector pi0_kspipipi0 = (pi0_p4.at(0));
1997 pip_kspipipi0.boost(-0.011, 0, 0);
1998 pim_kspipipi0.boost(-0.011, 0, 0);
1999 k0_kspipipi0.boost(-0.011, 0, 0);
2000 pi0_kspipipi0.boost(-0.011, 0, 0);
2001
2002 complex<double> D0, D0bar;
2003 vector<double> k0l;vector<double> pip;vector<double> pim; vector<double> pi0;
2004 k0l.clear();pip.clear();pim.clear();pi0.clear();
2005 k0l.push_back(k0_kspipipi0[0]);k0l.push_back(k0_kspipipi0[1]);k0l.push_back(k0_kspipipi0[2]);k0l.push_back(k0_kspipipi0[3]);
2006 pip.push_back(pip_kspipipi0[0]);pip.push_back(pip_kspipipi0[1]);pip.push_back(pip_kspipipi0[2]);pip.push_back(pip_kspipipi0[3]);
2007 pim.push_back(pim_kspipipi0[0]);pim.push_back(pim_kspipipi0[1]);pim.push_back(pim_kspipipi0[2]);pim.push_back(pim_kspipipi0[3]);
2008 pi0.push_back(pi0_kspipipi0[0]);pi0.push_back(pi0_kspipipi0[1]);pi0.push_back(pi0_kspipipi0[2]);pi0.push_back(pi0_kspipipi0[3]);
2009
2010 if(isKsPiPiPi0){
2011 D0ToKSpipipi0 tKSpipipi0;tKSpipipi0.init();
2012 D0 = tKSpipipi0.Amp(k0l,pip,pim,pi0);
2013
2014 vector<double> cpk0l;vector<double> cppip;vector<double> cppim;vector<double> cppi0;
2015 cpk0l.clear();cppip.clear();cppim.clear();cppi0.clear();
2016 cpk0l.push_back(-k0l[0]);cpk0l.push_back(-k0l[1]);cpk0l.push_back(-k0l[2]);cpk0l.push_back(k0l[3]);
2017 cppim.push_back(-pip[0]);cppim.push_back(-pip[1]);cppim.push_back(-pip[2]);cppim.push_back(pip[3]);
2018 cppip.push_back(-pim[0]);cppip.push_back(-pim[1]);cppip.push_back(-pim[2]);cppip.push_back(pim[3]);
2019 cppi0.push_back(-pi0[0]);cppi0.push_back(-pi0[1]);cppi0.push_back(-pi0[2]);cppi0.push_back(pi0[3]);
2020 D0bar = tKSpipipi0.Amp(cpk0l, cppip, cppim, cppi0);
2021 }else{
2022 D0ToKSpipipi0 tKSpipipi0;tKSpipipi0.init();
2023 D0 = tKSpipipi0.Amp(k0l,pip,pim,pi0);
2024
2025 vector<double> cpk0l;vector<double> cppip;vector<double> cppim;vector<double> cppi0;
2026 cpk0l.clear();cppip.clear();cppim.clear();cppi0.clear();
2027 cpk0l.push_back(-k0l[0]);cpk0l.push_back(-k0l[1]);cpk0l.push_back(-k0l[2]);cpk0l.push_back(k0l[3]);
2028 cppim.push_back(-pip[0]);cppim.push_back(-pip[1]);cppim.push_back(-pip[2]);cppim.push_back(pip[3]);
2029 cppip.push_back(-pim[0]);cppip.push_back(-pim[1]);cppip.push_back(-pim[2]);cppip.push_back(pim[3]);
2030 cppi0.push_back(-pi0[0]);cppi0.push_back(-pi0[1]);cppi0.push_back(-pi0[2]);cppi0.push_back(pi0[3]);
2031 D0bar = tKSpipipi0.Amp(cpk0l, cppip, cppim, cppi0);
2032 //D0ToKLpipi tKLpipi;tKLpipi.init();//be aware of the convention in D0ToKLpipi.cxx is not the same as used in other analyses
2033 //D0 = tKLpipi.Amp_PFT(k0l,pip,pim);
2034 //vector<double> cpk0l;vector<double> cppip;vector<double> cppim;
2035 //cpk0l.clear();pip.clear();pim.clear();
2036 //cpk0l.push_back(-k0l[0]);cpk0l.push_back(-k0l[1]);cpk0l.push_back(-k0l[2]);cpk0l.push_back(k0l[3]);
2037 //cppim.push_back(-piPlus[0]);cppim.push_back(-piPlus[1]);cppim.push_back(-piPlus[2]);cppim.push_back(piPlus[3]);
2038 //cppip.push_back(-piMinus[0]);cppip.push_back(-piMinus[1]);cppip.push_back(-piMinus[2]);cppip.push_back(piMinus[3]);
2039 //D0bar = tKLpipi.Amp_PFT(cpk0l, cppip, cppim);
2040 }
2041
2042 if( charm == 1){
2043
2044 r2D0 = std::norm(D0bar)/std::norm(D0);
2045 double temp = std::arg(D0)-std::arg(D0bar);
2046 log << MSG::INFO << "D0 calculation " << sqrt(r2D0) << " "<< temp << endmsg;
2047 while (temp < -TMath::Pi())
2048 {
2049 temp += 2.0 * TMath::Pi();
2050 }
2051 while (temp > TMath::Pi())
2052 {
2053 temp -= 2.0 * TMath::Pi();
2054 }
2055 deltaD0 = temp;
2056 deltaD0 = isKsPiPiPi0 ? (deltaD0+PI) : (deltaD0);
2057
2058 double cosd = cos(deltaD0);
2059 double sind = sin(deltaD0);
2060 if( mpn_gen - mnm_gen > 0. )
2061 {
2062 m_KSPiPiPi0Numer1 += -2. * sqrt( r2D0 ) * cosd ;
2063 m_KSPiPiPi0Numer2 += 2. * r2D0 * cosd * sind * m_x ;
2064 m_KSPiPiPi0Denom += -2. * r2D0 * cosd * cosd ;
2065 }
2066 }
2067 else {
2068 r2D0 = std::norm(D0)/std::norm(D0bar);
2069 double temp = std::arg(D0bar)-std::arg(D0);
2070 while (temp < -TMath::Pi())
2071 {
2072 temp += 2.0 * TMath::Pi();
2073 }
2074 while (temp > TMath::Pi())
2075 {
2076 temp -= 2.0 * TMath::Pi();
2077 }
2078 deltaD0 = temp;
2079 deltaD0 = isKsPiPiPi0 ? (deltaD0+PI) : (deltaD0);
2080
2081 double cosd = cos( deltaD0 ) ;
2082 double sind = sin( deltaD0 ) ;
2083 if( mpn_gen - mnm_gen < 0. )
2084 {
2085 m_KSPiPiPi0Numer1 += -2. * sqrt( r2D0 ) * cosd ;
2086 m_KSPiPiPi0Numer2 += 2. * r2D0 * cosd * sind * m_x ;
2087 m_KSPiPiPi0Denom += -2. * r2D0 * cosd * cosd ;
2088 }
2089 }
2090 }else if( decay_mode == kKKPiPi ){
2091
2092 ++m_nKKPiPi ;
2093
2094 complex<double> D0, D0bar;
2095 vector<double> kp; vector<double> km; vector<double> pip;vector<double> pim;
2096 double pdau[16];
2097 pip.clear();pim.clear();
2098 kp.clear();km.clear();
2099
2100 HepLorentzVector pip_kkpipi = (piPlus_4pi.at(0));
2101 HepLorentzVector pim_kkpipi = (piMinus_4pi.at(0));
2102 HepLorentzVector kp_kkpipi = (kPlus_kkpipi.at(0));
2103 HepLorentzVector km_kkpipi = (kMinus_kkpipi.at(0));
2104
2105 pip.push_back(pip_kkpipi[0]);pip.push_back(pip_kkpipi[1]);pip.push_back(pip_kkpipi[2]);pip.push_back(pip_kkpipi[3]);
2106 pim.push_back(pim_kkpipi[0]);pim.push_back(pim_kkpipi[1]);pim.push_back(pim_kkpipi[2]);pim.push_back(pim_kkpipi[3]);
2107 kp.push_back(kp_kkpipi[0]);kp.push_back(kp_kkpipi[1]);kp.push_back(kp_kkpipi[2]);kp.push_back(kp_kkpipi[3]);
2108 km.push_back(km_kkpipi[0]);km.push_back(km_kkpipi[1]);km.push_back(km_kkpipi[2]);km.push_back(km_kkpipi[3]);
2109
2110 pdau[0]=kp_kkpipi[0]; pdau[1]=kp_kkpipi[1]; pdau[2]=kp_kkpipi[2]; pdau[3]=kp_kkpipi[3];
2111 pdau[4]=km_kkpipi[0]; pdau[5]=km_kkpipi[1]; pdau[6]=km_kkpipi[2]; pdau[7]=km_kkpipi[3];
2112 pdau[8]=pip_kkpipi[0];pdau[9]=pip_kkpipi[1];pdau[10]=pip_kkpipi[2];pdau[11]=pip_kkpipi[3];
2113 pdau[12]=pim_kkpipi[0];pdau[13]=pim_kkpipi[1];pdau[14]=pim_kkpipi[2];pdau[15]=pim_kkpipi[3];
2114
2115 D0ToKKpipi tKKpipi;tKKpipi.init();
2116 D0 = tKKpipi.AMP(pdau,1);
2117
2118 pdau[0]=-km_kkpipi[0]; pdau[1]=-km_kkpipi[1]; pdau[2]=-km_kkpipi[2]; pdau[3]=km_kkpipi[3];
2119 pdau[4]=-kp_kkpipi[0]; pdau[5]=-kp_kkpipi[1]; pdau[6]=-kp_kkpipi[2]; pdau[7]=kp_kkpipi[3];
2120 pdau[8]=-pim_kkpipi[0]; pdau[9]=-pim_kkpipi[1]; pdau[10]=-pim_kkpipi[2];pdau[11]=pim_kkpipi[3];
2121 pdau[12]=-pip_kkpipi[0];pdau[13]=-pip_kkpipi[1];pdau[14]=-pip_kkpipi[2];pdau[15]=pip_kkpipi[3];
2122 D0bar = tKKpipi.AMP(pdau,1);
2123
2124 if( charm == 1){
2125
2126 r2D0 = std::norm(D0bar)/std::norm(D0);
2127 double temp = std::arg(D0)-std::arg(D0bar);
2128 log << MSG::INFO << "D0 calculation " << sqrt(r2D0) << " "<< temp << endmsg;
2129 while (temp < -TMath::Pi())
2130 {
2131 temp += 2.0 * TMath::Pi();
2132 }
2133 while (temp > TMath::Pi())
2134 {
2135 temp -= 2.0 * TMath::Pi();
2136 }
2137 deltaD0 = temp;
2138 }
2139 else {
2140 r2D0 = std::norm(D0)/std::norm(D0bar);
2141 double temp = std::arg(D0bar)-std::arg(D0);
2142 while (temp < -TMath::Pi())
2143 {
2144 temp += 2.0 * TMath::Pi();
2145 }
2146 while (temp > TMath::Pi())
2147 {
2148 temp -= 2.0 * TMath::Pi();
2149 }
2150 deltaD0 = temp;
2151 }
2152
2153 }else if( num[ kLeptonPlus ] == 1 )
2154 {
2155 decay_mode = kSLPlus ;
2156 ++m_nSL ;
2157
2158 r2D0 = 0. ;
2159 deltaD0 = 0. ;
2160 }
2161 else if( num[ kLeptonMinus ] == 1 )
2162 {
2163 decay_mode = kSLMinus ;
2164 ++m_nSL ;
2165
2166 r2D0 = 0. ;
2167 deltaD0 = 0. ;
2168 }
2169 //Check CP even modes
2170 else if(
2171 ( nDaughters == 2 &&
2172 ( ( num[ kPiPlus ] == 1 && num[ kPiMinus ] == 1 ) ||
2173 nUFP0 == 2 ||
2174 num[ kNeutralScalar ] == 2 ||
2175 ( nUFP0 == 1 && nUFV0 == 1 ) ||
2176 ( num[ kKPlus ] == 1 && num[ kKMinus ] == 1 ) ||
2177 num[ kK0S ] == 2 ||
2178 num[ kK0L ] == 2 ||
2179 ( num[ kK0L ] == 1 && nUFP0 == 1 ) ||
2180 ( num[ kK0S ] == 1 && num[ kNeutralScalar ] == 1 ) ||
2181 ( num[ kK0L ] == 1 && nUFV0 == 1 ) ||
2182 ( num[ kUFVPlus ] ==1 && num [ kPiMinus ] == 1 ) ||
2183 ( num[ kUFVMinus ] ==1 && num [ kPiPlus ] == 1 ) ) ) ||
2184 ( nDaughters == 3 &&
2185 ( ( num[ kPiPlus ] == 1 && num[ kPiMinus ] == 1 && num[ kPi0 ] == 1 ) ||
2186 ( nCPPlusEig == 1 && num[ kPi0 ] == 2 ) ||
2187 ( num[ kK0S ] == 3 ) || ( num[ kK0S ] == 1 && num[ kK0L ] == 2 ) ||
2188 ( num[ kK0S ] == 1 && num[ kK0L ] == 1 && nUFP0 == 1 ) ||
2189 ( num[ kK0S ] == 1 && nUFP0 == 2 )
2190 ) ) ||
2191 ( nDaughters == 4 &&
2192 ( num[ kK0L ] == 1 && nUFP0 == 3 )||
2193 ( ( num[ kK0S ] == 2 || num[ kK0L ] == 2 ) && nUFP0 == 2 ) )
2194 )
2195 {
2196 decay_mode = kCPPlus ;
2197 ++m_nCPPlus ;
2198
2199 r2D0 = 1. ;
2200 deltaD0 = 0.;
2201 }
2202 else if(( num[ kPiPlus ] == 1 && num[ kPiMinus ] == 1 && num[ kPi0 ] == 1 ) ){
2203
2204 decay_mode = kCPPlus ;
2205 ++m_nCPPlus ;
2206
2207 FCP = 0.973;
2208 r2D0 = 1. ;
2209 deltaD0 = 0.;
2210 }
2211 //Check CP odd modes
2212 else if(
2213 ( nDaughters == 2 &&
2214 ( ( num[ kK0S ] == 1 && nUFP0 == 1 ) ||
2215 ( num[ kK0L ] == 1 && num[ kNeutralScalar ] == 1 ) ||
2216 ( num[ kK0S ] == 1 && nUFV0 == 1 ) ||
2217 ( num[ kNeutralScalar ] == 1 && nUFV0 == 1 ) ||
2218 ( nUFP0 == 1 && num[ kNeutralScalar ] == 1 ) ) ) ||
2219 ( nDaughters == 3 &&
2220 ( ( nCPMinusEig == 3 && num[ kPi0 ] == 2 ) || // pi0 is subset of CP-
2221 ( num[ kPi0 ] == 3 ) ||
2222 ( num[ kK0L ] == 3 ) || ( num[ kK0S ] == 2 && num[ kK0L ] == 1 ) ||
2223 ( ( num[ kK0S ] == 2 || num[ kK0L ] == 2 ) && nUFP0 == 1 ) ||
2224 ( num[ kK0L ] == 1 && nUFP0 == 2 ) ) ) ||
2225 ( nDaughters == 4 &&
2226 ( num[ kK0S ] == 1 && num[ kPi0 ] == 3 )|| ( ( num[ kK0S ] == 1 && num[ kK0L ] == 1 ) && nUFP0 == 2 ) )
2227 )
2228 {
2229 decay_mode = kCPMinus ;
2230 ++m_nCPMinus ;
2231
2232 r2D0 = 1. ;
2233 deltaD0 = PI ;
2234 }
2235 else if( nStrange == nAntiStrange + 1 )
2236 {
2237 if(m_debug)cout<< nDaughters <<endl;
2238 if(m_debug)cout<< num[ kKMinus ]<<endl;
2239 if(m_debug)cout<< num[ kPiMinus ]<<endl;
2240 if(m_debug)cout<< num[ kPiPlus ]<<endl;
2241 if( charm == 1 )
2242 {
2243 if( ( num[ kKMinus ] == 1 && num[ kPiPlus ] == 2 && num[ kPiMinus ] == 1 ) && nDaughters == 4 ){
2244 if(m_debug)cout<< "True K-3Pi "<<endl;
2245 decay_mode = kFlavoredCF_3 ;
2246 ++m_nFlavoredCFD0_3 ;
2247
2248 complex<double> D0, D0bar;
2249 vector<double> kp; vector<double> pi1; vector<double> pi2;vector<double> pi3;
2250 double pdau[16];
2251 pi1.clear();pi2.clear();pi3.clear();
2252 kp.clear();
2253
2254 HepLorentzVector kp_kpipipi = (kMinus_kkpipi.at(0));
2255 HepLorentzVector pi1_kpipipi = (piPlus_4pi.at(0));
2256 HepLorentzVector pi2_kpipipi = (piPlus_4pi.at(1));
2257 HepLorentzVector pi3_kpipipi = (piMinus_4pi.at(0));
2258
2259 kp.push_back( kp_kpipipi[0]); kp.push_back(kp_kpipipi[1]); kp.push_back(kp_kpipipi[2]); kp.push_back(kp_kpipipi[3]); //k-
2260 pi1.push_back(pi1_kpipipi[0]);pi1.push_back(pi1_kpipipi[1]);pi1.push_back(pi1_kpipipi[2]);pi1.push_back(pi1_kpipipi[3]);//pi+
2261 pi2.push_back(pi2_kpipipi[0]);pi2.push_back(pi2_kpipipi[1]);pi2.push_back(pi2_kpipipi[2]);pi2.push_back(pi2_kpipipi[3]);//pi+
2262 pi3.push_back(pi3_kpipipi[0]);pi3.push_back(pi3_kpipipi[1]);pi3.push_back(pi3_kpipipi[2]);pi3.push_back(pi3_kpipipi[3]);//pi-
2263
2264 pdau[0]=kp[0]; pdau[1]=kp[1]; pdau[2]=kp[2] ; pdau[3]=kp[3];
2265 pdau[4]=pi1[0]; pdau[5]=pi1[1]; pdau[6]=pi1[2]; pdau[7]=pi1[3];
2266 pdau[8]=pi2[0]; pdau[9]=pi2[1]; pdau[10]=pi2[2];pdau[11]=pi2[3];
2267 pdau[12]=pi3[0];pdau[13]=pi3[1];pdau[14]=pi3[2];pdau[15]=pi3[3];
2268
2269 D0ToKpipipiLHCb tKpipipi;tKpipipi.init();
2270 D0 = tKpipipi.AMP(pdau,1);
2271
2272 D0TopiKpipi tpiKpipi;tpiKpipi.init();
2273 std::complex<double> dcs_offset = 0.0601387 * exp( std::complex<double>(0,1) * 1.04827 * M_PI / 180. );
2274 D0bar = dcs_offset*tpiKpipi.AMP(pdau,1);
2275 deltaD0 = ( std::arg(D0 * std::conj(D0bar)) + m_dCF_3 );
2276 r2D0 = std::norm(D0bar)/std::norm(D0);
2277
2278 }else if( ( num[ kKMinus ] == 1 && num[ kPiPlus ] == 1 && num[ kPi0 ] == 1 ) && nDaughters == 3 ){
2279 decay_mode = kFlavoredCF_1 ;
2280 ++m_nFlavoredCFD0_1 ;
2281 r2D0 = pow(m_rCF_1,2) ;
2282 deltaD0 = m_deltaCF_1 ;
2283 RD0 = m_RCF_1;
2284 }else{
2285 //if( ( num[ kKMinus ] == 1 && num[ kPiMinus ] == 1 ) && nDaughters == 2 )
2286 if(m_debug)cout<< "True K-Pi "<<endl;
2287 decay_mode = kFlavoredCF_0 ;
2288 ++m_nFlavoredCFD0_0 ;
2289 r2D0 = pow(m_rCF_0,2) ;
2290 deltaD0 = m_deltaCF_0 ;
2291 RD0 = m_RCF_0;
2292 }
2293 }
2294 else
2295 {
2296 if( ( num[ kKMinus ] == 1 && num[ kPiPlus ] == 2 && num[ kPiMinus ] == 1 ) && nDaughters == 4 ){
2297 decay_mode = kFlavoredCF_3 ;
2298 ++m_nFlavoredDCSD0_3 ;
2299
2300 complex<double> D0, D0bar;
2301 vector<double> kp; vector<double> pi1; vector<double> pi2;vector<double> pi3;
2302 double pdau[16];
2303 pi1.clear();pi2.clear();pi3.clear();
2304 kp.clear();
2305
2306 HepLorentzVector kp_kpipipi = (kMinus_kkpipi.at(0));
2307 HepLorentzVector pi1_kpipipi = (piPlus_4pi.at(0));
2308 HepLorentzVector pi2_kpipipi = (piPlus_4pi.at(1));
2309 HepLorentzVector pi3_kpipipi = (piMinus_4pi.at(0));
2310
2311 kp.push_back( kp_kpipipi[0]); kp.push_back(kp_kpipipi[1]); kp.push_back(kp_kpipipi[2]); kp.push_back(kp_kpipipi[3]); //k-
2312 pi1.push_back(pi1_kpipipi[0]);pi1.push_back(pi1_kpipipi[1]);pi1.push_back(pi1_kpipipi[2]);pi1.push_back(pi1_kpipipi[3]);//pi+
2313 pi2.push_back(pi2_kpipipi[0]);pi2.push_back(pi2_kpipipi[1]);pi2.push_back(pi2_kpipipi[2]);pi2.push_back(pi2_kpipipi[3]);//pi+
2314 pi3.push_back(pi3_kpipipi[0]);pi3.push_back(pi3_kpipipi[1]);pi3.push_back(pi3_kpipipi[2]);pi3.push_back(pi3_kpipipi[3]);//pi-
2315
2316 pdau[0]=kp[0]; pdau[1]=kp[1]; pdau[2]=kp[2] ; pdau[3]=kp[3];
2317 pdau[4]=pi1[0]; pdau[5]=pi1[1]; pdau[6]=pi1[2]; pdau[7]=pi1[3];
2318 pdau[8]=pi2[0]; pdau[9]=pi2[1]; pdau[10]=pi2[2];pdau[11]=pi2[3];
2319 pdau[12]=pi3[0];pdau[13]=pi3[1];pdau[14]=pi3[2];pdau[15]=pi3[3];
2320
2321 D0ToKpipipiLHCb tKpipipi;tKpipipi.init();
2322 D0 = tKpipipi.AMP(pdau,1);
2323
2324 D0TopiKpipi tpiKpipi;tpiKpipi.init();
2325 std::complex<double> dcs_offset = 0.0601387 * exp( std::complex<double>(0,1) * 1.04827 * M_PI / 180. );
2326 D0bar = dcs_offset*tpiKpipi.AMP(pdau,1);
2327 deltaD0 = - ( std::arg(D0 * std::conj(D0bar)) + m_dCF_3 );
2328 r2D0 = std::norm(D0)/std::norm(D0bar);
2329
2330 }else if( ( num[ kKMinus ] == 1 && num[ kPiPlus ] == 1 && num[ kPi0 ] == 1 ) && nDaughters == 3 ){
2331 decay_mode = kFlavoredCF_1 ;
2332 ++m_nFlavoredDCSD0_1 ;
2333 r2D0 = 1. / pow( m_rCF_1,2 ) ;
2334 deltaD0 = -m_deltaCF_1 ;
2335 RD0 = m_RCF_1;
2336 }else{
2337 //if( ( num[ kKMinus ] == 1 && num[ kPiPlus ] == 1 ) && nDaughters == 2 )
2338 decay_mode = kFlavoredCF_0 ;
2339 ++m_nFlavoredDCSD0_0 ;
2340 r2D0 = 1. / pow( m_rCF_0,2 ) ;
2341 deltaD0 = -m_deltaCF_0 ;
2342 RD0 = m_RCF_0;
2343 }
2344 }
2345 }
2346 else if( nAntiStrange == nStrange + 1 )
2347 {
2348 if( charm == -1 )
2349 {
2350 if( ( num[ kKPlus ] == 1 && num[ kPiMinus ] == 2 && num[ kPiPlus ] == 1 ) && nDaughters == 4 ){
2351 decay_mode = kFlavoredCFBar_3 ;
2352 ++m_nFlavoredCFD0_3 ;
2353
2354 complex<double> D0, D0bar;
2355 vector<double> kp; vector<double> pi1; vector<double> pi2;vector<double> pi3;
2356 double pdau[16];
2357 pi1.clear();pi2.clear();pi3.clear();
2358 kp.clear();
2359
2360 HepLorentzVector kp_kpipipi = (kPlus_kkpipi.at(0));
2361 HepLorentzVector pi1_kpipipi = (piMinus_4pi.at(0));
2362 HepLorentzVector pi2_kpipipi = (piMinus_4pi.at(1));
2363 HepLorentzVector pi3_kpipipi = (piPlus_4pi.at(0));
2364
2365 kp.push_back( -kp_kpipipi[0]); kp.push_back(-kp_kpipipi[1]); kp.push_back(-kp_kpipipi[2]); kp.push_back(kp_kpipipi[3]); //k-
2366 pi1.push_back(-pi1_kpipipi[0]);pi1.push_back(-pi1_kpipipi[1]);pi1.push_back(-pi1_kpipipi[2]);pi1.push_back(pi1_kpipipi[3]);//pi+
2367 pi2.push_back(-pi2_kpipipi[0]);pi2.push_back(-pi2_kpipipi[1]);pi2.push_back(-pi2_kpipipi[2]);pi2.push_back(pi2_kpipipi[3]);//pi+
2368 pi3.push_back(-pi3_kpipipi[0]);pi3.push_back(-pi3_kpipipi[1]);pi3.push_back(-pi3_kpipipi[2]);pi3.push_back(pi3_kpipipi[3]);//pi-
2369
2370 pdau[0]=kp[0]; pdau[1]=kp[1]; pdau[2]=kp[2] ; pdau[3]=kp[3];
2371 pdau[4]=pi1[0]; pdau[5]=pi1[1]; pdau[6]=pi1[2]; pdau[7]=pi1[3];
2372 pdau[8]=pi2[0]; pdau[9]=pi2[1]; pdau[10]=pi2[2];pdau[11]=pi2[3];
2373 pdau[12]=pi3[0];pdau[13]=pi3[1];pdau[14]=pi3[2];pdau[15]=pi3[3];
2374
2375 D0ToKpipipiLHCb tKpipipi;tKpipipi.init();
2376 D0 = tKpipipi.AMP(pdau,1);
2377
2378 D0TopiKpipi tpiKpipi;tpiKpipi.init();
2379 std::complex<double> dcs_offset = 0.0601387 * exp( std::complex<double>(0,1) * 1.04827 * M_PI / 180. );
2380 D0bar = dcs_offset*tpiKpipi.AMP(pdau,1);
2381 deltaD0 = ( std::arg(D0 * std::conj(D0bar)) + m_dCF_3 );
2382 r2D0 = std::norm(D0bar)/std::norm(D0);
2383
2384
2385 }else if( ( num[ kKPlus ] == 1 && num[ kPiMinus ] == 1 && num[ kPi0 ] == 1 ) && nDaughters == 3 ){
2386 decay_mode = kFlavoredCFBar_1 ;
2387 ++m_nFlavoredCFD0_1 ;
2388 r2D0 = pow(m_rCF_1,2) ;
2389 deltaD0 = m_deltaCF_1 ;
2390 RD0 = m_RCF_1;
2391 }else{
2392 //else if( ( num[ kKPlus ] == 1 && num[ kPiMinus ] == 1 ) && nDaughters == 2 )
2393 decay_mode = kFlavoredCFBar_0 ;
2394 ++m_nFlavoredCFD0_0 ;
2395 r2D0 = pow(m_rCF_0,2) ;
2396 deltaD0 = m_deltaCF_0 ;
2397 RD0 = m_RCF_0;
2398 }
2399 }
2400 else
2401 {
2402 if( ( num[ kKPlus ] == 1 && num[ kPiMinus ] == 2 && num[ kPiPlus ] == 1 ) && nDaughters == 4 ){
2403 decay_mode = kFlavoredCFBar_3 ;
2404 ++m_nFlavoredDCSD0_3 ;
2405
2406 complex<double> D0, D0bar;
2407 vector<double> kp; vector<double> pi1; vector<double> pi2;vector<double> pi3;
2408 double pdau[16];
2409 pi1.clear();pi2.clear();pi3.clear();
2410 kp.clear();
2411
2412 HepLorentzVector kp_kpipipi = (kPlus_kkpipi.at(0));
2413 HepLorentzVector pi1_kpipipi = (piMinus_4pi.at(0));
2414 HepLorentzVector pi2_kpipipi = (piMinus_4pi.at(1));
2415 HepLorentzVector pi3_kpipipi = (piPlus_4pi.at(0));
2416
2417 kp.push_back( -kp_kpipipi[0]); kp.push_back(-kp_kpipipi[1]); kp.push_back(-kp_kpipipi[2]); kp.push_back(kp_kpipipi[3]);//k-
2418 pi1.push_back(-pi1_kpipipi[0]);pi1.push_back(-pi1_kpipipi[1]);pi1.push_back(-pi1_kpipipi[2]);pi1.push_back(pi1_kpipipi[3]);//pi+
2419 pi2.push_back(-pi2_kpipipi[0]);pi2.push_back(-pi2_kpipipi[1]);pi2.push_back(-pi2_kpipipi[2]);pi2.push_back(pi2_kpipipi[3]);//pi+
2420 pi3.push_back(-pi3_kpipipi[0]);pi3.push_back(-pi3_kpipipi[1]);pi3.push_back(-pi3_kpipipi[2]);pi3.push_back(pi3_kpipipi[3]);//pi-
2421
2422 pdau[0]=kp[0]; pdau[1]=kp[1]; pdau[2]=kp[2] ; pdau[3]=kp[3];
2423 pdau[4]=pi1[0]; pdau[5]=pi1[1]; pdau[6]=pi1[2]; pdau[7]=pi1[3];
2424 pdau[8]=pi2[0]; pdau[9]=pi2[1]; pdau[10]=pi2[2];pdau[11]=pi2[3];
2425 pdau[12]=pi3[0];pdau[13]=pi3[1];pdau[14]=pi3[2];pdau[15]=pi3[3];
2426
2427 D0ToKpipipiLHCb tKpipipi;tKpipipi.init();
2428 D0 = tKpipipi.AMP(pdau,1);
2429
2430 D0TopiKpipi tpiKpipi;tpiKpipi.init();
2431 std::complex<double> dcs_offset = 0.0601387 * exp( std::complex<double>(0,1) * 1.04827 * M_PI / 180. );
2432 D0bar = dcs_offset*tpiKpipi.AMP(pdau,1);
2433 deltaD0 = - ( std::arg(D0 * std::conj(D0bar)) + m_dCF_3 );
2434 r2D0 = std::norm(D0)/std::norm(D0bar);
2435
2436 }else if( ( num[ kKPlus ] == 1 && num[ kPiMinus ] == 1 && num[ kPi0 ] == 1 ) && nDaughters == 3 ){
2437 decay_mode = kFlavoredCFBar_1 ;
2438 ++m_nFlavoredDCSD0_1 ;
2439 r2D0 = 1. / pow( m_rCF_1,2 ) ;
2440 deltaD0 = -m_deltaCF_1 ;
2441 RD0 = m_RCF_1;
2442 }else{
2443 //if( ( num[ kKMinus ] == 1 && num[ kPiMinus ] == 1 ) && nDaughters == 2 )
2444 decay_mode = kFlavoredCFBar_0 ;
2445 ++m_nFlavoredDCSD0_0 ;
2446 r2D0 = 1. / pow( m_rCF_0,2 ) ;
2447 deltaD0 = -m_deltaCF_0 ;
2448 RD0 = m_RCF_0;
2449 }
2450 }
2451 }
2452 else if( nStrange == nAntiStrange ) // Unflavored or Cabibbo-suppressed
2453 {
2454 if( ( num[ kKPlus ] > 0 &&
2455 ( num[ kKPlus ] == num[ kFVMinus ] ||
2456 num[ kKPlus ] == nFV0Bar ) ) || // for anti-K*0 K+ pi-
2457 ( nK0 > 0 &&
2458 nFV0 != nFV0Bar &&
2459 nK0 == nFV0Bar ) ||
2460 ( num[ kPiPlus ] > 0 &&
2461 num[ kPiPlus ] == num[ kUFVMinus ] ) )
2462 {
2463 decay_mode = kFlavoredCS ;
2464 ++m_nFlavoredCSD0 ;
2465
2466 if( charm == 1 )
2467 {
2468 r2D0 = pow( m_rCS, 2 ) ;
2469 deltaD0 = m_deltaCS ;
2470 }
2471 else
2472 {
2473 r2D0 = 1. / pow( m_rCS, 2 ) ;
2474 deltaD0 = -m_deltaCS ;
2475 }
2476 }
2477 else if( ( num[ kKMinus ] > 0 && ( num[ kKMinus ] == num[ kFVPlus ] || num[ kKMinus ] == nFV0 ) ) || // for K*0 K- pi+
2478 ( nK0 > 0 && nFV0 != nFV0Bar && nK0 == nFV0 ) ||
2479 ( num[ kPiMinus ] > 0 && num[ kPiMinus ] == num[ kUFVPlus ] ) )
2480 {
2481 decay_mode = kFlavoredCSBar ;
2482 ++m_nFlavoredCSD0 ;
2483
2484 if( charm == -1 )
2485 {
2486 r2D0 = pow( m_rCS, 2 ) ;
2487 deltaD0 = m_deltaCS ;
2488 }
2489 else
2490 {
2491 r2D0 = 1. / pow( m_rCS, 2 ) ;
2492 deltaD0 = -m_deltaCS ;
2493 }
2494 }
2495
2496 // Self-conjugate mixed-CP final states -- pick CP or non-CP at
2497 // random. If CP, pick + or - at random. If non-CP, pick
2498 // flavored or flavoredBar according to the appropriate value of r.
2499 else if( ( nDaughters >= 3 && num[ kPiPlus ] == num[ kPiMinus ] &&
2500 num[ kKPlus ] == num[ kKMinus ] && nChargedPiK + nCPEig == nDaughters ) ||
2501 nUFV0 == nDaughters ||
2502 ( num[ kKStar0 ] > 0 && num[ kKStar0 ] == num[ kKStar0Bar ] ) ||
2503 ( num[ kK10 ] > 0 && num[ kK10 ] == num[ kK10Bar ] ) ||
2504 ( num[ kUFVPlus ] == 1 && num[ kUFVMinus ] == 1 )
2505 // for rho+rho-; no a1+a1- and rhoa1 final states in DECAY.DEC
2506 )
2507 {
2508 log << MSG::DEBUG <<" [ Self-conjugate mixed-CP ]" << endmsg ;
2509
2510 if( RandFlat::shoot() > 0.5 )
2511 {
2512 if( RandFlat::shoot() > 0.5 )
2513 {
2514 decay_mode = kCPPlus ;
2515 ++m_nCPPlus ;
2516
2517 r2D0 = 1. ;
2518 deltaD0 = PI;
2519 }
2520 else
2521 {
2522 decay_mode = kCPMinus ;
2523 ++m_nCPMinus ;
2524
2525 r2D0 = 1. ;
2526 deltaD0 = 0. ;
2527 }
2528 }
2529 else
2530 {
2531 // Odd # of K0S or K0L = CF; even # of K0S or K0L = SC.
2532 if( nK0 % 2 )
2533 {
2534 // Nov 2007: correct for r2 -> RWS and BR != A^2
2535 double cutoff = m_rwsCF_0 / ( 1. + m_rwsCF_0 ) ;
2536
2537 if( charm == -1 )
2538 {
2539 cutoff = 1. - cutoff ;
2540 }
2541
2542 log << MSG::DEBUG <<" [ cutoff = " << cutoff << " ]" << endmsg ;
2543
2544 decay_mode = ( RandFlat::shoot() > cutoff ) ? kFlavoredCF_0 : kFlavoredCFBar_0 ;
2545
2546 if( ( decay_mode == kFlavoredCF_0 && charm == 1 ) ||
2547 ( decay_mode == kFlavoredCFBar_0 && charm == -1 ) )
2548 {
2549 ++m_nFlavoredCFD0_0 ;
2550
2551 r2D0 = sqrt( m_rCF_0 ) ;
2552 deltaD0 = m_deltaCF_0 ;
2553 }
2554 else
2555 {
2556 ++m_nFlavoredDCSD0_0 ;
2557
2558 r2D0 = 1. / sqrt( m_rCF_0 ) ;
2559 deltaD0 = -m_deltaCF_0 ;
2560 }
2561 }
2562 else
2563 {
2564 // Nov 2007: correct for r2 -> RWS and BR != A^2
2565 double cutoff = m_rwsCS / ( 1. + m_rwsCS ) ;
2566
2567 if( charm == -1 )
2568 {
2569 cutoff = 1. - cutoff ;
2570 }
2571
2572 log << MSG::DEBUG <<" [ cutoff = " << cutoff << " ]" << endmsg ;
2573
2574 decay_mode = ( RandFlat::shoot() > cutoff ) ?
2575 kFlavoredCS : kFlavoredCSBar ;
2576 ++m_nFlavoredCSD0 ;
2577
2578 if( ( decay_mode == kFlavoredCS && charm == 1 ) ||
2579 ( decay_mode == kFlavoredCSBar && charm == -1 ) )
2580 {
2581 r2D0 = sqrt( m_rCS ) ;
2582 deltaD0 = m_deltaCS ;
2583 }
2584 else
2585 {
2586 r2D0 = 1. / sqrt( m_rCS ) ;
2587 deltaD0 = -m_deltaCS ;
2588 }
2589 }
2590 }
2591 }
2592 }
2593
2594 if (num[kUnknown] >= 1)
2595 {
2596 cout << "decay mode " << decay_mode << endl ;
2597 cout << "D #" << charm << endl ;
2598 cout << "pi+: " << num[ kPiPlus ] << endl ;
2599 cout << "pi-: " << num[ kPiMinus ] << endl ;
2600 cout << "pi0: " << num[ kPi0 ] << endl ;
2601 cout << "eta: " << num[ kEta ] << endl ;
2602 cout << "eta': " << num[ kEtaPrime ] << endl ;
2603 cout << "f0/a0: " << num[ kNeutralScalar ] << endl ;
2604 cout << "UFV+: " << num[ kUFVPlus ] << endl ;
2605 cout << "UFV-: " << num[ kUFVMinus ] << endl ;
2606 cout << "rho0: " << num[ kRho0 ] << endl ;
2607 cout << "omega: " << num[ kOmega ] << endl ;
2608 cout << "phi: " << num[ kPhi ] << endl ;
2609 cout << "K+: " << num[ kKPlus ] << endl ;
2610 cout << "K-: " << num[ kKMinus ] << endl ;
2611 cout << "K0S: " << num[ kK0S ] << endl ;
2612 cout << "K0L: " << num[ kK0L ] << endl ;
2613 cout << "FV+: " << num[ kFVPlus ] << endl ;
2614 cout << "FV-: " << num[ kFVMinus ] << endl ;
2615 cout << "K*0: " << num[ kKStar0 ] << endl ;
2616 cout << "K*0b: " << num[ kKStar0Bar ] << endl ;
2617 cout << "K10: " << num[ kK10 ] << endl ;
2618 cout << "K10b: " << num[ kK10Bar ] << endl ;
2619 cout << "l+: " << num[ kLeptonPlus ] << endl ;
2620 cout << "l-: " << num[ kLeptonMinus ] << endl ;
2621 cout << "nu: " << num[ kNeutrino ] << endl ;
2622 cout << "gamma: " << num[ kGamma ] << endl ;
2623 cout << "Unknown: " << num[ 25 ] << endl ;
2624 }
2625
2626
2627 d0_info[0]=decay_mode;
2628 d0_info[1]=r2D0;
2629 d0_info[2]=deltaD0;
2630 d0_info[3]=RD0;
2631 d0_info[4]=FCP;
2632 d0_info[5]=mnm_gen;
2633 d0_info[6]=mpn_gen;
2634 d0_info[7]= double (isKsPiPi);
2635 d0_info[8]=mpm_gen;
2636 return d0_info;
2637}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
runNo
int runNo
Definition: DQA_TO_DB.cxx:12
TComplex
complex< double > TComplex
Definition: Dalitz.h:14
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition: Evt3Rank3C.cc:175
imag
double imag(const EvtComplex &c)
Definition: EvtComplex.hh:246
exp
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
bin
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition: FoamA.h:85
s
XmlRpcServer s
Definition: HelloServer.cpp:11
weight
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
Definition: KarFin.h:34
m_x
NTuple::Array< double > m_x
Definition: MdcHistItem.h:112
m_y
NTuple::Array< double > m_y
Definition: MdcHistItem.h:113
PI
const double PI
Definition: PipiJpsi.cxx:55
Pip2Pim2Denom_fil
double Pip2Pim2Denom_fil
Definition: QCMCFilter.cxx:203
Pip2Pim2Numer1_fil
double Pip2Pim2Numer1_fil
Definition: QCMCFilter.cxx:201
m_rwsCS
double m_rwsCS
Definition: QCMCFilter.cxx:221
xmpion
const double xmpion
Definition: QCMCFilter.cxx:71
m_PipPim2Pi0Numer2
double m_PipPim2Pi0Numer2
Definition: QCMCFilter.cxx:205
m_rwsCF_3
double m_rwsCF_3
Definition: QCMCFilter.cxx:220
m_nCPPlus
int m_nCPPlus
Definition: QCMCFilter.cxx:177
dalitzNumer2_fil
double dalitzNumer2_fil
Definition: QCMCFilter.cxx:196
m_n2Pip2Pim
int m_n2Pip2Pim
Definition: QCMCFilter.cxx:188
xmk0
const double xmk0
Definition: QCMCFilter.cxx:72
m_nCPMinus
int m_nCPMinus
Definition: QCMCFilter.cxx:178
m_deltaCF_1
double m_deltaCF_1
Definition: QCMCFilter.cxx:224
m_deltaCF
double m_deltaCF
Definition: QCMCFilter.cxx:222
m_nUnknownEvents
int m_nUnknownEvents
Definition: QCMCFilter.cxx:170
m_dalitzNumer1
double m_dalitzNumer1
Definition: QCMCFilter.cxx:192
m_2Pip2PimNumer1
double m_2Pip2PimNumer1
Definition: QCMCFilter.cxx:198
m_modeCounter
HepMatrix m_modeCounter(19, 19, 0)
m_deltaCF_0
double m_deltaCF_0
Definition: QCMCFilter.cxx:223
m_rwsCF
double m_rwsCF
Definition: QCMCFilter.cxx:217
m_nSL
int m_nSL
Definition: QCMCFilter.cxx:186
m_rwsCF_0
double m_rwsCF_0
Definition: QCMCFilter.cxx:218
m_nKKPiPi
int m_nKKPiPi
Definition: QCMCFilter.cxx:191
m_nDalitz
int m_nDalitz
Definition: QCMCFilter.cxx:187
PipPim2Pi0Denom_fil
double PipPim2Pi0Denom_fil
Definition: QCMCFilter.cxx:209
m_PipPim2Pi0Numer1
double m_PipPim2Pi0Numer1
Definition: QCMCFilter.cxx:204
m_nD0D0barDiscarded
int m_nD0D0barDiscarded
Definition: QCMCFilter.cxx:175
m_nDpDmEvents
int m_nDpDmEvents
Definition: QCMCFilter.cxx:174
PipPim2Pi0Numer1_fil
double PipPim2Pi0Numer1_fil
Definition: QCMCFilter.cxx:207
Pip2Pim2Numer2_fil
double Pip2Pim2Numer2_fil
Definition: QCMCFilter.cxx:202
m_dalitzDenom
double m_dalitzDenom
Definition: QCMCFilter.cxx:194
m_weights
HepSymMatrix m_weights(18, 0)
KSPiPiPi0Denom_fil
double KSPiPiPi0Denom_fil
Definition: QCMCFilter.cxx:215
m_nK0PiPiPi0
int m_nK0PiPiPi0
Definition: QCMCFilter.cxx:190
xmd0
const double xmd0
Definition: QCMCFilter.cxx:74
m_rho_flag
int m_rho_flag
Definition: QCMCFilter.cxx:230
m_nFlavoredCFD0_3
int m_nFlavoredCFD0_3
Definition: QCMCFilter.cxx:181
xmeta
const double xmeta
Definition: QCMCFilter.cxx:69
m_2Pip2PimNumer2
double m_2Pip2PimNumer2
Definition: QCMCFilter.cxx:199
m_PipPim2Pi0Denom
double m_PipPim2Pi0Denom
Definition: QCMCFilter.cxx:206
PI
const double PI
Definition: QCMCFilter.cxx:76
xmkaon
const double xmkaon
Definition: QCMCFilter.cxx:70
m_nD0bar
int m_nD0bar
Definition: QCMCFilter.cxx:173
m_nFlavoredCFD0_0
int m_nFlavoredCFD0_0
Definition: QCMCFilter.cxx:179
m_KSPiPiPi0Numer1
double m_KSPiPiPi0Numer1
Definition: QCMCFilter.cxx:210
KSPiPiPi0Numer1_fil
double KSPiPiPi0Numer1_fil
Definition: QCMCFilter.cxx:213
m_nDpDmDiscarded
int m_nDpDmDiscarded
Definition: QCMCFilter.cxx:176
PipPim2Pi0Numer2_fil
double PipPim2Pi0Numer2_fil
Definition: QCMCFilter.cxx:208
m_nFlavoredDCSD0_3
int m_nFlavoredDCSD0_3
Definition: QCMCFilter.cxx:185
m_nFlavoredDCSD0_0
int m_nFlavoredDCSD0_0
Definition: QCMCFilter.cxx:183
xmpi0
const double xmpi0
Definition: QCMCFilter.cxx:68
m_dalitzNumer2
double m_dalitzNumer2
Definition: QCMCFilter.cxx:193
m_KSPiPiPi0Numer2
double m_KSPiPiPi0Numer2
Definition: QCMCFilter.cxx:211
dalitzDenom_fil
double dalitzDenom_fil
Definition: QCMCFilter.cxx:197
m_rwsCF_1
double m_rwsCF_1
Definition: QCMCFilter.cxx:219
m_nFlavoredCSD0
int m_nFlavoredCSD0
Definition: QCMCFilter.cxx:182
m_nUnknownDecays
int m_nUnknownDecays
Definition: QCMCFilter.cxx:171
m_nFlavoredDCSD0_1
int m_nFlavoredDCSD0_1
Definition: QCMCFilter.cxx:184
dalitzNumer1_fil
double dalitzNumer1_fil
Definition: QCMCFilter.cxx:195
m_KSPiPiPi0Denom
double m_KSPiPiPi0Denom
Definition: QCMCFilter.cxx:212
m_nD0D0barEvents
int m_nD0D0barEvents
Definition: QCMCFilter.cxx:172
m_deltaCS
double m_deltaCS
Definition: QCMCFilter.cxx:226
m_deltaCF_3
double m_deltaCF_3
Definition: QCMCFilter.cxx:225
m_nPipPim2Pi0
int m_nPipPim2Pi0
Definition: QCMCFilter.cxx:189
m_nFlavoredCFD0_1
int m_nFlavoredCFD0_1
Definition: QCMCFilter.cxx:180
KSPiPiPi0Numer2_fil
double KSPiPiPi0Numer2_fil
Definition: QCMCFilter.cxx:214
xmrho
const double xmrho
Definition: QCMCFilter.cxx:73
m_2Pip2PimDenom
double m_2Pip2PimDenom
Definition: QCMCFilter.cxx:200
m_keptModeCounter
HepMatrix m_keptModeCounter(19, 19, 0)
xmpi0
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmpi0
Definition: RRes.h:32
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
M_PI
#define M_PI
Definition: TConstant.h:4
t
TTree * t
Definition: binning.cxx:23
D0To2pip2pim::Amp
complex< double > Amp(vector< double > Pip1, vector< double > Pim1, vector< double > Pip2, vector< double > Pim2)
Definition: D0To2pip2pim.cxx:1426
D0ToKKpipi::init
void init()
Definition: D0ToKKpipi.cxx:15
D0ToKKpipi::AMP
complex< double > AMP(double const *x0, const int &x1)
Definition: D0ToKKpipi.cxx:2656
D0ToKLpipi2018::Amp_PFT
complex< double > Amp_PFT(vector< double > k0l, vector< double > pip, vector< double > pim)
Definition: D0ToKLpipi2018.cxx:92
D0ToKSpipi2018::Amp_PFT
complex< double > Amp_PFT(vector< double > k0l, vector< double > pip, vector< double > pim)
Definition: D0ToKSpipi2018.cxx:94
D0ToKSpipipi0::Amp
complex< double > Amp(vector< double > ks0, vector< double > pip, vector< double > pim, vector< double > pi0)
Definition: D0ToKSpipipi0.cxx:1615
D0ToKpipipiLHCb::AMP
std::complex< double > AMP(const double *x0, const int &x1)
Definition: D0ToKpipipiLHCb.cxx:2976
D0TopiKpipi::AMP
complex< double > AMP(double const *x0, const int &x1)
Definition: D0TopiKpipi.cxx:1086
D0Topippim2pi0::Amp
complex< double > Amp(vector< double > Pip, vector< double > Pim, vector< double > Pi01, vector< double > Pi02)
Definition: D0Topippim2pi0.cxx:1384
Dalitz
Definition: Dalitz.h:30
Event::McParticle::mother
const McParticle & mother() const
access to the mother particle
Definition: McParticle.cxx:117
Event::McParticle::particleProperty
StdHepId particleProperty() const
Retrieve particle property.
Definition: McParticle.cxx:7
QCMCFilter::initialize
StatusCode initialize()
Definition: QCMCFilter.cxx:272
QCMCFilter::finalize
StatusCode finalize()
Definition: QCMCFilter.cxx:1174
QCMCFilter::findD0Decay
std::vector< double > findD0Decay(int charm)
Definition: QCMCFilter.cxx:1507
QCMCFilter::execute
StatusCode execute()
Definition: QCMCFilter.cxx:826
y
double y[1000]
Definition: draw_charge_space_coarse_II.cxx:8
m_debug
bool m_debug
Definition: draw_charge_space_coarse_II.cxx:3
EventModel::MC::McParticleCol
_EXTERN_ std::string McParticleCol
Definition: EventModel.h:41