BOSS 7.0.7
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
7#include "GaudiKernel/MsgStream.h"
8#include "GaudiKernel/AlgFactory.h"
9#include "GaudiKernel/ISvcLocator.h"
10#include "GaudiKernel/SmartDataPtr.h"
11#include "GaudiKernel/IDataProviderSvc.h"
12#include "GaudiKernel/PropertyMgr.h"
13#include "GaudiKernel/Bootstrap.h"
14#include "GaudiKernel/RegistryEntry.h"
15
16#include "TMath.h"
17
18#include <cmath>
19#include "HepPDT/ParticleDataTable.hh"
20#include "HepPDT/ParticleData.hh"
21//#include "PartPropSvc/PartPropSvc.h"
22#include "GaudiKernel/IPartPropSvc.h"
23
24// Monte-Carlo data
25#include "McTruth/McParticle.h"
26#include "McTruth/MdcMcHit.h"
27
28///////////////////
29#include "EventModel/EventModel.h"
30#include "EventModel/Event.h"
31#include "EventModel/EventHeader.h"
32
33#include "EvtRecEvent/EvtRecEvent.h"
34
35#include "CLHEP/Random/RandFlat.h"
36#include "CLHEP/Matrix/Vector.h"
37#include "CLHEP/Matrix/Matrix.h"
38#include "CLHEP/Matrix/SymMatrix.h"
39#include "CLHEP/Vector/ThreeVector.h"
40#include "CLHEP/Vector/LorentzVector.h"
41#include "CLHEP/Vector/TwoVector.h"
42using CLHEP::HepVector;
43using CLHEP::Hep3Vector;
44using CLHEP::Hep2Vector;
45using CLHEP::HepLorentzVector;
46
47#include <vector>
48
49//
50// constants, enums and typedefs
51//
52// Partilcle masses
53const double xmpi0 = 0.134976; // BOSS 6.4.1 pdt.table value
54const double xmeta = 0.54784; // PDG08 = 0.547853, BOSS 6.4.1 pdt.table = 0.54784
55const double xmkaon = 0.49368;
56const double xmpion = 0.13957;
57const double xmk0 = 0.49761;
58const double xmrho = 0.77549;
59const double xmd0 = 1.86484; // PDG08
60
61const double PI = 3.1415926; //pi
62
63// Antiparticles have negative ID.
64
65static const int kPsi3770ID = 30443 ;
66static const int kD0ID = 421 ;
67static const int kD0BarID = -421 ;
68static const int kDpID = 411 ;
69static const int kDmID = -411 ;
70static const int kGammaID = 22 ;
71static const int kGammaFSRID = -22 ;
72static const int kPiPlusID = 211 ;
73static const int kPiMinusID = -211 ;
74static const int kPi0ID = 111 ;
75static const int kEtaID = 221 ;
76static const int kEtaPrimeID = 331 ;
77static const int kF0ID = 9010221;
78static const int kFPrime0ID = 10221 ;
79static const int kF0m1500ID = 50221;
80static const int kF2ID = 225;
81static const int kA00ID = 10111 ;
82static const int kA0PlusID = 10211 ;
83static const int kA0MinusID = -10211 ;
84static const int kRhoPlusID = 213 ;
85static const int kRhoMinusID = -213 ;
86static const int kRho0ID = 113 ;
87static const int kRho2SPlusID = 30213 ;
88static const int kRho2SMinusID = -30213 ;
89static const int kRho2S0ID = 30113 ;
90static const int kA1PlusID = 20213 ;
91static const int kA1MinusID = -20213 ;
92static const int kA10ID = 20113 ;
93static const int kOmegaID = 223 ;
94static const int kPhiID = 333 ;
95static const int kKPlusID = 321 ;
96static const int kKMinusID = -321 ;
97static const int kK0SID = 310 ;
98static const int kK0LID = 130 ;
99static const int kK0ID = 311 ;
100static const int kK0BarID = -311 ;
101static const int kKStarPlusID = 323 ;
102static const int kKStarMinusID = -323 ;
103static const int kKStar0ID = 313 ;
104static const int kKStar0BarID = -313 ;
105static const int kK0Star0ID = 10311;
106static const int kK0Star0BarID = -10311;
107static const int kK0StarPlusID = 10321;
108static const int kK0StarMinusID = -10321;
109static const int kK1PlusID = 10323 ;
110static const int kK1MinusID = -10323 ;
111static const int kK10ID = 10313 ;
112static const int kK10BarID = -10313 ;
113static const int kK1PrimePlusID = 20323 ;
114static const int kK1PrimeMinusID = -20323 ;
115static const int kK1Prime0ID = 20313 ;
116static const int kK1Prime0BarID = -20313 ;
117static const int kK2StarPlusID = 325 ;
118static const int kK2StarMinusID = -325 ;
119static const int kK2Star0ID = 315 ;
120static const int kK2Star0BarID = -315 ;
121static const int kEMinusID = 11 ;
122static const int kEPlusID = -11 ;
123static const int kMuMinusID = 13 ;
124static const int kMuPlusID = -13 ;
125static const int kNuEID = 12 ;
126static const int kNuEBarID = -12 ;
127static const int kNuMuID = 14 ;
128static const int kNuMuBarID = -14 ;
129
130//D0 decay modes
131static const int kFlavoredCF = 0;
132static const int kFlavoredCFBar = 1;
133static const int kFlavoredCS = 2;
134static const int kFlavoredCSBar = 3;
135static const int kSLPlus = 4;
136static const int kSLMinus = 5;
137static const int kCPPlus = 6;
138static const int kCPMinus = 7;
139static const int kDalitz = 8;
140static const int kNDecayTypes = 9;
141
142
143
144//Varibales to keep track of the Events
145int m_nUnknownEvents = 0;
146int m_nUnknownDecays= 0;
147int m_nD0D0barEvents= 0;
148int m_nD0bar = 0;
149int m_nDpDmEvents = 0;
150int m_nD0D0barDiscarded = 0;
151int m_nDpDmDiscarded = 0;
152int m_nCPPlus = 0;
153int m_nCPMinus = 0;
154int m_nFlavoredCFD0 = 0;
155int m_nFlavoredCSD0 = 0;
156int m_nFlavoredDCSD0 = 0;
157int m_nSL = 0;
158int m_nDalitz = 0;
159double m_dalitzNumer1 = 0;
160double m_dalitzNumer2 = 0;
161double m_dalitzDenom = 0;
162double dalitzNumer1_fil = 0;
163double dalitzNumer2_fil = 0;
164double dalitzDenom_fil = 0;
165HepSymMatrix m_weights(10,0);
166double m_rwsCF=0.;
167double m_rwsCS=0.;
168double m_deltaCF=0.;
169double m_deltaCS=0.;
170HepMatrix m_modeCounter(10,10,0);
171HepMatrix m_keptModeCounter(10,10,0);
172
173int m_rho_flag;
174
175
176/////////////////////////////////////////////////////////////////
177
178QCMCFilter::QCMCFilter(const std::string& name, ISvcLocator* pSvcLocator) :
179 Algorithm(name, pSvcLocator){
180 //Declare the properties
181 declareProperty("x", m_x=0.);
182 declareProperty("y", m_y= 0.0); //For values outside of (-.11 , .06) measured y value may vary from input y
183 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)
184 declareProperty("zForCFModes", m_zCF= 2.0);
185 declareProperty("wSignForCFModes", m_wCFSign= true);
186 declareProperty("rForCSModes", m_rCS= 1.0);
187 declareProperty("zForCSModes", m_zCS= -2.0);
188 declareProperty("wSignForCSModes", m_wCSSign= true);
189 declareProperty("DplusDminusWeight", m_dplusDminusWeight= 0.5);
190 declareProperty("dalitzModel", m_dalitzModel= 0); // 0 CLEO-c, 1 Babar, 2 Belle
191 declareProperty("UseNewWeights", m_useNewWeights= true);
192 declareProperty("InvertSelection", m_invertSelection= false);
193}
194
195//***************************************************************
196
197StatusCode QCMCFilter::initialize() {
198 MsgStream log(msgSvc(), name());
199
200 log << MSG::INFO << "in initialize()" << endmsg;
201
202 // Get the Particle Properties Service
203 IPartPropSvc* p_PartPropSvc;
204 static const bool CREATEIFNOTTHERE(true);
205 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
206 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc)
207 {
208 log << MSG::ERROR << " Could not initialize Particle Properties Service" << endmsg;
209 return PartPropStatus;
210 }
211 m_particleTable = p_PartPropSvc->PDT();
212
213
214 //Calculates parameters here based on Parameter Input by User
215 // (e.g. run expensive algorithms that are based on parameters
216 // specified by user at run-time)
217
218 m_y = (m_y - 0.008)/0.292;
219 //corrects the y input so that we get the desired measurable y from CP:Semi-leptonic doubletag method
220 // This equation is a result of solving equations for y dependence.
221 //If parent MC Branching ratio's this equation will need resolved. January 2013 -DA
222 double x = m_x ;
223 double y = m_y ;
224 double rCF = m_rCF ;
225 double zCF = m_zCF ;
226 double rCF2 = rCF * rCF ;
227 double zCF2 = zCF * zCF ;
228 double rCS = m_rCS ;
229 double zCS = m_zCS ;
230 double rCS2 = rCS * rCS ;
231 double zCS2 = zCS * zCS ;
232
233 double rmix = ( x * x + y * y ) / 2. ;
234 double wCF = ( m_wCFSign ? 1. : -1. ) * sqrt( 4. - zCF2 ) ;
235 double wCS = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 ) ;
236 double vCFCSPlus = ( zCF * zCS + wCF * wCS ) / 2. ;
237 double vCFCSMinus = ( zCF * zCS - wCF * wCS ) / 2. ;
238
239 m_deltaCF = acos( zCF / 2. ) * ( m_wCFSign ? 1. : -1. ) ;
240 m_deltaCS = acos( zCS / 2. ) * ( m_wCSSign ? 1. : -1. ) ;
241
242 double bCF = 1. + 0.5 * rCF * ( zCF * y + wCF * x ) + 0.5 * ( y * y - x * x );//second order y term added -DA
243 double bCS = 1. + 0.5 * rCS * ( zCS * y + wCS * x ) + 0.5 * ( y * y - x * x );
244 double bCPPlus = 1. - y ;
245 double bCPMinus = 1. + y ;
246
247 if( !m_useNewWeights )
248 {
249 bCF = 1. ;
250 m_rwsCF = rCF2 ;
251 bCS = 1. ;
252 m_rwsCS = rCS2 ;
253 bCPPlus = 1. ;
254 bCPMinus = 1. ;
255 }
256
257 m_rwsCF = ( rCF2 + 0.5 * rCF * ( zCF * y - wCF * x ) + rmix ) / bCF;//though division by b not mentioned in memo it is correct
258 m_rwsCS = ( rCS2 + 0.5 * rCS * ( zCS * y - wCS * x ) + rmix ) / bCS;
259
260
261 // Calculate reweighting factors, normalized to the largest factor.
262 // Note:A fake Dalitz weight is added in to just to make initial predictions
263 // These Dalitz weights are not used in actual code
264
265 // FlavoredCF
266
267 // First calculate rate factors.
268 m_weights[ kFlavoredCF ][ kFlavoredCFBar ] =
269 1. + rCF2 * ( 2. - zCF2 ) + rCF2 * rCF2 ;
270 m_weights[ kFlavoredCF ][ kFlavoredCF ] =
271 rmix * m_weights[ kFlavoredCF ][ kFlavoredCFBar ] ;
272 m_weights[ kFlavoredCF ][ kFlavoredCSBar ] =
273 1. + rCF2 * rCS2 - rCF * rCS * vCFCSMinus ;
274 m_weights[ kFlavoredCF ][ kFlavoredCS ] =
275 rCF2 + rCS2 - rCF * rCS * vCFCSPlus +
276 rmix * m_weights[ kFlavoredCF ][ kFlavoredCSBar ] ;
277
278 // Nov 2007: correct for r2 -> RWS and BR != A^2
279 m_weights[ kFlavoredCF ][ kFlavoredCF ] /= m_rwsCF * bCF * bCF ;
280 m_weights[ kFlavoredCF ][ kFlavoredCFBar ] /=
281 ( 1. + m_rwsCF * m_rwsCF ) * bCF * bCF ;
282 m_weights[ kFlavoredCF ][ kFlavoredCS ] /=
283 ( m_rwsCF + m_rwsCS ) * bCF * bCS ;
284 m_weights[ kFlavoredCF ][ kFlavoredCSBar ] /=
285 ( 1. + m_rwsCF * m_rwsCS ) * bCF * bCS ;
286
287 m_weights[ kFlavoredCF ][ kSLPlus ] = 1. / bCF ;
288 m_weights[ kFlavoredCF ][ kSLMinus ] = 1. / bCF ;
289
290 m_weights[ kFlavoredCF ][ kCPPlus ] =
291 ( 1. + rCF2 + rCF * zCF ) / ( 1. + m_rwsCF ) / bCF / bCPPlus ;
292 m_weights[ kFlavoredCF ][ kCPMinus ] =
293 ( 1. + rCF2 - rCF * zCF ) / ( 1. + m_rwsCF ) / bCF / bCPMinus ;
294 m_weights[ kFlavoredCF ][ kDalitz ] = 1. / bCF;
295
296 // FlavoredCFBar
297 m_weights[ kFlavoredCFBar ][ kFlavoredCFBar ] = m_weights[ kFlavoredCF ][ kFlavoredCF ] ;
298 m_weights[ kFlavoredCFBar ][ kFlavoredCS ] = m_weights[ kFlavoredCF ][ kFlavoredCSBar ] ;
299 m_weights[ kFlavoredCFBar ][ kFlavoredCSBar ] = m_weights[ kFlavoredCF ][ kFlavoredCS ] ;
300 m_weights[ kFlavoredCFBar ][ kSLPlus ] = 1. / bCF ;
301 m_weights[ kFlavoredCFBar ][ kSLMinus ] = 1. / bCF ;
302 m_weights[ kFlavoredCFBar ][ kCPPlus ] = m_weights[ kFlavoredCF ][ kCPPlus ] ;
303 m_weights[ kFlavoredCFBar ][ kCPMinus ] = m_weights[ kFlavoredCF ][ kCPMinus ] ;
304 m_weights[ kFlavoredCFBar ][ kDalitz ] = 1. / bCF;
305
306 // FlavoredCS
307
308 // First calculate rate factors.
309 m_weights[ kFlavoredCS ][ kFlavoredCSBar ] = 1. + rCS2 * ( 2. - zCS2 ) + rCS2 * rCS2 ;
310 m_weights[ kFlavoredCS ][ kFlavoredCS ] = rmix * m_weights[ kFlavoredCS ][ kFlavoredCSBar ] ;
311
312 // Nov 2007: correct for r2 -> RWS and BR != A^2
313 m_weights[ kFlavoredCS ][ kFlavoredCS ] /= m_rwsCS * bCS * bCS ;
314 m_weights[ kFlavoredCS ][ kFlavoredCSBar ] /= ( 1. + m_rwsCS * m_rwsCS ) * bCS * bCS ;
315
316 m_weights[ kFlavoredCS ][ kSLPlus ] = 1. / bCS ;
317 m_weights[ kFlavoredCS ][ kSLMinus ] = 1. / bCS ;
318
319 m_weights[ kFlavoredCS ][ kCPPlus ] = ( 1. + rCS2 + rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPPlus ;
320 m_weights[ kFlavoredCS ][ kCPMinus ] = ( 1. + rCS2 - rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPMinus ;
321
322 m_weights[ kFlavoredCS ][ kDalitz ] = 1. / bCS;
323
324 // FlavoredCSBar
325
326 m_weights[ kFlavoredCSBar ][ kFlavoredCSBar ] = m_weights[ kFlavoredCS ][ kFlavoredCS ] ;
327 m_weights[ kFlavoredCSBar ][ kSLPlus ] = 1. / bCS ;
328 m_weights[ kFlavoredCSBar ][ kSLMinus ] = 1. / bCS ;
329 m_weights[ kFlavoredCSBar ][ kCPPlus ] = m_weights[ kFlavoredCS ][ kCPPlus ] ;
330 m_weights[ kFlavoredCSBar ][ kCPMinus ] = m_weights[ kFlavoredCS ][ kCPMinus ] ;
331 m_weights[ kFlavoredCSBar ][ kDalitz ] = 1. / bCS;
332
333 // SLPlus
334
335 // SLPlus/SLPlus should have rate factor = Rmix, but none of these are
336 // generated in uncorrelated MC.
337 m_weights[ kSLPlus ][ kSLPlus ] = 0. ;
338 m_weights[ kSLPlus ][ kSLMinus ] = 1. ;
339 m_weights[ kSLPlus ][ kCPPlus ] = 1. / bCPPlus ;
340 m_weights[ kSLPlus ][ kCPMinus ] = 1. / bCPMinus ;
341 m_weights[ kSLPlus ][ kDalitz ] = 1. ;
342
343 // SLMinus
344
345 m_weights[ kSLMinus ][ kSLMinus ] = 0. ;
346 m_weights[ kSLMinus ][ kCPPlus ] = 1. / bCPPlus ;
347 m_weights[ kSLMinus ][ kCPMinus ] = 1. / bCPMinus ;
348 m_weights[ kSLMinus ][ kDalitz ] = 1. ;
349
350 // CPPlus
351
352 m_weights[ kCPPlus ][ kCPPlus ] = 0. ;
353 m_weights[ kCPPlus ][ kCPMinus ] = 2. / bCPPlus / bCPMinus ;
354 m_weights[ kCPPlus ][ kDalitz ] = 1. / bCPPlus;
355
356 // CPMinus
357
358 m_weights[ kCPMinus ][ kCPMinus ] = 0. ;
359 m_weights[ kCPMinus ][ kDalitz ] = 1. / bCPMinus;
360
361 //Dalitz estimate
362 m_weights[ kDalitz ][ kDalitz ] = 1;
363
364 log << MSG::FATAL << "Weight matrix" << m_weights << endmsg ;
365
366 //Boss 6.6.2 MC Branching fractions
367 HepVector fractionsD0( kNDecayTypes+1, 0 ) ;
368 fractionsD0[ kFlavoredCF ] = 0.531 ; //.510(CLEO-c BR for comparison)
369 fractionsD0[ kFlavoredCFBar ] = 0.0082 ; //.0103
370 fractionsD0[ kFlavoredCS ] = 0.031 ; //.0259
371 fractionsD0[ kFlavoredCSBar ] = 0.031 ; //.0259
372 fractionsD0[ kSLPlus ] = 0.125 ; //.129
373 fractionsD0[ kSLMinus ] = 0. ;
374 fractionsD0[ kCPPlus ] = 0.113 ; //.123
375 fractionsD0[ kCPMinus ] = 0.102 ; //.110
376 fractionsD0[ kDalitz ] = 0.0588 ; //.0578
377 HepVector scales = m_weights * fractionsD0 ;
378 log << MSG::INFO << "Integrated scale factors " <<scales << endmsg ;
379
380
381 HepVector fractionsD0bar( kNDecayTypes+1, 0 ) ;
382 fractionsD0bar[ kFlavoredCF ] = 0.0082 ;
383 fractionsD0bar[ kFlavoredCFBar ] = 0.531 ;
384 fractionsD0bar[ kFlavoredCS ] = 0.031 ;
385 fractionsD0bar[ kFlavoredCSBar ] = 0.031 ;
386 fractionsD0bar[ kSLPlus ] = 0. ;
387 fractionsD0bar[ kSLMinus ] = 0.125 ;
388 fractionsD0bar[ kCPPlus ] = 0.113 ;
389 fractionsD0bar[ kCPMinus ] = 0.102 ;
390 fractionsD0bar[ kDalitz ] = 0.0588 ;
391 HepVector weight_norm = fractionsD0bar.T() * scales;
392 log << MSG::INFO << "Overall scale factor " << fractionsD0bar.T() * scales << endmsg ;
393
394
395
396 //Original MC
397 //Branching fraction predictions after weight
398 double brCF = ( fractionsD0bar[ kFlavoredCFBar ] * scales[ kFlavoredCFBar ] ) / weight_norm[0] ;
399 double brCS = ( fractionsD0bar[ kFlavoredCS ] * scales[ kFlavoredCS ] + fractionsD0bar[ kFlavoredCSBar ] * scales[ kFlavoredCSBar ] ) / weight_norm[0] ;
400 double brDCS = ( fractionsD0bar[ kFlavoredCF ] * scales[ kFlavoredCF ] ) / weight_norm[0] ;
401 double brCPPlus = ( fractionsD0bar[ kCPPlus ] * scales[ kCPPlus ] ) / weight_norm[0] ;
402 double brCPMinus = ( fractionsD0bar[ kCPMinus ] * scales[ kCPMinus ] ) / weight_norm[0] ;
403
404 //y=0.0 dalitz example used for prediction
405 double dalitz_y_num = -0.000177536;
406 double dalitz_y_den = -0.0150846;
407
408 double numFactCF =
409 -rCF * zCF * ( 1. - 0.5 * rCF * wCF * m_x ) ;
410 double numFactCS =
411 -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x ) ;
412 double denFactCF = 0.5 * rCF2 * zCF2 ;
413 double denFactCS = 0.5 * rCS2 * zCS2 ;
414
415 double num_pre =
416 brCPPlus - brCPMinus + brCF * numFactCF + 0.5 * brCS * numFactCS + dalitz_y_num;
417 double den_pre =
418 1. - brCPPlus - brCPMinus - brCF * denFactCF - 0.5 * brCS * denFactCS + dalitz_y_den ;
419 double y_pre = num_pre / den_pre ;
420 double num =
421 fractionsD0bar[ kCPPlus ] - fractionsD0bar[ kCPMinus ] + fractionsD0bar[ kFlavoredCFBar ] * numFactCF + fractionsD0bar[ kFlavoredCS ] * numFactCS + dalitz_y_num;
422 double den =
423 1. - fractionsD0bar[ kCPPlus ] - fractionsD0bar[ kCPMinus ] - fractionsD0bar[ kFlavoredCFBar ] * denFactCF - fractionsD0bar[ kFlavoredCS ] * denFactCS + dalitz_y_den;
424 double y_prematrix = num / den ;
425 double y_test1 = 0.25 * ( ( ( scales[ kCPMinus ] * m_weights[ kCPPlus ][ kSLPlus ] ) / ( scales[ kCPPlus ] * m_weights[ kCPMinus ][ kSLPlus ] ) ) -
426 ( ( scales[ kCPPlus ] * m_weights[ kCPMinus ][ kSLPlus ] ) / ( scales[ kCPMinus ] * m_weights[ kCPPlus ][ kSLPlus ] ) ) ) ;
427
428 log << MSG::INFO << "An Input Y of " << m_y << endmsg ;
429 log << MSG::INFO << "Parent MC has an intrinsic value of y as " << y_prematrix << endmsg ;
430 log << MSG::INFO << "The BR method predics a y of " << y_pre << endmsg ;
431 log << MSG::INFO << "The CP-SL double tag method predicts y of " <<y_test1 << endmsg ;
432
433 // Renormalize
434 m_largestWeight = 0. ;
435 for( int i = 0 ; i < kNDecayTypes ; ++i )
436 {
437 for( int j = 0 ; j < kNDecayTypes ; ++j )
438 {
439 if( m_weights[ i ][ j ] < 0 ) m_weights[ i ][ j ] = 0; //gets rid of negative weights in early calculations
440 if( m_weights[ i ][ j ] > m_largestWeight )
441 {
442 m_largestWeight = m_weights[ i ][ j ] ;
443 }
444 }
445 }
446 m_weights /= m_largestWeight ;
447
448 log << MSG::INFO <<"Renormalized weight matrix " << m_weights << endmsg ;
449
450
451 //Set up weight parameters
452 Dalitz t(8);
453 double D0Sum[ 10 ], D0barSum[ 10 ] ; // last index is for sum over bins
454 TComplex D0D0barSum[ 10 ] ;
455 int points[ 10 ] ;
456
457 for( int i = 0 ; i < 10 ; ++i )
458 {
459 D0Sum[ i ] = 0. ;
460 D0barSum[ i ] = 0. ;
461 D0D0barSum[ i ] = TComplex( 0., 0. ) ;
462 points[ i ] = 0 ;
463 }
464
465 double m2min = 0. ;
466 double m2max = 3. ;
467 int nsteps = 1000 ;
468 double stepsize = ( m2max - m2min ) / ( double ) nsteps ;
469
470 for( int i = 0 ; i < nsteps ; ++i )
471 {
472 double m2i = m2min + ( ( double ) i + 0.5 ) * stepsize ;
473
474 for( int j = 0 ; j < nsteps ; ++j )
475 {
476 double m2j = m2min + ( ( double ) j + 0.5 ) * stepsize ;
477
478 if( t.Point_on_DP( m2i, m2j ) )
479 {
480 double m2k = 1.86450*1.86450 + 0.497648*0.497648 +
481 0.139570*0.139570 + 0.139570*0.139570 - m2i - m2j ;
482
483 TComplex D0, D0bar;
484 if (m_dalitzModel == 0) { //Cleo model
485 D0 = t.CLEO_Amplitude( m2i, m2j, m2k ) ;
486 D0bar = t.CLEO_Amplitude( m2j, m2i, m2k ) ;}
487 if (m_dalitzModel == 1) { //Babar model
488 D0 = t.Babar_Amplitude( m2i, m2j, m2k ) ;
489 D0bar = t.Babar_Amplitude( m2j, m2i, m2k ) ;}
490 if (m_dalitzModel == 2) { //Belle model
491 D0 = t.Amplitude( m2i, m2j, m2k ) ;
492 D0bar = t.Amplitude( m2j, m2i, m2k ) ;}
493
494 if ( j <= i ) {// lower half only
495 int bin = t.getBin( m2i, m2j, m2k ) ;
496 if( bin == -1 ) bin = 8 ;
497
498 ++points[ bin ] ;
499 D0Sum[ bin ] += norm( D0 ) ;
500 D0barSum[ bin ] += norm( D0bar ) ;
501 D0D0barSum[ bin ] += D0 * conj( D0bar ) ;
502
503 ++points[ 9 ] ;
504 D0Sum[ 9 ] += norm( D0 ) ;
505 D0barSum[ 9 ] += norm( D0bar ) ;
506 D0D0barSum[ 9 ] += D0 * conj( D0bar ) ;
507 }
508
509 }
510 }
511 }
512
513 for( int i = 0 ; i < 10 ; ++i )
514 {
515 double r2 = D0barSum[ i ] / D0Sum[ i ] ;
516 double c = real( D0D0barSum[ i ] ) / sqrt( D0barSum[ i ] ) / sqrt( D0Sum[ i ] ) ;
517 double s = imag( D0D0barSum[ i ] ) / sqrt( D0barSum[ i ] ) / sqrt( D0Sum[ i ] ) ;
518
519 cout << "BIN " << i << ": "
520 << points[ i ] << " pts "
521 << r2 << " " << c << " " << s << " " << c*c+s*s
522 << endmsg ;
523 }
524
525 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
526 return StatusCode::SUCCESS;
527}
528
529//********************************************************************
530StatusCode QCMCFilter::execute() {
531
532 //interface to event data service
533 ISvcLocator* svcLocator = Gaudi::svcLocator();
534 StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
535 if (sc.isFailure())
536 std::cout<<"Could not access EventDataSvc!"<<std::endl;
537
538 setFilterPassed(false);
539
540 MsgStream log(msgSvc(), name());
541 log<< MSG::INFO << "in execute()" << endmsg;
542
543 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
544
545 if(!eventHeader)
546 {
547 cout<<" eventHeader "<<endl;
548 return StatusCode::FAILURE;
549 }
550 long runNo = eventHeader->runNumber();
551 long event = eventHeader->eventNumber();
552 log << MSG::DEBUG << "run, evtnum = "
553 << runNo << " , "
554 << event << endmsg;
555
556 // Get McParticle collection
557 SmartDataPtr<Event::McParticleCol> mcParticles(eventSvc(), EventModel::MC::McParticleCol);
558
559 //Check if event is Psi(3770)
560 int psipp_daughter = 0;
561 int d0_count = 0;
562 int chd_count = 0;
563 for( Event::McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end(); it++)
564 {
565 int pdg_code = (*it)->particleProperty();
566 if ((*it)->primaryParticle()) continue;
567 Event::McParticle it2 = (*it)->mother();
568 int mother_pdg = 0;
569 //finds if the particle is daughter of Psi(3770)
570 mother_pdg = it2.particleProperty();
571 if (mother_pdg == kPsi3770ID)
572 {
573 psipp_daughter++; //Found daughter
574 if (abs(pdg_code) == kD0ID) d0_count++;
575 if (abs(pdg_code) == kDpID) chd_count++;
576 if (pdg_code == kD0BarID) m_nD0bar++;
577 }
578 }
579
580 if( psipp_daughter == 0 ) //didn't find Psi(3770)
581 {
582 ++m_nUnknownEvents ;
583
584 if( m_invertSelection )
585 {
586 log<< MSG::DEBUG << "Discarding event -- did not find psi(3770). " << endmsg ;
587 return StatusCode::SUCCESS; //discard event
588 }
589 else
590 {
591 // -------- Write to root file
592 setFilterPassed(true);
593
594 return StatusCode::SUCCESS; //kept event
595 }
596 }
597
598 log<< MSG::DEBUG << "Found psi(3770) -->" << endmsg ;
599
600 // Check that top particle has only two daughters + possible FSR photon
601 if ( psipp_daughter > 3 )
602 {
603 ++m_nUnknownEvents ;
604
605 if( m_invertSelection )
606 {
607 log<< MSG::DEBUG << "Discarding event -- psi(3770) has >3 daughters." << endmsg ;
608 return StatusCode::SUCCESS; //discard event
609 }
610 else
611 {
612 // -------- Write to root file
613 setFilterPassed(true);
614
615
616 return StatusCode::SUCCESS; //kept event
617 }
618 }
619
620 //Check if D+D- state
621 if (chd_count == 2) { // D+D- event
622 ++m_nDpDmEvents ;
623
624 // D+D- must be rewieghted otherwise too many D+D- events relative to D0D0bar.
625 double random = RandFlat::shoot() ;
626
627 if( ( random < m_dplusDminusWeight && !m_invertSelection ) ||
628 ( random > m_dplusDminusWeight && m_invertSelection ) )
629 {
630 // -------- Write to root file
631 setFilterPassed(true);
632
633
634 return StatusCode::SUCCESS; //event kept
635 }
636 else
637 {
638 log<< MSG::DEBUG << "Discarding event -- D+D- event failed reweighting." << endmsg ;
639 ++m_nDpDmDiscarded ;
640 return StatusCode::SUCCESS; //event discarded
641 }
642 }
643
644 //Check if D0D0Bar
645 if (d0_count != 2) {
646 if( !m_invertSelection )
647 {
648 log<< MSG::DEBUG << "Discarding event -- non-D0-D0bar event." << endmsg ;
649 return StatusCode::SUCCESS; //discard event
650 }
651 else
652 {
653 // -------- Write to root file
654 setFilterPassed(true);
655
656
657 return StatusCode::SUCCESS; //kept event
658 }
659 }
660
661 log<< MSG::INFO << "D0-D0bar event." << endmsg ;
662 m_nD0D0barEvents++ ;
663
664 m_rho_flag = 0;
665
666 //Get D0-D0bar modes and info
667 vector<double> d0_decay = findD0Decay(1);
668 vector<double> d0bar_decay = findD0Decay(-1);
669
670 int d0mode = d0_decay[0];
671 int d0barmode = d0bar_decay[0];
672 m_modeCounter[d0mode][d0barmode]++;
673
674
675
676 //if any event unidentified remove
677 if (d0_decay[0] == 9 || d0bar_decay[0] == 9 ) {
678 if( !m_invertSelection )
679 {
680 log<< MSG::DEBUG << "Discarding event -- unknown D0-D0bar event." << endmsg ;
681 return StatusCode::SUCCESS; //discard event
682 }
683 else
684 {
685 // -------- Write to root file
686 setFilterPassed(true);
687
688
689
690 return StatusCode::SUCCESS; //kept event
691 }
692 }
693
694 // Loop over D's
695 double r2D0 = d0_decay[1] ;
696 double deltaD0 = d0_decay[2] ;
697 double r2D0bar = d0bar_decay[1] ;
698 double deltaD0bar = d0bar_decay[2] ;
699
700
701 //Weight based on DDbar combination
702 double weight ;
703 if( d0_decay[0] == kDalitz || d0bar_decay[0] == kDalitz )
704 {
705 double r2prod = r2D0 * r2D0bar ;
706 double x = m_x ;
707 double y = m_y ;
708
709 double bD0 = 1. + sqrt( r2D0 ) * ( cos( deltaD0 ) * y + sin( deltaD0 ) * x ) ;
710 double rwsD0 = ( r2D0 + sqrt( r2D0 ) * ( cos( deltaD0 ) * y - sin( deltaD0 ) * x ) ) / bD0 ;
711 double bD0bar = 1. + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y + sin( deltaD0bar ) * x ) ;
712 double rwsD0bar = ( r2D0bar + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y - sin( deltaD0bar ) * x ) ) / bD0bar ;
713
714 weight = 1. + r2prod - 2. * sqrt( r2prod ) * cos( deltaD0 + deltaD0bar ) ;
715 weight /= ( 1. + rwsD0 * rwsD0bar ) * bD0 * bD0bar ;
716 weight /= m_largestWeight ;
717
718 //get dalitz graph information before cut
719 int k = -10;
720 double m2i= 0;
721 double m2j= 0;
722 double m2k= 0;
723 if( d0_decay[0] == kDalitz) {
724 k = d0bar_decay[0];
725 m2i= d0_decay[3];
726 m2j= d0_decay[4];
727 m2k= d0_decay[6];}
728 if( d0bar_decay[0] == kDalitz) {
729 k = d0_decay[0];
730 m2i= d0bar_decay[3];
731 m2j= d0bar_decay[4];
732 m2k= d0bar_decay[6];}
733
734 }
735 else
736 {
737 weight = m_weights[ d0_decay[0] ][ d0bar_decay[0] ] ;
738 }
739
740 double random = RandFlat::shoot() ;
741 log<< MSG::DEBUG << "type: " << d0_decay[0] << " " << d0bar_decay[0] << ", weight " <<
742 weight << " <? random " << random << endmsg ;
743
744 if( ( random < weight && !m_invertSelection ) || ( random > weight && m_invertSelection ) )
745 {
746 // -------- Write to root file
747 setFilterPassed(true);
748
749 //get dalitz graph information after filter
750 int k = -10;
751 double m2i= 0;
752 double m2j= 0;
753 double m2k= 0;
754 double cosd = cos( deltaD0 ) ;
755 double sind = sin( deltaD0 ) ;
756 if( d0_decay[0] == kDalitz) {
757 k = d0bar_decay[0];
758 m2i= d0_decay[3];
759 m2j= d0_decay[4];
760 m2k= d0_decay[6];
761 if ( m2j - m2i > 0. ) { //avoids doublecounting the branching fraction
762 dalitzNumer1_fil += -2. * sqrt( r2D0 ) * cosd ;
763 dalitzNumer2_fil += 2. * r2D0 * cosd * sind * m_x ;
764 dalitzDenom_fil += -2. * r2D0 * cosd * cosd ; }
765 }
766 if( d0bar_decay[0] == kDalitz) {
767 k = d0_decay[0];
768 m2i= d0bar_decay[3];
769 m2j= d0bar_decay[4];
770 m2k= d0bar_decay[6];
771 cosd = cos( deltaD0bar ) ;
772 sind = sin( deltaD0bar ) ;
773 if( m2j - m2i < 0. ) { //avoids doublecounting the branching fraction
774 dalitzNumer1_fil += -2. * sqrt( r2D0bar ) * cosd ;
775 dalitzNumer2_fil += 2. * r2D0bar * cosd * sind * m_x ;
776 dalitzDenom_fil += -2. * r2D0bar * cosd * cosd ; }
777 }
778
779 m_keptModeCounter[d0mode][d0barmode]++;
780 return StatusCode::SUCCESS; //event kept
781 }
782 else
783 {
784 log << MSG::DEBUG << "Discarding event -- failed reweighting." << endmsg ;
785
786 ++m_nD0D0barDiscarded ;
787 return StatusCode::SUCCESS; //event discarded
788 }
789}
790
791
792//*********************************************************************
793StatusCode QCMCFilter::finalize() {
794 MsgStream log(msgSvc(), name());
795 log << MSG::INFO << "in finalize()" << endmsg;
796 log << MSG::INFO << "Number of unknown events: " << m_nUnknownEvents << endmsg ;
797 log << MSG::INFO << "Number of D+D- events seen: " << m_nDpDmEvents << endmsg ;
798 log << MSG::INFO << "Number of D+D- events discarded: " << m_nDpDmDiscarded << endmsg ;
799 if( m_nDpDmEvents > 0 )
800 {
801 log << MSG::INFO <<"Fraction discarded: " << double( m_nDpDmDiscarded ) / double( m_nDpDmEvents ) << endmsg ;
802 }
803
804 log << MSG::INFO << "Number of D0D0bar events seen: " << m_nD0D0barEvents << " " << m_nD0bar << endmsg ;
805 log << MSG::INFO << "Number of D0D0bar events discarded: " << m_nD0D0barDiscarded << endmsg ;
806 if( m_nD0D0barEvents > 0 && m_nCPPlus > 0 && m_nCPMinus > 0)
807 {
808 log << MSG::INFO << "Fraction discarded: " << double( m_nD0D0barDiscarded ) / double( m_nD0D0barEvents ) << endl <<
809 "Fraction unknown decays: " << 0.5 * double( m_nUnknownDecays ) / double( m_nD0D0barEvents ) <<endmsg ;
810
811 double nD0 = 2. * m_nD0D0barEvents ;
812 double brSL = m_nSL / nD0 ;
813 double dBrSL = sqrt( brSL * ( 1. - brSL ) / nD0 ) ;
814 double brCF = m_nFlavoredCFD0 / nD0 ;
815 double dBrCF = sqrt( brCF * ( 1. - brCF ) / nD0 ) ;
816 double brCS = m_nFlavoredCSD0 / nD0 ;
817 double dBrCS = sqrt( brCS * ( 1. - brCS ) / nD0 ) ;
818 double brDCS = m_nFlavoredDCSD0 / nD0 ;
819 double dBrDCS = sqrt( brDCS * ( 1. - brDCS ) / nD0 ) ;
820 double brCPPlus = m_nCPPlus / nD0 ;
821 double dBrCPPlus = sqrt( brCPPlus * ( 1. - brCPPlus ) / nD0 ) ;
822 double brCPMinus = m_nCPMinus / nD0 ;
823 double dBrCPMinus = sqrt( brCPMinus * ( 1. - brCPMinus ) / nD0 ) ;
824 double brDalitz = m_nDalitz / nD0 ;
825 double dBrDalitz = sqrt( brDalitz * ( 1. - brDalitz ) / nD0 ) ;
826 double fdBrDalitz = dBrDalitz / brDalitz ;
827 double dalitzNumer1Norm = m_dalitzNumer1 / nD0 ;
828 double dDalitzNumer1Norm = dalitzNumer1Norm * fdBrDalitz ;
829 double dalitzNumer2Norm = m_dalitzNumer2 / nD0 ;
830 double dDalitzNumer2Norm = dalitzNumer2Norm * fdBrDalitz ;
831 double dalitzDenomNorm = m_dalitzDenom / nD0 ;
832 double dDalitzDenomNorm = dalitzDenomNorm * fdBrDalitz ;
833
834 //after the filter was applied
835 double fil_SL = 0, fil_FlavoredCFD0 = 0, fil_FlavoredCSD0 = 0, fil_FlavoredDCSD0 =0, fil_CPPlus = 0, fil_CPMinus = 0, fil_Dalitz = 0;
836 for( int i = 0 ; i < 9 ; ++i )
837 {
838 fil_SL += m_keptModeCounter[i][4] + m_keptModeCounter[4][i] + m_keptModeCounter[i][5] + m_keptModeCounter[5][i] ;
839 fil_FlavoredCFD0 += m_keptModeCounter[i][1] + m_keptModeCounter[0][i];
840 fil_FlavoredCSD0 += m_keptModeCounter[i][2] + m_keptModeCounter[2][i] + m_keptModeCounter[i][3] + m_keptModeCounter[3][i];
841 fil_FlavoredDCSD0 += m_keptModeCounter[i][0] + m_keptModeCounter[1][i];
842 fil_CPPlus += m_keptModeCounter[i][6] + m_keptModeCounter[6][i];
843 fil_CPMinus += m_keptModeCounter[i][7] + m_keptModeCounter[7][i];
844 fil_Dalitz += m_keptModeCounter[i][8] + m_keptModeCounter[8][i];
845 }
846 double nD0_fil = 2. * ( m_nD0D0barEvents - m_nD0D0barDiscarded ) ;
847 double brSL_fil = fil_SL / nD0_fil ;
848 double dBrSL_fil = sqrt( brSL_fil * ( 1. - brSL_fil ) / nD0_fil ) ;
849 double brCF_fil = fil_FlavoredCFD0 / nD0_fil ;
850 double dBrCF_fil = sqrt( brCF_fil * ( 1. - brCF_fil ) / nD0_fil ) ;
851 double brCS_fil = fil_FlavoredCSD0 / nD0_fil ;
852 double dBrCS_fil = sqrt( brCS_fil * ( 1. - brCS_fil ) / nD0_fil ) ;
853 double brDCS_fil = fil_FlavoredDCSD0 / nD0_fil ;
854 double dBrDCS_fil = sqrt( brDCS_fil * ( 1. - brDCS_fil ) / nD0_fil ) ;
855 double brCPPlus_fil = fil_CPPlus / nD0_fil ;
856 double dBrCPPlus_fil = sqrt( brCPPlus_fil * ( 1. - brCPPlus_fil ) / nD0_fil ) ;
857 double brCPMinus_fil = fil_CPMinus / nD0_fil ;
858 double dBrCPMinus_fil = sqrt( brCPMinus_fil * ( 1. - brCPMinus_fil ) / nD0_fil ) ;
859 double brDalitz_fil = fil_Dalitz / nD0_fil ;
860 double dBrDalitz_fil = sqrt( brDalitz_fil * ( 1. - brDalitz_fil ) / nD0_fil ) ;
861 double fdBrDalitz_fil = dBrDalitz_fil / brDalitz_fil ;
862 double dalitzNumer1Norm_fil = dalitzNumer1_fil / nD0_fil ;
863 double dDalitzNumer1Norm_fil = dalitzNumer1Norm_fil * fdBrDalitz_fil ;
864 double dalitzNumer2Norm_fil = dalitzNumer2_fil / nD0_fil ;
865 double dDalitzNumer2Norm_fil = dalitzNumer2Norm_fil * fdBrDalitz_fil ;
866 double dalitzDenomNorm_fil = dalitzDenom_fil / nD0_fil ;
867 double dDalitzDenomNorm_fil = dalitzDenomNorm_fil * fdBrDalitz_fil ;
868
869
870 log << MSG::INFO <<
871 "Original number of SL decays: " << m_nSL <<" ==> Br(SL) = " << brSL << " +- " << dBrSL << endl <<
872 "Filtered number of SL decays: " << fil_SL <<" ==> Br(SL) = " << brSL_fil << " +- " << dBrSL_fil << endl <<
873 "Original number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: " <<m_nFlavoredCFD0 << "/" << m_nFlavoredCSD0 << "/" << m_nFlavoredDCSD0 << endl <<
874 " ==> Br(CF) = " << brCF << " +- " << dBrCF << endl <<
875 " ==> Br(CS) = " << brCS << " +- " << dBrCS << endl <<
876 " ==> Br(DCS) = " << brDCS << " +- " << dBrDCS << endl <<
877 "Filtered number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: " <<fil_FlavoredCFD0 << "/" << fil_FlavoredCSD0 << "/" << fil_FlavoredDCSD0 << endl <<
878 " ==> Br(CF) = " << brCF_fil << " +- " << dBrCF_fil << endl <<
879 " ==> Br(CS) = " << brCS_fil << " +- " << dBrCS_fil << endl <<
880 " ==> Br(DCS) = " << brDCS_fil << " +- " << dBrDCS_fil << endl <<
881
882 "Original number of CP+/-: " << m_nCPPlus << "/" << m_nCPMinus <<endl <<
883 " ==> Br(CP+) = " << brCPPlus << " +- " << dBrCPPlus <<endl <<
884 " ==> Br(CP-) = " << brCPMinus << " +- " << dBrCPMinus << endl <<
885 "Filtered number of CP+/-: " << fil_CPPlus << "/" << fil_CPMinus <<endl <<
886 " ==> Br(CP+) = " << brCPPlus_fil << " +- " << dBrCPPlus_fil <<endl <<
887 " ==> Br(CP-) = " << brCPMinus_fil << " +- " << dBrCPMinus_fil << endl <<
888
889 "Original number of Dalitz: " << m_nDalitz << endl <<
890 " ==> Br(Dalitz) = " << brDalitz << " +- " << dBrDalitz <<endl <<
891 " y contrib. numer1 " << dalitzNumer1Norm << " +- " << dDalitzNumer1Norm << endl <<
892 " numer2 " << dalitzNumer2Norm << " +- " << dDalitzNumer2Norm << endl <<
893 " denom " << dalitzDenomNorm << " +- " << dDalitzDenomNorm << endmsg <<
894 "Filtered number of Dalitz: " << fil_Dalitz << endl <<
895 " ==> Br(Dalitz) = " << brDalitz_fil << " +- " << dBrDalitz_fil <<endl <<
896 " y contrib. numer1 " << dalitzNumer1Norm_fil << " +- " << dDalitzNumer1Norm_fil << endl <<
897 " numer2 " << dalitzNumer2Norm_fil << " +- " << dDalitzNumer2Norm_fil << endl <<
898 " denom " << dalitzDenomNorm_fil << " +- " << dDalitzDenomNorm_fil << endmsg ;
899
900 cout<<""<<endl;
901
902 HepMatrix differencetoWeight(10,10,0);
903 for( int i = 0 ; i < 9 ; ++i ) {
904 for( int j = 0 ; j < 9 ; ++j ) {
905 if (m_modeCounter[i][j] != 0) {
906 differencetoWeight[i][j] = m_keptModeCounter[i][j] / m_modeCounter[i][j] - m_weights[i][j]; } } }
907
908 log << MSG::INFO << "Weight matrix" << m_weights << endmsg ;
909 log << MSG::INFO << "D0 Modes before filter" << m_modeCounter << endmsg ;
910 log << MSG::INFO << "D0 Modes after filter" << m_keptModeCounter << endmsg ;
911 log << MSG::INFO << "Compare percentage difference to weights " << differencetoWeight << endmsg ;
912
913 //Calculating and comparing Y before/after the filter
914 // To calculate y, we want half the CS BF
915 brCS /= 2. ;
916 dBrCS /= 2. ;
917 brCS_fil /= 2 ;
918 dBrCS_fil /= 2 ;
919
920 double y, dy, y_fil, dy_fil ;
921 double y_cpsl, dy_cpsl, y_fil_cpsl, dy_fil_cpsl ;
922
923 //Original MC
924 double rCF = m_rCF ;
925 double zCF = m_zCF ;
926 double rCF2 = rCF * rCF ;
927 double zCF2 = zCF * zCF ;
928 double rCS = m_rCS ;
929 double zCS = m_zCS ;
930 double rCS2 = rCS * rCS ;
931 double zCS2 = zCS * zCS ;
932
933 double wCF = ( m_wCFSign ? 1. : -1. ) * sqrt( 4. - zCF2 ) ;
934 double wCS = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 ) ;
935
936 double numFactCF =
937 -rCF * zCF * ( 1. - 0.5 * rCF * wCF * m_x ) ;
938 double numFactCS =
939 -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x ) ;
940 double denFactCF = 0.5 * rCF2 * zCF2 ;
941 double denFactCS = 0.5 * rCS2 * zCS2 ;
942 double dalitzNumerNorm = dalitzNumer1Norm + dalitzNumer2Norm ;
943
944 double num =
945 brCPPlus - brCPMinus + brCF * numFactCF + brCS * numFactCS +
946 dalitzNumerNorm ;
947 double den =
948 1. - brCPPlus - brCPMinus - brCF * denFactCF - brCS * denFactCS +
949 dalitzDenomNorm ;
950 y = num / den ;
951 dy = sqrt(
952 ( num + den ) * ( num + den ) * dBrCPPlus * dBrCPPlus +
953 ( num - den ) * ( num - den ) * dBrCPMinus * dBrCPMinus +
954 ( numFactCF * den + denFactCF * num ) * dBrCF *
955 ( numFactCF * den + denFactCF * num ) * dBrCF +
956 ( numFactCS * den + denFactCS * num ) * dBrCS *
957 ( numFactCS * den + denFactCS * num ) * dBrCS +
958 ( dalitzNumerNorm * den + dalitzDenomNorm * num ) * fdBrDalitz *
959 ( dalitzNumerNorm * den + dalitzDenomNorm * num ) * fdBrDalitz
960 ) / ( den * den ) ;
961
962 double n_cpplussl = m_modeCounter[6][4] + m_modeCounter[6][5] + m_modeCounter[4][6] + m_modeCounter[5][6] ;
963 double n_cpminussl = m_modeCounter[7][4] + m_modeCounter[7][5] + m_modeCounter[4][7] + m_modeCounter[5][7] ;
964 double a_cpsl = ( n_cpplussl * m_nCPMinus ) / ( n_cpminussl * m_nCPPlus ) ;
965 double b_cpsl = ( n_cpminussl * m_nCPPlus ) / ( n_cpplussl * m_nCPMinus ) ;
966 y_cpsl = 0.25 * ( a_cpsl - b_cpsl ) ;
967 dy_cpsl = 0.25 * ( a_cpsl + b_cpsl ) * sqrt( ( 1/n_cpplussl ) + ( 1/n_cpminussl ) + ( 1/m_nCPPlus ) + ( 1/m_nCPMinus ) ) ;
968
969
970 //Filtered MC
971 double dalitzNumerNorm_fil = dalitzNumer1Norm_fil + dalitzNumer2Norm_fil ;
972 double num_fil =
973 brCPPlus_fil - brCPMinus_fil + brCF_fil * numFactCF + brCS_fil * numFactCS +
974 dalitzNumerNorm_fil ;
975 double den_fil =
976 1. - brCPPlus_fil - brCPMinus_fil - brCF_fil * denFactCF - brCS_fil * denFactCS +
977 dalitzDenomNorm_fil ;
978 y_fil = num_fil / den_fil ;
979 dy_fil = sqrt(
980 ( num_fil + den_fil ) * ( num_fil + den_fil ) * dBrCPPlus_fil * dBrCPPlus_fil +
981 ( num_fil - den_fil ) * ( num_fil - den_fil ) * dBrCPMinus_fil * dBrCPMinus_fil +
982 ( numFactCF * den_fil + denFactCF * num_fil ) * dBrCF_fil *
983 ( numFactCF * den_fil + denFactCF * num_fil ) * dBrCF_fil +
984 ( numFactCS * den_fil + denFactCS * num_fil ) * dBrCS_fil *
985 ( numFactCS * den_fil + denFactCS * num_fil ) * dBrCS_fil +
986 ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil * num_fil ) * fdBrDalitz_fil *
987 ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil * num_fil ) * fdBrDalitz_fil
988 ) / ( den_fil * den_fil ) ;
989
990 double fil_cpplussl = m_keptModeCounter[6][4] + m_keptModeCounter[6][5] + m_keptModeCounter[4][6] + m_keptModeCounter[5][6] ;
991 double fil_cpminussl = m_keptModeCounter[7][4] + m_keptModeCounter[7][5] + m_keptModeCounter[4][7] + m_keptModeCounter[5][7] ;
992 a_cpsl = ( fil_cpplussl * fil_CPMinus ) / ( fil_cpminussl * fil_CPPlus ) ;
993 b_cpsl = ( fil_cpminussl * fil_CPPlus ) / ( fil_cpplussl * fil_CPMinus ) ;
994 y_fil_cpsl = 0.25 * ( a_cpsl - b_cpsl ) ;
995 dy_fil_cpsl = 0.25 * ( a_cpsl + b_cpsl ) * sqrt( ( 1/n_cpplussl ) + ( 1/n_cpminussl ) + ( 1/fil_CPPlus ) + ( 1/fil_CPMinus ) ) ;
996
997
998
999 log << MSG::INFO <<"BR Method : Parent MC ==> y = " << y << " +- " << dy << endmsg ;
1000 log << MSG::INFO <<"BR Method : Filtered MC ==> y = " << y_fil << " +- " << dy_fil << endmsg ;
1001 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 corrilated"<< endmsg ;
1002 log << MSG::INFO <<"CP-SL doubletag : Filtered MC ==> y = " << y_fil_cpsl << " +- " << dy_fil_cpsl << endmsg ;
1003
1004 }
1005
1006 return StatusCode::SUCCESS;
1007}
1008
1009
1010///////////////////////////////////////////////////////////////////////////////////////////
1011//////////////////////////////////////////////////////////////////////////////////////////
1012vector<double> QCMCFilter::findD0Decay(int charm)
1013{
1014 MsgStream log(msgSvc(), name());
1015 log << MSG::INFO <<" In findD0Decay " << endmsg ;
1016
1017 vector<double> d0_info(7);
1018 for(int i=0; i< 7; i++) d0_info[i]=0;
1019 double decay_mode = 9; //d0_info[0] = decay mode;
1020 double r2D0 = -1; //d0_info[1] = r2D0;
1021 double deltaD0 = -1; //d0_info[2] = deltaD0;
1022
1023 double num[26];
1024 for(int i=0; i< 26; i++) num[i]=0;
1025 // Number of each partile type
1026 int kPiPlus = 0;// num[0] = Pi+, a0+
1027 int kPiMinus = 1;// num[1] = Pi-, a0-
1028 int kPi0 = 2;// num[2] = Pi0
1029 int kEta = 3;// num[3] = Eta
1030 int kEtaPrime = 4;// num[4] = Eta Prime
1031 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)
1032 int kUFVPlus = 6;// num[6] = unflavored positive (axial) vectors: rho+, a1+, rho(2S)+
1033 int kUFVMinus = 7;// num[7] = unflavored negative (axial) vectors: rho-, a1-, rho(2S)-
1034 int kRho0 = 8;// num[8] = Rho0, Rho(2S)0
1035 int kOmega = 9;// num[9] = Omega
1036 int kPhi = 10;// num[10] = Phi
1037 int kKPlus = 11;// num[11] = K+
1038 int kKMinus = 12;// num[12] = K-
1039 int kK0S = 13;// num[13] = K short
1040 int kK0L = 14;// num[14] = K long
1041 int kFVPlus = 15;// num[15] = flavored positive (axial) vectors: K*+, K_1+, K_2*+
1042 int kFVMinus = 16;// num[16] = flavored negative (axial) vectors: K*-, K_1-, K_2*-
1043 int kKStar0 = 17;// num[17] = K*0
1044 int kKStar0Bar = 18;// num[18] = anti-K*0
1045 int kK10 = 19;// num[19] = K1_0
1046 int kK10Bar = 20;// num[20] = anti-K1_0
1047 int kLeptonPlus = 21;// num[21] = LeptonPlus: e+, mu+
1048 int kLeptonMinus = 22;// num[22] = LeptonMinus: e-, mu-
1049 int kNeutrino = 23;// num[23] = Neutrino: nu, nubar
1050 int kGamma = 24;// num[24] = gamma
1051 int kUnknown = 25;// num[25] = Not identified particle
1052 //int kGammaFSR GammaFSR is ignored as we are interested in the PreFSR quantum #
1053
1054 // Get McParticle collection
1055 SmartDataPtr<Event::McParticleCol> mcParticles(eventSvc(), EventModel::MC::McParticleCol);
1056 int d0_pdg = charm == 1 ? kD0ID : kD0BarID;
1057
1058 HepLorentzVector piPlus, piMinus, k0;
1059
1060 for( Event::McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end(); it++)
1061 {
1062 int pdg_code = (*it)->particleProperty();
1063 if ((*it)->primaryParticle()) continue;
1064 Event::McParticle it2 = (*it)->mother();
1065 int mother_pdg = 0;
1066
1067 //finds if the particle is related to the d0 we need
1068 mother_pdg = it2.particleProperty();
1069
1070 //special cases
1071 if (pdg_code == kKStar0ID || pdg_code == kKStar0BarID || pdg_code == kK0ID || pdg_code == kK0BarID
1072 || pdg_code == kK10ID || pdg_code == kK10BarID
1073 || pdg_code == kK0Star0ID || pdg_code == kK0Star0BarID || pdg_code == kK0StarPlusID || pdg_code == kK0BarID ) continue; //don't count special intermediate states
1074
1075 if (mother_pdg == kK0ID || mother_pdg == kK0BarID) {//Looks another step up if mother was k0 or k0bar
1076 it2 = it2.mother();
1077 mother_pdg = it2.particleProperty(); }
1078
1079 if (mother_pdg == kKStar0ID || mother_pdg == kKStar0BarID) {//Looks another step up if mother was k*0 or k*0bar
1080 //if code is found to be K*0 -> K+Pi- we save it as a KStar0 and K*0Bar -> K-Pi+ as KStar0Bar
1081 //if is it a neutral decay K*0(Bar)->K0(bar) Pi0 or Gamma. We save them as K0(bar) and Pi0 or Gamma.
1082 if (pdg_code == kPiPlusID || pdg_code == kPiMinusID) continue;
1083 if (mother_pdg == kKStar0ID && pdg_code == kKPlusID) pdg_code = kKStar0ID;
1084 if (mother_pdg == kKStar0BarID && pdg_code == kKMinusID) pdg_code = kKStar0BarID;
1085 it2 = it2.mother();
1086 mother_pdg = it2.particleProperty();
1087 }
1088
1089 if ( mother_pdg == kK0Star0ID || mother_pdg == kK0Star0BarID ) { // Looks another step up if mother was K_0
1090 it2 = it2.mother();
1091 mother_pdg = it2.particleProperty();
1092 }
1093
1094 if (mother_pdg == kK10ID || mother_pdg == kK10BarID) {//Looks another step up if mother was k_10 or k_10bar
1095 it2 = it2.mother();
1096 mother_pdg = it2.particleProperty();
1097 }
1098
1099
1100
1101 //Identifies what particles are the daughters
1102 if (mother_pdg == d0_pdg) //Found daughter
1103 {
1104 if (pdg_code == kPiPlusID || pdg_code == kA0PlusID ) num[0]++;
1105 else if (pdg_code == kPiMinusID || pdg_code == kA0MinusID ) num[1]++;
1106 else if (pdg_code == kPi0ID ) num[2]++;
1107 else if (pdg_code == kEtaID ) num[3]++;
1108 else if (pdg_code == kEtaPrimeID ) num[4]++;
1109 else if (pdg_code == kF0ID || pdg_code == kA00ID
1110 || pdg_code == kFPrime0ID|| pdg_code == kF0m1500ID
1111 || pdg_code == kF2ID) num[5]++;
1112 else if (pdg_code == kRhoPlusID || pdg_code == kRho2SPlusID
1113 || pdg_code == kA1PlusID ) num[6]++;
1114 else if (pdg_code == kRhoMinusID || pdg_code == kRho2SMinusID
1115 || pdg_code == kA1MinusID) num[7]++;
1116 else if (pdg_code == kRho0ID || pdg_code == kRho2S0ID) num[8]++;
1117 else if (pdg_code == kOmegaID ) num[9]++;
1118 else if (pdg_code == kPhiID ) num[10]++;
1119 else if (pdg_code == kKPlusID ) num[11]++;
1120 else if (pdg_code == kKMinusID ) num[12]++;
1121 else if (pdg_code == kK0SID ) num[13]++;
1122 else if (pdg_code == kK0LID ) num[14]++;
1123 else if (pdg_code == kKStarPlusID || pdg_code == kK1PlusID
1124 || pdg_code == kK2StarPlusID || pdg_code == kK1PrimePlusID
1125 || pdg_code == kK0StarPlusID) num[15]++;
1126 else if (pdg_code == kKStarMinusID || pdg_code == kK1MinusID
1127 || pdg_code == kK2StarMinusID || pdg_code == kK1PrimeMinusID
1128 || pdg_code == kK0StarMinusID ) num[16]++;
1129 else if (pdg_code == kKStar0ID ) num[17]++;
1130 else if (pdg_code == kKStar0BarID ) num[18]++;
1131 else if (pdg_code == kK10ID || pdg_code == kK1Prime0ID) num[19]++;
1132 else if (pdg_code == kK10BarID || pdg_code == kK1Prime0BarID) num[20]++;
1133 else if (pdg_code == kEPlusID || pdg_code == kMuPlusID ) num[21]++;
1134 else if (pdg_code == kEMinusID || pdg_code == kMuMinusID ) num[22]++;
1135 else if (pdg_code == kNuEID || pdg_code == kNuEBarID
1136 || pdg_code == kNuMuID || pdg_code == kNuMuBarID ) num[23]++;
1137 else if (pdg_code == kGammaID ) num[24]++;
1138 else if (pdg_code == kGammaFSRID ) continue;
1139 else {
1140 num[25]++;
1141 cout <<"Unknown particle: "<< pdg_code << endl;
1142 }
1143
1144 //if (pdg_code == kA10ID ) num[0]++; // not present
1145
1146 //Enter Momentums for Dalitz canidates
1147 //It would be prefered if we could get PreFSR Momentum
1148 if (pdg_code == kPiPlusID) piPlus.setVectM((*it)->initialFourMomentum().vect(), xmpion);
1149 if (pdg_code == kPiMinusID) piMinus.setVectM((*it)->initialFourMomentum().vect(), xmpion);
1150 if (pdg_code == kK0SID) k0.setVectM((*it)->initialFourMomentum().vect(), xmk0);
1151 if (pdg_code == kK0LID) k0.setVectM((*it)->initialFourMomentum().vect(), xmk0);
1152 }
1153 }
1154
1155 // D decay daughters have been tabulated. Now, classify decay.
1156 int nDaughters = 0 ;
1157 for( int i = 0 ; i < 26 ; i++ )
1158 {
1159 nDaughters += num[ i ] ;
1160 }
1161
1162 int nUFP0 = num[ kPi0 ] + num[ kEta ] + num[ kEtaPrime ] ;
1163 int nUFV0 = num[ kRho0 ] + num[ kPhi ] + num[ kOmega ] + num[ kGamma ] ;
1164 int nFV0 = num[ kKStar0 ] + num[ kK10 ] ;
1165 int nFV0Bar = num[ kKStar0Bar ] + num[ kK10Bar ] ;
1166 int nStrange = num[ kKMinus ] + num[ kFVMinus ] + nFV0Bar ;
1167 int nAntiStrange = num[ kKPlus ] + num[ kFVPlus ] + nFV0 ;
1168 int nCPPlusEig = num[ kNeutralScalar ] + num[ kRho0 ] + num[ kOmega ] + num[ kPhi ] + num[ kK0S ] + num[ kGamma ] ;
1169 int nCPMinusEig = num[ kPi0 ] + num[ kEta ] + num[ kEtaPrime ] + num[ kK0L ] ;
1170 int nCPEig = nCPPlusEig + nCPMinusEig ;
1171 int nChargedPiK = num[ kPiPlus ] + num[ kPiMinus ] + num[ kKPlus ] + num[ kKMinus ] ;
1172 int nK0 = num[ kK0S ] + num[ kK0L ] ;
1173
1174
1175 //check Dalitz modes: KsPi+Pi- or KlPi+Pi-
1176 double mnm_gen = 0. ;
1177 double mpn_gen = 0. ;
1178 double mpm_gen = 0. ;
1179 bool isKsPiPi = false ;
1180
1181 if( nK0 == 1 && num[ kPiPlus ] == 1 && num[ kPiMinus ] && nDaughters == 3 )
1182 {
1183 decay_mode = kDalitz ;
1184 k0.boost(-0.011, 0, 0);
1185 piMinus.boost(-0.011, 0, 0);
1186 piPlus.boost(-0.011, 0, 0);
1187 mnm_gen = (k0 + piMinus).m2();
1188 mpn_gen = (k0 + piPlus).m2();
1189 mpm_gen = (piPlus + piMinus).m2();
1190 if( num[ kK0S ] == 1) isKsPiPi = true ;
1191 }
1192
1193 // The order of if-statements below matters.
1194 if( decay_mode == kDalitz )
1195 {
1196 ++m_nDalitz ;
1197
1198 Dalitz t(8) ;
1199 TComplex D0, D0bar;
1200 if (m_dalitzModel == 0) { //Cleo model
1201 D0 = t.CLEO_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1202 D0bar = t.CLEO_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1203 if (m_dalitzModel == 1) { //Babar model
1204 D0 = t.Babar_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1205 D0bar = t.Babar_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1206 if (m_dalitzModel == 2) { //Belle model
1207 D0 = t.Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
1208 D0bar = t.Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
1209
1210 if( charm == 1){
1211 r2D0 = norm( D0bar ) / norm( D0 ) ;
1212 deltaD0 = arg( D0 * conj( D0bar ) ) + ( isKsPiPi ? PI : 0. );
1213
1214 double cosd = cos( deltaD0 ) ;
1215 double sind = sin( deltaD0 ) ;
1216 if( mpn_gen - mnm_gen > 0. )
1217 {
1218 m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd ;
1219 m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x ;
1220 m_dalitzDenom += -2. * r2D0 * cosd * cosd ;
1221 }
1222 }
1223 else {
1224 r2D0 = norm( D0 ) / norm( D0bar ) ;
1225 deltaD0 = arg( D0bar * conj( D0 ) ) + ( isKsPiPi ? PI : 0. ) ;
1226
1227 double cosd = cos( deltaD0 ) ;
1228 double sind = sin( deltaD0 ) ;
1229 if( mpn_gen - mnm_gen < 0. )
1230 {
1231 m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd ;
1232 m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x ;
1233 m_dalitzDenom += -2. * r2D0 * cosd * cosd ;
1234 }
1235 }
1236 }
1237 else if( num[ kLeptonPlus ] == 1 )
1238 {
1239 decay_mode = kSLPlus ;
1240 ++m_nSL ;
1241
1242 r2D0 = 0. ;
1243 deltaD0 = 0. ;
1244 }
1245 else if( num[ kLeptonMinus ] == 1 )
1246 {
1247 decay_mode = kSLMinus ;
1248 ++m_nSL ;
1249
1250 r2D0 = 0. ;
1251 deltaD0 = 0. ;
1252 }
1253 //Check CP even modes
1254 else if(
1255 ( nDaughters == 2 &&
1256 ( ( num[ kPiPlus ] == 1 && num[ kPiMinus ] == 1 ) ||
1257 nUFP0 == 2 ||
1258 num[ kNeutralScalar ] == 2 ||
1259 ( nUFP0 == 1 && nUFV0 == 1 ) ||
1260 ( num[ kNeutralScalar ] == 1 && nUFV0 == 1 ) ||
1261 ( num[ kKPlus ] == 1 && num[ kKMinus ] == 1 ) ||
1262 num[ kK0S ] == 2 ||
1263 num[ kK0L ] == 2 ||
1264 ( num[ kK0L ] == 1 && nUFP0 == 1 ) ||
1265 ( num[ kK0S ] == 1 && num[ kNeutralScalar ] == 1 ) ||
1266 ( num[ kK0L ] == 1 && nUFV0 == 1 ) ) ) ||
1267 ( nDaughters == 3 &&
1268 ( ( nCPPlusEig == 1 && num[ kPi0 ] == 2 ) ||
1269 ( num[ kK0S ] == 3 ) ) ) ||
1270 ( nDaughters == 4 &&
1271 num[ kK0L ] == 1 && num[ kPi0 ] == 3 )
1272 )
1273 {
1274 decay_mode = kCPPlus ;
1275 ++m_nCPPlus ;
1276
1277 r2D0 = 1. ;
1278 deltaD0 = PI;
1279 }
1280 //Check CP odd modes
1281 else if(
1282 ( nDaughters == 2 &&
1283 ( ( num[ kK0S ] == 1 && nUFP0 == 1 ) ||
1284 ( num[ kK0L ] == 1 && num[ kNeutralScalar ] == 1 ) ||
1285 ( num[ kK0S ] == 1 && nUFV0 == 1 ) ||
1286 ( nUFP0 == 1 && num[ kNeutralScalar ] == 1 ) ) ) ||
1287 ( nDaughters == 3 &&
1288 ( ( nCPMinusEig == 3 && num[ kPi0 ] == 2 ) || // pi0 is subset of CP-
1289 ( num[ kPi0 ] == 3 ) ||
1290 ( num[ kK0L ] == 3 ) ) ) ||
1291 ( nDaughters == 4 &&
1292 num[ kK0S ] == 1 && num[ kPi0 ] == 3 )
1293 )
1294 {
1295 decay_mode = kCPMinus ;
1296 ++m_nCPMinus ;
1297
1298 r2D0 = 1. ;
1299 deltaD0 = 0. ;
1300 }
1301 else if( nStrange == nAntiStrange + 1 )
1302 {
1303 decay_mode = kFlavoredCF ;
1304
1305 if( charm == 1 )
1306 {
1307 ++m_nFlavoredCFD0 ;
1308 r2D0 = pow(m_rCF,2) ;
1309 deltaD0 = m_deltaCF ;
1310 }
1311 else
1312 {
1313 ++m_nFlavoredDCSD0 ;
1314 r2D0 = 1. / pow( m_rCF,2 ) ;
1315 deltaD0 = -m_deltaCF ;
1316 }
1317 }
1318 else if( nAntiStrange == nStrange + 1 )
1319 {
1320 decay_mode = kFlavoredCFBar ;
1321
1322 if( charm == -1 )
1323 {
1324 ++m_nFlavoredCFD0 ;
1325 r2D0 = pow( m_rCF, 2) ;
1326 deltaD0 = m_deltaCF ;
1327 }
1328 else
1329 {
1330 ++m_nFlavoredDCSD0 ;
1331 r2D0 = 1. / pow( m_rCF, 2) ;
1332 deltaD0 = -m_deltaCF ;
1333 }
1334 }
1335 else if( nStrange == nAntiStrange ) // Unflavored or Cabibbo-suppressed
1336 {
1337 if( ( num[ kKPlus ] > 0 &&
1338 ( num[ kKPlus ] == num[ kFVMinus ] ||
1339 num[ kKPlus ] == nFV0Bar ) ) || // for anti-K*0 K+ pi-
1340 ( nK0 > 0 &&
1341 nFV0 != nFV0Bar &&
1342 nK0 == nFV0Bar ) ||
1343 ( num[ kPiPlus ] > 0 &&
1344 num[ kPiPlus ] == num[ kUFVMinus ] ) )
1345 {
1346 decay_mode = kFlavoredCS ;
1347 ++m_nFlavoredCSD0 ;
1348
1349 if( charm == 1 )
1350 {
1351 r2D0 = pow( m_rCS, 2 ) ;
1352 deltaD0 = m_deltaCS ;
1353 }
1354 else
1355 {
1356 r2D0 = 1. / pow( m_rCS, 2 ) ;
1357 deltaD0 = -m_deltaCS ;
1358 }
1359 }
1360 else if( ( num[ kKMinus ] > 0 && ( num[ kKMinus ] == num[ kFVPlus ] || num[ kKMinus ] == nFV0 ) ) || // for K*0 K- pi+
1361 ( nK0 > 0 && nFV0 != nFV0Bar && nK0 == nFV0 ) ||
1362 ( num[ kPiMinus ] > 0 && num[ kPiMinus ] == num[ kUFVPlus ] ) )
1363 {
1364 decay_mode = kFlavoredCSBar ;
1365 ++m_nFlavoredCSD0 ;
1366
1367 if( charm == -1 )
1368 {
1369 r2D0 = pow( m_rCS, 2 ) ;
1370 deltaD0 = m_deltaCS ;
1371 }
1372 else
1373 {
1374 r2D0 = 1. / pow( m_rCS, 2 ) ;
1375 deltaD0 = -m_deltaCS ;
1376 }
1377 }
1378
1379 // Self-conjugate mixed-CP final states -- pick CP or non-CP at
1380 // random. If CP, pick + or - at random. If non-CP, pick
1381 // flavored or flavoredBar according to the appropriate value of r.
1382 else if( ( nDaughters >= 3 && num[ kPiPlus ] == num[ kPiMinus ] &&
1383 num[ kKPlus ] == num[ kKMinus ] && nChargedPiK + nCPEig == nDaughters ) ||
1384 nUFV0 == nDaughters ||
1385 ( num[ kKStar0 ] > 0 && num[ kKStar0 ] == num[ kKStar0Bar ] ) ||
1386 ( num[ kK10 ] > 0 && num[ kK10 ] == num[ kK10Bar ] ) ||
1387 ( num[ kUFVPlus ] == 1 && num[ kUFVMinus ] == 1 )
1388 // for rho+rho-; no a1+a1- and rhoa1 final states in DECAY.DEC
1389 )
1390 {
1391 log << MSG::DEBUG <<" [ Self-conjugate mixed-CP ]" << endmsg ;
1392
1393 if( RandFlat::shoot() > 0.5 )
1394 {
1395 if( RandFlat::shoot() > 0.5 )
1396 {
1397 decay_mode = kCPPlus ;
1398 ++m_nCPPlus ;
1399
1400 r2D0 = 1. ;
1401 deltaD0 = PI;
1402 }
1403 else
1404 {
1405 decay_mode = kCPMinus ;
1406 ++m_nCPMinus ;
1407
1408 r2D0 = 1. ;
1409 deltaD0 = 0. ;
1410 }
1411 }
1412 else
1413 {
1414 // Odd # of K0S or K0L = CF; even # of K0S or K0L = SC.
1415 if( nK0 % 2 )
1416 {
1417 // Nov 2007: correct for r2 -> RWS and BR != A^2
1418 double cutoff = m_rwsCF / ( 1. + m_rwsCF ) ;
1419
1420 if( charm == -1 )
1421 {
1422 cutoff = 1. - cutoff ;
1423 }
1424
1425 log << MSG::DEBUG <<" [ cutoff = " << cutoff << " ]" << endmsg ;
1426
1427 decay_mode = ( RandFlat::shoot() > cutoff ) ? kFlavoredCF : kFlavoredCFBar ;
1428
1429 if( ( decay_mode == kFlavoredCF && charm == 1 ) ||
1430 ( decay_mode == kFlavoredCFBar && charm == -1 ) )
1431 {
1432 ++m_nFlavoredCFD0 ;
1433
1434 r2D0 = sqrt( m_rCF ) ;
1435 deltaD0 = m_deltaCF ;
1436 }
1437 else
1438 {
1439 ++m_nFlavoredDCSD0 ;
1440
1441 r2D0 = 1. / sqrt( m_rCF ) ;
1442 deltaD0 = -m_deltaCF ;
1443 }
1444 }
1445 else
1446 {
1447 // Nov 2007: correct for r2 -> RWS and BR != A^2
1448 double cutoff = m_rwsCS / ( 1. + m_rwsCS ) ;
1449
1450 if( charm == -1 )
1451 {
1452 cutoff = 1. - cutoff ;
1453 }
1454
1455 log << MSG::DEBUG <<" [ cutoff = " << cutoff << " ]" << endmsg ;
1456
1457 decay_mode = ( RandFlat::shoot() > cutoff ) ?
1458 kFlavoredCS : kFlavoredCSBar ;
1459 ++m_nFlavoredCSD0 ;
1460
1461 if( ( decay_mode == kFlavoredCS && charm == 1 ) ||
1462 ( decay_mode == kFlavoredCSBar && charm == -1 ) )
1463 {
1464 r2D0 = sqrt( m_rCS ) ;
1465 deltaD0 = m_deltaCS ;
1466 }
1467 else
1468 {
1469 r2D0 = 1. / sqrt( m_rCS ) ;
1470 deltaD0 = -m_deltaCS ;
1471 }
1472 }
1473 }
1474 }
1475 }
1476
1477 if (num[kUnknown] >= 1)
1478 {
1479 cout << "decay mode " << decay_mode << endl ;
1480 cout << "D #" << charm << endl ;
1481 cout << "pi+: " << num[ kPiPlus ] << endl ;
1482 cout << "pi-: " << num[ kPiMinus ] << endl ;
1483 cout << "pi0: " << num[ kPi0 ] << endl ;
1484 cout << "eta: " << num[ kEta ] << endl ;
1485 cout << "eta': " << num[ kEtaPrime ] << endl ;
1486 cout << "f0/a0: " << num[ kNeutralScalar ] << endl ;
1487 cout << "UFV+: " << num[ kUFVPlus ] << endl ;
1488 cout << "UFV-: " << num[ kUFVMinus ] << endl ;
1489 cout << "rho0: " << num[ kRho0 ] << endl ;
1490 cout << "omega: " << num[ kOmega ] << endl ;
1491 cout << "phi: " << num[ kPhi ] << endl ;
1492 cout << "K+: " << num[ kKPlus ] << endl ;
1493 cout << "K-: " << num[ kKMinus ] << endl ;
1494 cout << "K0S: " << num[ kK0S ] << endl ;
1495 cout << "K0L: " << num[ kK0L ] << endl ;
1496 cout << "FV+: " << num[ kFVPlus ] << endl ;
1497 cout << "FV-: " << num[ kFVMinus ] << endl ;
1498 cout << "K*0: " << num[ kKStar0 ] << endl ;
1499 cout << "K*0b: " << num[ kKStar0Bar ] << endl ;
1500 cout << "K10: " << num[ kK10 ] << endl ;
1501 cout << "K10b: " << num[ kK10Bar ] << endl ;
1502 cout << "l+: " << num[ kLeptonPlus ] << endl ;
1503 cout << "l-: " << num[ kLeptonMinus ] << endl ;
1504 cout << "nu: " << num[ kNeutrino ] << endl ;
1505 cout << "gamma: " << num[ kGamma ] << endl ;
1506 cout << "Unknown: " << num[ 25 ] << endl ;
1507 }
1508
1509
1510 d0_info[0]=decay_mode;
1511 d0_info[1]=r2D0;
1512 d0_info[2]=deltaD0;
1513 d0_info[3]=mnm_gen;
1514 d0_info[4]=mpn_gen;
1515 d0_info[5]= double (isKsPiPi);
1516 d0_info[6]=mpm_gen;
1517 return d0_info;
1518}
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
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition: Evt3Rank3C.cc:175
imag
double imag(const EvtComplex &c)
Definition: EvtComplex.hh:246
arg
double arg(const EvtComplex &c)
Definition: EvtComplex.hh:227
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
PI
const double PI
Definition: PipiJpsi.cxx:55
m_rwsCS
double m_rwsCS
Definition: QCMCFilter.cxx:167
xmpion
const double xmpion
Definition: QCMCFilter.cxx:56
m_nFlavoredCFD0
int m_nFlavoredCFD0
Definition: QCMCFilter.cxx:154
m_nCPPlus
int m_nCPPlus
Definition: QCMCFilter.cxx:152
dalitzNumer2_fil
double dalitzNumer2_fil
Definition: QCMCFilter.cxx:163
xmk0
const double xmk0
Definition: QCMCFilter.cxx:57
m_nCPMinus
int m_nCPMinus
Definition: QCMCFilter.cxx:153
m_deltaCF
double m_deltaCF
Definition: QCMCFilter.cxx:168
m_nUnknownEvents
int m_nUnknownEvents
Definition: QCMCFilter.cxx:145
m_dalitzNumer1
double m_dalitzNumer1
Definition: QCMCFilter.cxx:159
m_weights
HepSymMatrix m_weights(10, 0)
m_rwsCF
double m_rwsCF
Definition: QCMCFilter.cxx:166
m_nSL
int m_nSL
Definition: QCMCFilter.cxx:157
m_nDalitz
int m_nDalitz
Definition: QCMCFilter.cxx:158
m_modeCounter
HepMatrix m_modeCounter(10, 10, 0)
m_nD0D0barDiscarded
int m_nD0D0barDiscarded
Definition: QCMCFilter.cxx:150
m_nDpDmEvents
int m_nDpDmEvents
Definition: QCMCFilter.cxx:149
m_dalitzDenom
double m_dalitzDenom
Definition: QCMCFilter.cxx:161
xmd0
const double xmd0
Definition: QCMCFilter.cxx:59
m_rho_flag
int m_rho_flag
Definition: QCMCFilter.cxx:173
xmeta
const double xmeta
Definition: QCMCFilter.cxx:54
PI
const double PI
Definition: QCMCFilter.cxx:61
xmkaon
const double xmkaon
Definition: QCMCFilter.cxx:55
m_nD0bar
int m_nD0bar
Definition: QCMCFilter.cxx:148
m_keptModeCounter
HepMatrix m_keptModeCounter(10, 10, 0)
m_nDpDmDiscarded
int m_nDpDmDiscarded
Definition: QCMCFilter.cxx:151
xmpi0
const double xmpi0
Definition: QCMCFilter.cxx:53
m_dalitzNumer2
double m_dalitzNumer2
Definition: QCMCFilter.cxx:160
m_nFlavoredDCSD0
int m_nFlavoredDCSD0
Definition: QCMCFilter.cxx:156
dalitzDenom_fil
double dalitzDenom_fil
Definition: QCMCFilter.cxx:164
m_nFlavoredCSD0
int m_nFlavoredCSD0
Definition: QCMCFilter.cxx:155
m_nUnknownDecays
int m_nUnknownDecays
Definition: QCMCFilter.cxx:146
dalitzNumer1_fil
double dalitzNumer1_fil
Definition: QCMCFilter.cxx:162
m_nD0D0barEvents
int m_nD0D0barEvents
Definition: QCMCFilter.cxx:147
m_deltaCS
double m_deltaCS
Definition: QCMCFilter.cxx:169
xmrho
const double xmrho
Definition: QCMCFilter.cxx:58
nD0
int nD0
Definition: SD0Tag.cxx:53
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
t
TTree * t
Definition: binning.cxx:23
Dalitz
Definition: Dalitz.h:30
Event::McParticle::mother
const McParticle & mother() const
access to the mother particle
Definition: McParticle.cxx:117
Event::McParticle::particleProperty
StdHepId particleProperty() const
Retrieve particle property.
Definition: McParticle.cxx:7
QCMCFilter::initialize
StatusCode initialize()
Definition: QCMCFilter.cxx:197
QCMCFilter::finalize
StatusCode finalize()
Definition: QCMCFilter.cxx:793
QCMCFilter::QCMCFilter
QCMCFilter(const std::string &name, ISvcLocator *pSvcLocator)
Definition: QCMCFilter.cxx:178
QCMCFilter::findD0Decay
std::vector< double > findD0Decay(int charm)
Definition: QCMCFilter.cxx:1012
QCMCFilter::execute
StatusCode execute()
Definition: QCMCFilter.cxx:530
y
double y[1000]
Definition: draw_charge_space_coarse_II.cxx:8
x
double x[1000]
Definition: draw_charge_space_coarse_II.cxx:8
EventModel::MC::McParticleCol
_EXTERN_ std::string McParticleCol
Definition: EventModel.h:41