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