00001
00002
00003
00004 #include "QCMCFilterAlg/QCMCFilter.h"
00005 #include "QCMCFilterAlg/Dalitz.h"
00006
00007 #include "GaudiKernel/MsgStream.h"
00008 #include "GaudiKernel/AlgFactory.h"
00009 #include "GaudiKernel/ISvcLocator.h"
00010 #include "GaudiKernel/SmartDataPtr.h"
00011 #include "GaudiKernel/IDataProviderSvc.h"
00012 #include "GaudiKernel/PropertyMgr.h"
00013 #include "GaudiKernel/Bootstrap.h"
00014 #include "GaudiKernel/RegistryEntry.h"
00015
00016 #include "TMath.h"
00017
00018 #include <cmath>
00019 #include "HepPDT/ParticleDataTable.hh"
00020 #include "HepPDT/ParticleData.hh"
00021
00022 #include "GaudiKernel/IPartPropSvc.h"
00023
00024
00025 #include "McTruth/McParticle.h"
00026 #include "McTruth/MdcMcHit.h"
00027
00029 #include "EventModel/EventModel.h"
00030 #include "EventModel/Event.h"
00031 #include "EventModel/EventHeader.h"
00032
00033 #include "EvtRecEvent/EvtRecEvent.h"
00034
00035 #include "CLHEP/Random/RandFlat.h"
00036 #include "CLHEP/Matrix/Vector.h"
00037 #include "CLHEP/Matrix/Matrix.h"
00038 #include "CLHEP/Matrix/SymMatrix.h"
00039 #include "CLHEP/Vector/ThreeVector.h"
00040 #include "CLHEP/Vector/LorentzVector.h"
00041 #include "CLHEP/Vector/TwoVector.h"
00042 using CLHEP::HepVector;
00043 using CLHEP::Hep3Vector;
00044 using CLHEP::Hep2Vector;
00045 using CLHEP::HepLorentzVector;
00046
00047 #include <vector>
00048
00049
00050
00051
00052
00053 const double xmpi0 = 0.134976;
00054 const double xmeta = 0.54784;
00055 const double xmkaon = 0.49368;
00056 const double xmpion = 0.13957;
00057 const double xmk0 = 0.49761;
00058 const double xmrho = 0.77549;
00059 const double xmd0 = 1.86484;
00060
00061 const double PI = 3.1415926;
00062
00063
00064
00065 static const int kPsi3770ID = 30443 ;
00066 static const int kD0ID = 421 ;
00067 static const int kD0BarID = -421 ;
00068 static const int kDpID = 411 ;
00069 static const int kDmID = -411 ;
00070 static const int kGammaID = 22 ;
00071 static const int kGammaFSRID = -22 ;
00072 static const int kPiPlusID = 211 ;
00073 static const int kPiMinusID = -211 ;
00074 static const int kPi0ID = 111 ;
00075 static const int kEtaID = 221 ;
00076 static const int kEtaPrimeID = 331 ;
00077 static const int kF0ID = 9010221;
00078 static const int kFPrime0ID = 10221 ;
00079 static const int kF0m1500ID = 50221;
00080 static const int kF2ID = 225;
00081 static const int kA00ID = 10111 ;
00082 static const int kA0PlusID = 10211 ;
00083 static const int kA0MinusID = -10211 ;
00084 static const int kRhoPlusID = 213 ;
00085 static const int kRhoMinusID = -213 ;
00086 static const int kRho0ID = 113 ;
00087 static const int kRho2SPlusID = 30213 ;
00088 static const int kRho2SMinusID = -30213 ;
00089 static const int kRho2S0ID = 30113 ;
00090 static const int kA1PlusID = 20213 ;
00091 static const int kA1MinusID = -20213 ;
00092 static const int kA10ID = 20113 ;
00093 static const int kOmegaID = 223 ;
00094 static const int kPhiID = 333 ;
00095 static const int kKPlusID = 321 ;
00096 static const int kKMinusID = -321 ;
00097 static const int kK0SID = 310 ;
00098 static const int kK0LID = 130 ;
00099 static const int kK0ID = 311 ;
00100 static const int kK0BarID = -311 ;
00101 static const int kKStarPlusID = 323 ;
00102 static const int kKStarMinusID = -323 ;
00103 static const int kKStar0ID = 313 ;
00104 static const int kKStar0BarID = -313 ;
00105 static const int kK0Star0ID = 10311;
00106 static const int kK0Star0BarID = -10311;
00107 static const int kK0StarPlusID = 10321;
00108 static const int kK0StarMinusID = -10321;
00109 static const int kK1PlusID = 10323 ;
00110 static const int kK1MinusID = -10323 ;
00111 static const int kK10ID = 10313 ;
00112 static const int kK10BarID = -10313 ;
00113 static const int kK1PrimePlusID = 20323 ;
00114 static const int kK1PrimeMinusID = -20323 ;
00115 static const int kK1Prime0ID = 20313 ;
00116 static const int kK1Prime0BarID = -20313 ;
00117 static const int kK2StarPlusID = 325 ;
00118 static const int kK2StarMinusID = -325 ;
00119 static const int kK2Star0ID = 315 ;
00120 static const int kK2Star0BarID = -315 ;
00121 static const int kEMinusID = 11 ;
00122 static const int kEPlusID = -11 ;
00123 static const int kMuMinusID = 13 ;
00124 static const int kMuPlusID = -13 ;
00125 static const int kNuEID = 12 ;
00126 static const int kNuEBarID = -12 ;
00127 static const int kNuMuID = 14 ;
00128 static const int kNuMuBarID = -14 ;
00129
00130
00131 static const int kFlavoredCF = 0;
00132 static const int kFlavoredCFBar = 1;
00133 static const int kFlavoredCS = 2;
00134 static const int kFlavoredCSBar = 3;
00135 static const int kSLPlus = 4;
00136 static const int kSLMinus = 5;
00137 static const int kCPPlus = 6;
00138 static const int kCPMinus = 7;
00139 static const int kDalitz = 8;
00140 static const int kNDecayTypes = 9;
00141
00142
00143
00144
00145 int m_nUnknownEvents = 0;
00146 int m_nUnknownDecays= 0;
00147 int m_nD0D0barEvents= 0;
00148 int m_nD0bar = 0;
00149 int m_nDpDmEvents = 0;
00150 int m_nD0D0barDiscarded = 0;
00151 int m_nDpDmDiscarded = 0;
00152 int m_nCPPlus = 0;
00153 int m_nCPMinus = 0;
00154 int m_nFlavoredCFD0 = 0;
00155 int m_nFlavoredCSD0 = 0;
00156 int m_nFlavoredDCSD0 = 0;
00157 int m_nSL = 0;
00158 int m_nDalitz = 0;
00159 double m_dalitzNumer1 = 0;
00160 double m_dalitzNumer2 = 0;
00161 double m_dalitzDenom = 0;
00162 double dalitzNumer1_fil = 0;
00163 double dalitzNumer2_fil = 0;
00164 double dalitzDenom_fil = 0;
00165 HepSymMatrix m_weights(10,0);
00166 double m_rwsCF=0.;
00167 double m_rwsCS=0.;
00168 double m_deltaCF=0.;
00169 double m_deltaCS=0.;
00170 HepMatrix m_modeCounter(10,10,0);
00171 HepMatrix m_keptModeCounter(10,10,0);
00172
00173 int m_rho_flag;
00174
00175
00177
00178 QCMCFilter::QCMCFilter(const std::string& name, ISvcLocator* pSvcLocator) :
00179 Algorithm(name, pSvcLocator){
00180
00181 declareProperty("x", m_x=0.);
00182 declareProperty("y", m_y= 0.0);
00183 declareProperty("rForCFModes", m_rCF= 0.0621);
00184 declareProperty("zForCFModes", m_zCF= 2.0);
00185 declareProperty("wSignForCFModes", m_wCFSign= true);
00186 declareProperty("rForCSModes", m_rCS= 1.0);
00187 declareProperty("zForCSModes", m_zCS= -2.0);
00188 declareProperty("wSignForCSModes", m_wCSSign= true);
00189 declareProperty("DplusDminusWeight", m_dplusDminusWeight= 0.5);
00190 declareProperty("dalitzModel", m_dalitzModel= 0);
00191 declareProperty("UseNewWeights", m_useNewWeights= true);
00192 declareProperty("InvertSelection", m_invertSelection= false);
00193 }
00194
00195
00196
00197 StatusCode QCMCFilter::initialize() {
00198 MsgStream log(msgSvc(), name());
00199
00200 log << MSG::INFO << "in initialize()" << endmsg;
00201
00202
00203 IPartPropSvc* p_PartPropSvc;
00204 static const bool CREATEIFNOTTHERE(true);
00205 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
00206 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc)
00207 {
00208 log << MSG::ERROR << " Could not initialize Particle Properties Service" << endmsg;
00209 return PartPropStatus;
00210 }
00211 m_particleTable = p_PartPropSvc->PDT();
00212
00213
00214
00215
00216
00217
00218 m_y = (m_y - 0.008)/0.292;
00219
00220
00221
00222 double x = m_x ;
00223 double y = m_y ;
00224 double rCF = m_rCF ;
00225 double zCF = m_zCF ;
00226 double rCF2 = rCF * rCF ;
00227 double zCF2 = zCF * zCF ;
00228 double rCS = m_rCS ;
00229 double zCS = m_zCS ;
00230 double rCS2 = rCS * rCS ;
00231 double zCS2 = zCS * zCS ;
00232
00233 double rmix = ( x * x + y * y ) / 2. ;
00234 double wCF = ( m_wCFSign ? 1. : -1. ) * sqrt( 4. - zCF2 ) ;
00235 double wCS = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 ) ;
00236 double vCFCSPlus = ( zCF * zCS + wCF * wCS ) / 2. ;
00237 double vCFCSMinus = ( zCF * zCS - wCF * wCS ) / 2. ;
00238
00239 m_deltaCF = acos( zCF / 2. ) * ( m_wCFSign ? 1. : -1. ) ;
00240 m_deltaCS = acos( zCS / 2. ) * ( m_wCSSign ? 1. : -1. ) ;
00241
00242 double bCF = 1. + 0.5 * rCF * ( zCF * y + wCF * x ) + 0.5 * ( y * y - x * x );
00243 double bCS = 1. + 0.5 * rCS * ( zCS * y + wCS * x ) + 0.5 * ( y * y - x * x );
00244 double bCPPlus = 1. - y ;
00245 double bCPMinus = 1. + y ;
00246
00247 if( !m_useNewWeights )
00248 {
00249 bCF = 1. ;
00250 m_rwsCF = rCF2 ;
00251 bCS = 1. ;
00252 m_rwsCS = rCS2 ;
00253 bCPPlus = 1. ;
00254 bCPMinus = 1. ;
00255 }
00256
00257 m_rwsCF = ( rCF2 + 0.5 * rCF * ( zCF * y - wCF * x ) + rmix ) / bCF;
00258 m_rwsCS = ( rCS2 + 0.5 * rCS * ( zCS * y - wCS * x ) + rmix ) / bCS;
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268 m_weights[ kFlavoredCF ][ kFlavoredCFBar ] =
00269 1. + rCF2 * ( 2. - zCF2 ) + rCF2 * rCF2 ;
00270 m_weights[ kFlavoredCF ][ kFlavoredCF ] =
00271 rmix * m_weights[ kFlavoredCF ][ kFlavoredCFBar ] ;
00272 m_weights[ kFlavoredCF ][ kFlavoredCSBar ] =
00273 1. + rCF2 * rCS2 - rCF * rCS * vCFCSMinus ;
00274 m_weights[ kFlavoredCF ][ kFlavoredCS ] =
00275 rCF2 + rCS2 - rCF * rCS * vCFCSPlus +
00276 rmix * m_weights[ kFlavoredCF ][ kFlavoredCSBar ] ;
00277
00278
00279 m_weights[ kFlavoredCF ][ kFlavoredCF ] /= m_rwsCF * bCF * bCF ;
00280 m_weights[ kFlavoredCF ][ kFlavoredCFBar ] /=
00281 ( 1. + m_rwsCF * m_rwsCF ) * bCF * bCF ;
00282 m_weights[ kFlavoredCF ][ kFlavoredCS ] /=
00283 ( m_rwsCF + m_rwsCS ) * bCF * bCS ;
00284 m_weights[ kFlavoredCF ][ kFlavoredCSBar ] /=
00285 ( 1. + m_rwsCF * m_rwsCS ) * bCF * bCS ;
00286
00287 m_weights[ kFlavoredCF ][ kSLPlus ] = 1. / bCF ;
00288 m_weights[ kFlavoredCF ][ kSLMinus ] = 1. / bCF ;
00289
00290 m_weights[ kFlavoredCF ][ kCPPlus ] =
00291 ( 1. + rCF2 + rCF * zCF ) / ( 1. + m_rwsCF ) / bCF / bCPPlus ;
00292 m_weights[ kFlavoredCF ][ kCPMinus ] =
00293 ( 1. + rCF2 - rCF * zCF ) / ( 1. + m_rwsCF ) / bCF / bCPMinus ;
00294 m_weights[ kFlavoredCF ][ kDalitz ] = 1. / bCF;
00295
00296
00297 m_weights[ kFlavoredCFBar ][ kFlavoredCFBar ] = m_weights[ kFlavoredCF ][ kFlavoredCF ] ;
00298 m_weights[ kFlavoredCFBar ][ kFlavoredCS ] = m_weights[ kFlavoredCF ][ kFlavoredCSBar ] ;
00299 m_weights[ kFlavoredCFBar ][ kFlavoredCSBar ] = m_weights[ kFlavoredCF ][ kFlavoredCS ] ;
00300 m_weights[ kFlavoredCFBar ][ kSLPlus ] = 1. / bCF ;
00301 m_weights[ kFlavoredCFBar ][ kSLMinus ] = 1. / bCF ;
00302 m_weights[ kFlavoredCFBar ][ kCPPlus ] = m_weights[ kFlavoredCF ][ kCPPlus ] ;
00303 m_weights[ kFlavoredCFBar ][ kCPMinus ] = m_weights[ kFlavoredCF ][ kCPMinus ] ;
00304 m_weights[ kFlavoredCFBar ][ kDalitz ] = 1. / bCF;
00305
00306
00307
00308
00309 m_weights[ kFlavoredCS ][ kFlavoredCSBar ] = 1. + rCS2 * ( 2. - zCS2 ) + rCS2 * rCS2 ;
00310 m_weights[ kFlavoredCS ][ kFlavoredCS ] = rmix * m_weights[ kFlavoredCS ][ kFlavoredCSBar ] ;
00311
00312
00313 m_weights[ kFlavoredCS ][ kFlavoredCS ] /= m_rwsCS * bCS * bCS ;
00314 m_weights[ kFlavoredCS ][ kFlavoredCSBar ] /= ( 1. + m_rwsCS * m_rwsCS ) * bCS * bCS ;
00315
00316 m_weights[ kFlavoredCS ][ kSLPlus ] = 1. / bCS ;
00317 m_weights[ kFlavoredCS ][ kSLMinus ] = 1. / bCS ;
00318
00319 m_weights[ kFlavoredCS ][ kCPPlus ] = ( 1. + rCS2 + rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPPlus ;
00320 m_weights[ kFlavoredCS ][ kCPMinus ] = ( 1. + rCS2 - rCS * zCS ) / ( 1. + m_rwsCS ) / bCS / bCPMinus ;
00321
00322 m_weights[ kFlavoredCS ][ kDalitz ] = 1. / bCS;
00323
00324
00325
00326 m_weights[ kFlavoredCSBar ][ kFlavoredCSBar ] = m_weights[ kFlavoredCS ][ kFlavoredCS ] ;
00327 m_weights[ kFlavoredCSBar ][ kSLPlus ] = 1. / bCS ;
00328 m_weights[ kFlavoredCSBar ][ kSLMinus ] = 1. / bCS ;
00329 m_weights[ kFlavoredCSBar ][ kCPPlus ] = m_weights[ kFlavoredCS ][ kCPPlus ] ;
00330 m_weights[ kFlavoredCSBar ][ kCPMinus ] = m_weights[ kFlavoredCS ][ kCPMinus ] ;
00331 m_weights[ kFlavoredCSBar ][ kDalitz ] = 1. / bCS;
00332
00333
00334
00335
00336
00337 m_weights[ kSLPlus ][ kSLPlus ] = 0. ;
00338 m_weights[ kSLPlus ][ kSLMinus ] = 1. ;
00339 m_weights[ kSLPlus ][ kCPPlus ] = 1. / bCPPlus ;
00340 m_weights[ kSLPlus ][ kCPMinus ] = 1. / bCPMinus ;
00341 m_weights[ kSLPlus ][ kDalitz ] = 1. ;
00342
00343
00344
00345 m_weights[ kSLMinus ][ kSLMinus ] = 0. ;
00346 m_weights[ kSLMinus ][ kCPPlus ] = 1. / bCPPlus ;
00347 m_weights[ kSLMinus ][ kCPMinus ] = 1. / bCPMinus ;
00348 m_weights[ kSLMinus ][ kDalitz ] = 1. ;
00349
00350
00351
00352 m_weights[ kCPPlus ][ kCPPlus ] = 0. ;
00353 m_weights[ kCPPlus ][ kCPMinus ] = 2. / bCPPlus / bCPMinus ;
00354 m_weights[ kCPPlus ][ kDalitz ] = 1. / bCPPlus;
00355
00356
00357
00358 m_weights[ kCPMinus ][ kCPMinus ] = 0. ;
00359 m_weights[ kCPMinus ][ kDalitz ] = 1. / bCPMinus;
00360
00361
00362 m_weights[ kDalitz ][ kDalitz ] = 1;
00363
00364 log << MSG::FATAL << "Weight matrix" << m_weights << endmsg ;
00365
00366
00367 HepVector fractionsD0( kNDecayTypes+1, 0 ) ;
00368 fractionsD0[ kFlavoredCF ] = 0.531 ;
00369 fractionsD0[ kFlavoredCFBar ] = 0.0082 ;
00370 fractionsD0[ kFlavoredCS ] = 0.031 ;
00371 fractionsD0[ kFlavoredCSBar ] = 0.031 ;
00372 fractionsD0[ kSLPlus ] = 0.125 ;
00373 fractionsD0[ kSLMinus ] = 0. ;
00374 fractionsD0[ kCPPlus ] = 0.113 ;
00375 fractionsD0[ kCPMinus ] = 0.102 ;
00376 fractionsD0[ kDalitz ] = 0.0588 ;
00377 HepVector scales = m_weights * fractionsD0 ;
00378 log << MSG::INFO << "Integrated scale factors " <<scales << endmsg ;
00379
00380
00381 HepVector fractionsD0bar( kNDecayTypes+1, 0 ) ;
00382 fractionsD0bar[ kFlavoredCF ] = 0.0082 ;
00383 fractionsD0bar[ kFlavoredCFBar ] = 0.531 ;
00384 fractionsD0bar[ kFlavoredCS ] = 0.031 ;
00385 fractionsD0bar[ kFlavoredCSBar ] = 0.031 ;
00386 fractionsD0bar[ kSLPlus ] = 0. ;
00387 fractionsD0bar[ kSLMinus ] = 0.125 ;
00388 fractionsD0bar[ kCPPlus ] = 0.113 ;
00389 fractionsD0bar[ kCPMinus ] = 0.102 ;
00390 fractionsD0bar[ kDalitz ] = 0.0588 ;
00391 HepVector weight_norm = fractionsD0bar.T() * scales;
00392 log << MSG::INFO << "Overall scale factor " << fractionsD0bar.T() * scales << endmsg ;
00393
00394
00395
00396
00397
00398 double brCF = ( fractionsD0bar[ kFlavoredCFBar ] * scales[ kFlavoredCFBar ] ) / weight_norm[0] ;
00399 double brCS = ( fractionsD0bar[ kFlavoredCS ] * scales[ kFlavoredCS ] + fractionsD0bar[ kFlavoredCSBar ] * scales[ kFlavoredCSBar ] ) / weight_norm[0] ;
00400 double brDCS = ( fractionsD0bar[ kFlavoredCF ] * scales[ kFlavoredCF ] ) / weight_norm[0] ;
00401 double brCPPlus = ( fractionsD0bar[ kCPPlus ] * scales[ kCPPlus ] ) / weight_norm[0] ;
00402 double brCPMinus = ( fractionsD0bar[ kCPMinus ] * scales[ kCPMinus ] ) / weight_norm[0] ;
00403
00404
00405 double dalitz_y_num = -0.000177536;
00406 double dalitz_y_den = -0.0150846;
00407
00408 double numFactCF =
00409 -rCF * zCF * ( 1. - 0.5 * rCF * wCF * m_x ) ;
00410 double numFactCS =
00411 -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x ) ;
00412 double denFactCF = 0.5 * rCF2 * zCF2 ;
00413 double denFactCS = 0.5 * rCS2 * zCS2 ;
00414
00415 double num_pre =
00416 brCPPlus - brCPMinus + brCF * numFactCF + 0.5 * brCS * numFactCS + dalitz_y_num;
00417 double den_pre =
00418 1. - brCPPlus - brCPMinus - brCF * denFactCF - 0.5 * brCS * denFactCS + dalitz_y_den ;
00419 double y_pre = num_pre / den_pre ;
00420 double num =
00421 fractionsD0bar[ kCPPlus ] - fractionsD0bar[ kCPMinus ] + fractionsD0bar[ kFlavoredCFBar ] * numFactCF + fractionsD0bar[ kFlavoredCS ] * numFactCS + dalitz_y_num;
00422 double den =
00423 1. - fractionsD0bar[ kCPPlus ] - fractionsD0bar[ kCPMinus ] - fractionsD0bar[ kFlavoredCFBar ] * denFactCF - fractionsD0bar[ kFlavoredCS ] * denFactCS + dalitz_y_den;
00424 double y_prematrix = num / den ;
00425 double y_test1 = 0.25 * ( ( ( scales[ kCPMinus ] * m_weights[ kCPPlus ][ kSLPlus ] ) / ( scales[ kCPPlus ] * m_weights[ kCPMinus ][ kSLPlus ] ) ) -
00426 ( ( scales[ kCPPlus ] * m_weights[ kCPMinus ][ kSLPlus ] ) / ( scales[ kCPMinus ] * m_weights[ kCPPlus ][ kSLPlus ] ) ) ) ;
00427
00428 log << MSG::INFO << "An Input Y of " << m_y << endmsg ;
00429 log << MSG::INFO << "Parent MC has an intrinsic value of y as " << y_prematrix << endmsg ;
00430 log << MSG::INFO << "The BR method predics a y of " << y_pre << endmsg ;
00431 log << MSG::INFO << "The CP-SL double tag method predicts y of " <<y_test1 << endmsg ;
00432
00433
00434 m_largestWeight = 0. ;
00435 for( int i = 0 ; i < kNDecayTypes ; ++i )
00436 {
00437 for( int j = 0 ; j < kNDecayTypes ; ++j )
00438 {
00439 if( m_weights[ i ][ j ] < 0 ) m_weights[ i ][ j ] = 0;
00440 if( m_weights[ i ][ j ] > m_largestWeight )
00441 {
00442 m_largestWeight = m_weights[ i ][ j ] ;
00443 }
00444 }
00445 }
00446 m_weights /= m_largestWeight ;
00447
00448 log << MSG::INFO <<"Renormalized weight matrix " << m_weights << endmsg ;
00449
00450
00451
00452 Dalitz t(8);
00453 double D0Sum[ 10 ], D0barSum[ 10 ] ;
00454 TComplex D0D0barSum[ 10 ] ;
00455 int points[ 10 ] ;
00456
00457 for( int i = 0 ; i < 10 ; ++i )
00458 {
00459 D0Sum[ i ] = 0. ;
00460 D0barSum[ i ] = 0. ;
00461 D0D0barSum[ i ] = TComplex( 0., 0. ) ;
00462 points[ i ] = 0 ;
00463 }
00464
00465 double m2min = 0. ;
00466 double m2max = 3. ;
00467 int nsteps = 1000 ;
00468 double stepsize = ( m2max - m2min ) / ( double ) nsteps ;
00469
00470 for( int i = 0 ; i < nsteps ; ++i )
00471 {
00472 double m2i = m2min + ( ( double ) i + 0.5 ) * stepsize ;
00473
00474 for( int j = 0 ; j < nsteps ; ++j )
00475 {
00476 double m2j = m2min + ( ( double ) j + 0.5 ) * stepsize ;
00477
00478 if( t.Point_on_DP( m2i, m2j ) )
00479 {
00480 double m2k = 1.86450*1.86450 + 0.497648*0.497648 +
00481 0.139570*0.139570 + 0.139570*0.139570 - m2i - m2j ;
00482
00483 TComplex D0, D0bar;
00484 if (m_dalitzModel == 0) {
00485 D0 = t.CLEO_Amplitude( m2i, m2j, m2k ) ;
00486 D0bar = t.CLEO_Amplitude( m2j, m2i, m2k ) ;}
00487 if (m_dalitzModel == 1) {
00488 D0 = t.Babar_Amplitude( m2i, m2j, m2k ) ;
00489 D0bar = t.Babar_Amplitude( m2j, m2i, m2k ) ;}
00490 if (m_dalitzModel == 2) {
00491 D0 = t.Amplitude( m2i, m2j, m2k ) ;
00492 D0bar = t.Amplitude( m2j, m2i, m2k ) ;}
00493
00494 if ( j <= i ) {
00495 int bin = t.getBin( m2i, m2j, m2k ) ;
00496 if( bin == -1 ) bin = 8 ;
00497
00498 ++points[ bin ] ;
00499 D0Sum[ bin ] += norm( D0 ) ;
00500 D0barSum[ bin ] += norm( D0bar ) ;
00501 D0D0barSum[ bin ] += D0 * conj( D0bar ) ;
00502
00503 ++points[ 9 ] ;
00504 D0Sum[ 9 ] += norm( D0 ) ;
00505 D0barSum[ 9 ] += norm( D0bar ) ;
00506 D0D0barSum[ 9 ] += D0 * conj( D0bar ) ;
00507 }
00508
00509 }
00510 }
00511 }
00512
00513 for( int i = 0 ; i < 10 ; ++i )
00514 {
00515 double r2 = D0barSum[ i ] / D0Sum[ i ] ;
00516 double c = real( D0D0barSum[ i ] ) / sqrt( D0barSum[ i ] ) / sqrt( D0Sum[ i ] ) ;
00517 double s = imag( D0D0barSum[ i ] ) / sqrt( D0barSum[ i ] ) / sqrt( D0Sum[ i ] ) ;
00518
00519 cout << "BIN " << i << ": "
00520 << points[ i ] << " pts "
00521 << r2 << " " << c << " " << s << " " << c*c+s*s
00522 << endmsg ;
00523 }
00524
00525 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00526 return StatusCode::SUCCESS;
00527 }
00528
00529
00530 StatusCode QCMCFilter::execute() {
00531
00532
00533 ISvcLocator* svcLocator = Gaudi::svcLocator();
00534 StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
00535 if (sc.isFailure())
00536 std::cout<<"Could not access EventDataSvc!"<<std::endl;
00537
00538 setFilterPassed(false);
00539
00540 MsgStream log(msgSvc(), name());
00541 log<< MSG::INFO << "in execute()" << endmsg;
00542
00543 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00544
00545 if(!eventHeader)
00546 {
00547 cout<<" eventHeader "<<endl;
00548 return StatusCode::FAILURE;
00549 }
00550 long runNo = eventHeader->runNumber();
00551 long event = eventHeader->eventNumber();
00552 log << MSG::DEBUG << "run, evtnum = "
00553 << runNo << " , "
00554 << event << endmsg;
00555
00556
00557 SmartDataPtr<Event::McParticleCol> mcParticles(eventSvc(), EventModel::MC::McParticleCol);
00558
00559
00560 int psipp_daughter = 0;
00561 int d0_count = 0;
00562 int chd_count = 0;
00563 for( Event::McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end(); it++)
00564 {
00565 int pdg_code = (*it)->particleProperty();
00566 if ((*it)->primaryParticle()) continue;
00567 Event::McParticle it2 = (*it)->mother();
00568 int mother_pdg = 0;
00569
00570 mother_pdg = it2.particleProperty();
00571 if (mother_pdg == kPsi3770ID)
00572 {
00573 psipp_daughter++;
00574 if (abs(pdg_code) == kD0ID) d0_count++;
00575 if (abs(pdg_code) == kDpID) chd_count++;
00576 if (pdg_code == kD0BarID) m_nD0bar++;
00577 }
00578 }
00579
00580 if( psipp_daughter == 0 )
00581 {
00582 ++m_nUnknownEvents ;
00583
00584 if( m_invertSelection )
00585 {
00586 log<< MSG::DEBUG << "Discarding event -- did not find psi(3770). " << endmsg ;
00587 return StatusCode::SUCCESS;
00588 }
00589 else
00590 {
00591
00592 setFilterPassed(true);
00593
00594 return StatusCode::SUCCESS;
00595 }
00596 }
00597
00598 log<< MSG::DEBUG << "Found psi(3770) -->" << endmsg ;
00599
00600
00601 if ( psipp_daughter > 3 )
00602 {
00603 ++m_nUnknownEvents ;
00604
00605 if( m_invertSelection )
00606 {
00607 log<< MSG::DEBUG << "Discarding event -- psi(3770) has >3 daughters." << endmsg ;
00608 return StatusCode::SUCCESS;
00609 }
00610 else
00611 {
00612
00613 setFilterPassed(true);
00614
00615
00616 return StatusCode::SUCCESS;
00617 }
00618 }
00619
00620
00621 if (chd_count == 2) {
00622 ++m_nDpDmEvents ;
00623
00624
00625 double random = RandFlat::shoot() ;
00626
00627 if( ( random < m_dplusDminusWeight && !m_invertSelection ) ||
00628 ( random > m_dplusDminusWeight && m_invertSelection ) )
00629 {
00630
00631 setFilterPassed(true);
00632
00633
00634 return StatusCode::SUCCESS;
00635 }
00636 else
00637 {
00638 log<< MSG::DEBUG << "Discarding event -- D+D- event failed reweighting." << endmsg ;
00639 ++m_nDpDmDiscarded ;
00640 return StatusCode::SUCCESS;
00641 }
00642 }
00643
00644
00645 if (d0_count != 2) {
00646 if( !m_invertSelection )
00647 {
00648 log<< MSG::DEBUG << "Discarding event -- non-D0-D0bar event." << endmsg ;
00649 return StatusCode::SUCCESS;
00650 }
00651 else
00652 {
00653
00654 setFilterPassed(true);
00655
00656
00657 return StatusCode::SUCCESS;
00658 }
00659 }
00660
00661 log<< MSG::INFO << "D0-D0bar event." << endmsg ;
00662 m_nD0D0barEvents++ ;
00663
00664 m_rho_flag = 0;
00665
00666
00667 vector<double> d0_decay = findD0Decay(1);
00668 vector<double> d0bar_decay = findD0Decay(-1);
00669
00670 int d0mode = d0_decay[0];
00671 int d0barmode = d0bar_decay[0];
00672 m_modeCounter[d0mode][d0barmode]++;
00673
00674
00675
00676
00677 if (d0_decay[0] == 9 || d0bar_decay[0] == 9 ) {
00678 if( !m_invertSelection )
00679 {
00680 log<< MSG::DEBUG << "Discarding event -- unknown D0-D0bar event." << endmsg ;
00681 return StatusCode::SUCCESS;
00682 }
00683 else
00684 {
00685
00686 setFilterPassed(true);
00687
00688
00689
00690 return StatusCode::SUCCESS;
00691 }
00692 }
00693
00694
00695 double r2D0 = d0_decay[1] ;
00696 double deltaD0 = d0_decay[2] ;
00697 double r2D0bar = d0bar_decay[1] ;
00698 double deltaD0bar = d0bar_decay[2] ;
00699
00700
00701
00702 double weight ;
00703 if( d0_decay[0] == kDalitz || d0bar_decay[0] == kDalitz )
00704 {
00705 double r2prod = r2D0 * r2D0bar ;
00706 double x = m_x ;
00707 double y = m_y ;
00708
00709 double bD0 = 1. + sqrt( r2D0 ) * ( cos( deltaD0 ) * y + sin( deltaD0 ) * x ) ;
00710 double rwsD0 = ( r2D0 + sqrt( r2D0 ) * ( cos( deltaD0 ) * y - sin( deltaD0 ) * x ) ) / bD0 ;
00711 double bD0bar = 1. + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y + sin( deltaD0bar ) * x ) ;
00712 double rwsD0bar = ( r2D0bar + sqrt( r2D0bar ) * ( cos( deltaD0bar ) * y - sin( deltaD0bar ) * x ) ) / bD0bar ;
00713
00714 weight = 1. + r2prod - 2. * sqrt( r2prod ) * cos( deltaD0 + deltaD0bar ) ;
00715 weight /= ( 1. + rwsD0 * rwsD0bar ) * bD0 * bD0bar ;
00716 weight /= m_largestWeight ;
00717
00718
00719 int k = -10;
00720 double m2i= 0;
00721 double m2j= 0;
00722 double m2k= 0;
00723 if( d0_decay[0] == kDalitz) {
00724 k = d0bar_decay[0];
00725 m2i= d0_decay[3];
00726 m2j= d0_decay[4];
00727 m2k= d0_decay[6];}
00728 if( d0bar_decay[0] == kDalitz) {
00729 k = d0_decay[0];
00730 m2i= d0bar_decay[3];
00731 m2j= d0bar_decay[4];
00732 m2k= d0bar_decay[6];}
00733
00734 }
00735 else
00736 {
00737 weight = m_weights[ d0_decay[0] ][ d0bar_decay[0] ] ;
00738 }
00739
00740 double random = RandFlat::shoot() ;
00741 log<< MSG::DEBUG << "type: " << d0_decay[0] << " " << d0bar_decay[0] << ", weight " <<
00742 weight << " <? random " << random << endmsg ;
00743
00744 if( ( random < weight && !m_invertSelection ) || ( random > weight && m_invertSelection ) )
00745 {
00746
00747 setFilterPassed(true);
00748
00749
00750 int k = -10;
00751 double m2i= 0;
00752 double m2j= 0;
00753 double m2k= 0;
00754 double cosd = cos( deltaD0 ) ;
00755 double sind = sin( deltaD0 ) ;
00756 if( d0_decay[0] == kDalitz) {
00757 k = d0bar_decay[0];
00758 m2i= d0_decay[3];
00759 m2j= d0_decay[4];
00760 m2k= d0_decay[6];
00761 if ( m2j - m2i > 0. ) {
00762 dalitzNumer1_fil += -2. * sqrt( r2D0 ) * cosd ;
00763 dalitzNumer2_fil += 2. * r2D0 * cosd * sind * m_x ;
00764 dalitzDenom_fil += -2. * r2D0 * cosd * cosd ; }
00765 }
00766 if( d0bar_decay[0] == kDalitz) {
00767 k = d0_decay[0];
00768 m2i= d0bar_decay[3];
00769 m2j= d0bar_decay[4];
00770 m2k= d0bar_decay[6];
00771 cosd = cos( deltaD0bar ) ;
00772 sind = sin( deltaD0bar ) ;
00773 if( m2j - m2i < 0. ) {
00774 dalitzNumer1_fil += -2. * sqrt( r2D0bar ) * cosd ;
00775 dalitzNumer2_fil += 2. * r2D0bar * cosd * sind * m_x ;
00776 dalitzDenom_fil += -2. * r2D0bar * cosd * cosd ; }
00777 }
00778
00779 m_keptModeCounter[d0mode][d0barmode]++;
00780 return StatusCode::SUCCESS;
00781 }
00782 else
00783 {
00784 log << MSG::DEBUG << "Discarding event -- failed reweighting." << endmsg ;
00785
00786 ++m_nD0D0barDiscarded ;
00787 return StatusCode::SUCCESS;
00788 }
00789 }
00790
00791
00792
00793 StatusCode QCMCFilter::finalize() {
00794 MsgStream log(msgSvc(), name());
00795 log << MSG::INFO << "in finalize()" << endmsg;
00796 log << MSG::INFO << "Number of unknown events: " << m_nUnknownEvents << endmsg ;
00797 log << MSG::INFO << "Number of D+D- events seen: " << m_nDpDmEvents << endmsg ;
00798 log << MSG::INFO << "Number of D+D- events discarded: " << m_nDpDmDiscarded << endmsg ;
00799 if( m_nDpDmEvents > 0 )
00800 {
00801 log << MSG::INFO <<"Fraction discarded: " << double( m_nDpDmDiscarded ) / double( m_nDpDmEvents ) << endmsg ;
00802 }
00803
00804 log << MSG::INFO << "Number of D0D0bar events seen: " << m_nD0D0barEvents << " " << m_nD0bar << endmsg ;
00805 log << MSG::INFO << "Number of D0D0bar events discarded: " << m_nD0D0barDiscarded << endmsg ;
00806 if( m_nD0D0barEvents > 0 && m_nCPPlus > 0 && m_nCPMinus > 0)
00807 {
00808 log << MSG::INFO << "Fraction discarded: " << double( m_nD0D0barDiscarded ) / double( m_nD0D0barEvents ) << endl <<
00809 "Fraction unknown decays: " << 0.5 * double( m_nUnknownDecays ) / double( m_nD0D0barEvents ) <<endmsg ;
00810
00811 double nD0 = 2. * m_nD0D0barEvents ;
00812 double brSL = m_nSL / nD0 ;
00813 double dBrSL = sqrt( brSL * ( 1. - brSL ) / nD0 ) ;
00814 double brCF = m_nFlavoredCFD0 / nD0 ;
00815 double dBrCF = sqrt( brCF * ( 1. - brCF ) / nD0 ) ;
00816 double brCS = m_nFlavoredCSD0 / nD0 ;
00817 double dBrCS = sqrt( brCS * ( 1. - brCS ) / nD0 ) ;
00818 double brDCS = m_nFlavoredDCSD0 / nD0 ;
00819 double dBrDCS = sqrt( brDCS * ( 1. - brDCS ) / nD0 ) ;
00820 double brCPPlus = m_nCPPlus / nD0 ;
00821 double dBrCPPlus = sqrt( brCPPlus * ( 1. - brCPPlus ) / nD0 ) ;
00822 double brCPMinus = m_nCPMinus / nD0 ;
00823 double dBrCPMinus = sqrt( brCPMinus * ( 1. - brCPMinus ) / nD0 ) ;
00824 double brDalitz = m_nDalitz / nD0 ;
00825 double dBrDalitz = sqrt( brDalitz * ( 1. - brDalitz ) / nD0 ) ;
00826 double fdBrDalitz = dBrDalitz / brDalitz ;
00827 double dalitzNumer1Norm = m_dalitzNumer1 / nD0 ;
00828 double dDalitzNumer1Norm = dalitzNumer1Norm * fdBrDalitz ;
00829 double dalitzNumer2Norm = m_dalitzNumer2 / nD0 ;
00830 double dDalitzNumer2Norm = dalitzNumer2Norm * fdBrDalitz ;
00831 double dalitzDenomNorm = m_dalitzDenom / nD0 ;
00832 double dDalitzDenomNorm = dalitzDenomNorm * fdBrDalitz ;
00833
00834
00835 double fil_SL = 0, fil_FlavoredCFD0 = 0, fil_FlavoredCSD0 = 0, fil_FlavoredDCSD0 =0, fil_CPPlus = 0, fil_CPMinus = 0, fil_Dalitz = 0;
00836 for( int i = 0 ; i < 9 ; ++i )
00837 {
00838 fil_SL += m_keptModeCounter[i][4] + m_keptModeCounter[4][i] + m_keptModeCounter[i][5] + m_keptModeCounter[5][i] ;
00839 fil_FlavoredCFD0 += m_keptModeCounter[i][1] + m_keptModeCounter[0][i];
00840 fil_FlavoredCSD0 += m_keptModeCounter[i][2] + m_keptModeCounter[2][i] + m_keptModeCounter[i][3] + m_keptModeCounter[3][i];
00841 fil_FlavoredDCSD0 += m_keptModeCounter[i][0] + m_keptModeCounter[1][i];
00842 fil_CPPlus += m_keptModeCounter[i][6] + m_keptModeCounter[6][i];
00843 fil_CPMinus += m_keptModeCounter[i][7] + m_keptModeCounter[7][i];
00844 fil_Dalitz += m_keptModeCounter[i][8] + m_keptModeCounter[8][i];
00845 }
00846 double nD0_fil = 2. * ( m_nD0D0barEvents - m_nD0D0barDiscarded ) ;
00847 double brSL_fil = fil_SL / nD0_fil ;
00848 double dBrSL_fil = sqrt( brSL_fil * ( 1. - brSL_fil ) / nD0_fil ) ;
00849 double brCF_fil = fil_FlavoredCFD0 / nD0_fil ;
00850 double dBrCF_fil = sqrt( brCF_fil * ( 1. - brCF_fil ) / nD0_fil ) ;
00851 double brCS_fil = fil_FlavoredCSD0 / nD0_fil ;
00852 double dBrCS_fil = sqrt( brCS_fil * ( 1. - brCS_fil ) / nD0_fil ) ;
00853 double brDCS_fil = fil_FlavoredDCSD0 / nD0_fil ;
00854 double dBrDCS_fil = sqrt( brDCS_fil * ( 1. - brDCS_fil ) / nD0_fil ) ;
00855 double brCPPlus_fil = fil_CPPlus / nD0_fil ;
00856 double dBrCPPlus_fil = sqrt( brCPPlus_fil * ( 1. - brCPPlus_fil ) / nD0_fil ) ;
00857 double brCPMinus_fil = fil_CPMinus / nD0_fil ;
00858 double dBrCPMinus_fil = sqrt( brCPMinus_fil * ( 1. - brCPMinus_fil ) / nD0_fil ) ;
00859 double brDalitz_fil = fil_Dalitz / nD0_fil ;
00860 double dBrDalitz_fil = sqrt( brDalitz_fil * ( 1. - brDalitz_fil ) / nD0_fil ) ;
00861 double fdBrDalitz_fil = dBrDalitz_fil / brDalitz_fil ;
00862 double dalitzNumer1Norm_fil = dalitzNumer1_fil / nD0_fil ;
00863 double dDalitzNumer1Norm_fil = dalitzNumer1Norm_fil * fdBrDalitz_fil ;
00864 double dalitzNumer2Norm_fil = dalitzNumer2_fil / nD0_fil ;
00865 double dDalitzNumer2Norm_fil = dalitzNumer2Norm_fil * fdBrDalitz_fil ;
00866 double dalitzDenomNorm_fil = dalitzDenom_fil / nD0_fil ;
00867 double dDalitzDenomNorm_fil = dalitzDenomNorm_fil * fdBrDalitz_fil ;
00868
00869
00870 log << MSG::INFO <<
00871 "Original number of SL decays: " << m_nSL <<" ==> Br(SL) = " << brSL << " +- " << dBrSL << endl <<
00872 "Filtered number of SL decays: " << fil_SL <<" ==> Br(SL) = " << brSL_fil << " +- " << dBrSL_fil << endl <<
00873 "Original number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: " <<m_nFlavoredCFD0 << "/" << m_nFlavoredCSD0 << "/" << m_nFlavoredDCSD0 << endl <<
00874 " ==> Br(CF) = " << brCF << " +- " << dBrCF << endl <<
00875 " ==> Br(CS) = " << brCS << " +- " << dBrCS << endl <<
00876 " ==> Br(DCS) = " << brDCS << " +- " << dBrDCS << endl <<
00877 "Filtered number of D0->CF/CS/DCS or D0bar->CFbar/CSBar/DCSbar: " <<fil_FlavoredCFD0 << "/" << fil_FlavoredCSD0 << "/" << fil_FlavoredDCSD0 << endl <<
00878 " ==> Br(CF) = " << brCF_fil << " +- " << dBrCF_fil << endl <<
00879 " ==> Br(CS) = " << brCS_fil << " +- " << dBrCS_fil << endl <<
00880 " ==> Br(DCS) = " << brDCS_fil << " +- " << dBrDCS_fil << endl <<
00881
00882 "Original number of CP+/-: " << m_nCPPlus << "/" << m_nCPMinus <<endl <<
00883 " ==> Br(CP+) = " << brCPPlus << " +- " << dBrCPPlus <<endl <<
00884 " ==> Br(CP-) = " << brCPMinus << " +- " << dBrCPMinus << endl <<
00885 "Filtered number of CP+/-: " << fil_CPPlus << "/" << fil_CPMinus <<endl <<
00886 " ==> Br(CP+) = " << brCPPlus_fil << " +- " << dBrCPPlus_fil <<endl <<
00887 " ==> Br(CP-) = " << brCPMinus_fil << " +- " << dBrCPMinus_fil << endl <<
00888
00889 "Original number of Dalitz: " << m_nDalitz << endl <<
00890 " ==> Br(Dalitz) = " << brDalitz << " +- " << dBrDalitz <<endl <<
00891 " y contrib. numer1 " << dalitzNumer1Norm << " +- " << dDalitzNumer1Norm << endl <<
00892 " numer2 " << dalitzNumer2Norm << " +- " << dDalitzNumer2Norm << endl <<
00893 " denom " << dalitzDenomNorm << " +- " << dDalitzDenomNorm << endmsg <<
00894 "Filtered number of Dalitz: " << fil_Dalitz << endl <<
00895 " ==> Br(Dalitz) = " << brDalitz_fil << " +- " << dBrDalitz_fil <<endl <<
00896 " y contrib. numer1 " << dalitzNumer1Norm_fil << " +- " << dDalitzNumer1Norm_fil << endl <<
00897 " numer2 " << dalitzNumer2Norm_fil << " +- " << dDalitzNumer2Norm_fil << endl <<
00898 " denom " << dalitzDenomNorm_fil << " +- " << dDalitzDenomNorm_fil << endmsg ;
00899
00900 cout<<""<<endl;
00901
00902 HepMatrix differencetoWeight(10,10,0);
00903 for( int i = 0 ; i < 9 ; ++i ) {
00904 for( int j = 0 ; j < 9 ; ++j ) {
00905 if (m_modeCounter[i][j] != 0) {
00906 differencetoWeight[i][j] = m_keptModeCounter[i][j] / m_modeCounter[i][j] - m_weights[i][j]; } } }
00907
00908 log << MSG::INFO << "Weight matrix" << m_weights << endmsg ;
00909 log << MSG::INFO << "D0 Modes before filter" << m_modeCounter << endmsg ;
00910 log << MSG::INFO << "D0 Modes after filter" << m_keptModeCounter << endmsg ;
00911 log << MSG::INFO << "Compare percentage difference to weights " << differencetoWeight << endmsg ;
00912
00913
00914
00915 brCS /= 2. ;
00916 dBrCS /= 2. ;
00917 brCS_fil /= 2 ;
00918 dBrCS_fil /= 2 ;
00919
00920 double y, dy, y_fil, dy_fil ;
00921 double y_cpsl, dy_cpsl, y_fil_cpsl, dy_fil_cpsl ;
00922
00923
00924 double rCF = m_rCF ;
00925 double zCF = m_zCF ;
00926 double rCF2 = rCF * rCF ;
00927 double zCF2 = zCF * zCF ;
00928 double rCS = m_rCS ;
00929 double zCS = m_zCS ;
00930 double rCS2 = rCS * rCS ;
00931 double zCS2 = zCS * zCS ;
00932
00933 double wCF = ( m_wCFSign ? 1. : -1. ) * sqrt( 4. - zCF2 ) ;
00934 double wCS = ( m_wCSSign ? 1. : -1. ) * sqrt( 4. - zCS2 ) ;
00935
00936 double numFactCF =
00937 -rCF * zCF * ( 1. - 0.5 * rCF * wCF * m_x ) ;
00938 double numFactCS =
00939 -rCS * zCS * ( 1. - 0.5 * rCS * wCS * m_x ) ;
00940 double denFactCF = 0.5 * rCF2 * zCF2 ;
00941 double denFactCS = 0.5 * rCS2 * zCS2 ;
00942 double dalitzNumerNorm = dalitzNumer1Norm + dalitzNumer2Norm ;
00943
00944 double num =
00945 brCPPlus - brCPMinus + brCF * numFactCF + brCS * numFactCS +
00946 dalitzNumerNorm ;
00947 double den =
00948 1. - brCPPlus - brCPMinus - brCF * denFactCF - brCS * denFactCS +
00949 dalitzDenomNorm ;
00950 y = num / den ;
00951 dy = sqrt(
00952 ( num + den ) * ( num + den ) * dBrCPPlus * dBrCPPlus +
00953 ( num - den ) * ( num - den ) * dBrCPMinus * dBrCPMinus +
00954 ( numFactCF * den + denFactCF * num ) * dBrCF *
00955 ( numFactCF * den + denFactCF * num ) * dBrCF +
00956 ( numFactCS * den + denFactCS * num ) * dBrCS *
00957 ( numFactCS * den + denFactCS * num ) * dBrCS +
00958 ( dalitzNumerNorm * den + dalitzDenomNorm * num ) * fdBrDalitz *
00959 ( dalitzNumerNorm * den + dalitzDenomNorm * num ) * fdBrDalitz
00960 ) / ( den * den ) ;
00961
00962 double n_cpplussl = m_modeCounter[6][4] + m_modeCounter[6][5] + m_modeCounter[4][6] + m_modeCounter[5][6] ;
00963 double n_cpminussl = m_modeCounter[7][4] + m_modeCounter[7][5] + m_modeCounter[4][7] + m_modeCounter[5][7] ;
00964 double a_cpsl = ( n_cpplussl * m_nCPMinus ) / ( n_cpminussl * m_nCPPlus ) ;
00965 double b_cpsl = ( n_cpminussl * m_nCPPlus ) / ( n_cpplussl * m_nCPMinus ) ;
00966 y_cpsl = 0.25 * ( a_cpsl - b_cpsl ) ;
00967 dy_cpsl = 0.25 * ( a_cpsl + b_cpsl ) * sqrt( ( 1/n_cpplussl ) + ( 1/n_cpminussl ) + ( 1/m_nCPPlus ) + ( 1/m_nCPMinus ) ) ;
00968
00969
00970
00971 double dalitzNumerNorm_fil = dalitzNumer1Norm_fil + dalitzNumer2Norm_fil ;
00972 double num_fil =
00973 brCPPlus_fil - brCPMinus_fil + brCF_fil * numFactCF + brCS_fil * numFactCS +
00974 dalitzNumerNorm_fil ;
00975 double den_fil =
00976 1. - brCPPlus_fil - brCPMinus_fil - brCF_fil * denFactCF - brCS_fil * denFactCS +
00977 dalitzDenomNorm_fil ;
00978 y_fil = num_fil / den_fil ;
00979 dy_fil = sqrt(
00980 ( num_fil + den_fil ) * ( num_fil + den_fil ) * dBrCPPlus_fil * dBrCPPlus_fil +
00981 ( num_fil - den_fil ) * ( num_fil - den_fil ) * dBrCPMinus_fil * dBrCPMinus_fil +
00982 ( numFactCF * den_fil + denFactCF * num_fil ) * dBrCF_fil *
00983 ( numFactCF * den_fil + denFactCF * num_fil ) * dBrCF_fil +
00984 ( numFactCS * den_fil + denFactCS * num_fil ) * dBrCS_fil *
00985 ( numFactCS * den_fil + denFactCS * num_fil ) * dBrCS_fil +
00986 ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil * num_fil ) * fdBrDalitz_fil *
00987 ( dalitzNumerNorm_fil * den_fil + dalitzDenomNorm_fil * num_fil ) * fdBrDalitz_fil
00988 ) / ( den_fil * den_fil ) ;
00989
00990 double fil_cpplussl = m_keptModeCounter[6][4] + m_keptModeCounter[6][5] + m_keptModeCounter[4][6] + m_keptModeCounter[5][6] ;
00991 double fil_cpminussl = m_keptModeCounter[7][4] + m_keptModeCounter[7][5] + m_keptModeCounter[4][7] + m_keptModeCounter[5][7] ;
00992 a_cpsl = ( fil_cpplussl * fil_CPMinus ) / ( fil_cpminussl * fil_CPPlus ) ;
00993 b_cpsl = ( fil_cpminussl * fil_CPPlus ) / ( fil_cpplussl * fil_CPMinus ) ;
00994 y_fil_cpsl = 0.25 * ( a_cpsl - b_cpsl ) ;
00995 dy_fil_cpsl = 0.25 * ( a_cpsl + b_cpsl ) * sqrt( ( 1/n_cpplussl ) + ( 1/n_cpminussl ) + ( 1/fil_CPPlus ) + ( 1/fil_CPMinus ) ) ;
00996
00997
00998
00999 log << MSG::INFO <<"BR Method : Parent MC ==> y = " << y << " +- " << dy << endmsg ;
01000 log << MSG::INFO <<"BR Method : Filtered MC ==> y = " << y_fil << " +- " << dy_fil << endmsg ;
01001 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 ;
01002 log << MSG::INFO <<"CP-SL doubletag : Filtered MC ==> y = " << y_fil_cpsl << " +- " << dy_fil_cpsl << endmsg ;
01003
01004 }
01005
01006 return StatusCode::SUCCESS;
01007 }
01008
01009
01012 vector<double> QCMCFilter::findD0Decay(int charm)
01013 {
01014 MsgStream log(msgSvc(), name());
01015 log << MSG::INFO <<" In findD0Decay " << endmsg ;
01016
01017 vector<double> d0_info(7);
01018 for(int i=0; i< 7; i++) d0_info[i]=0;
01019 double decay_mode = 9;
01020 double r2D0 = -1;
01021 double deltaD0 = -1;
01022
01023 double num[26];
01024 for(int i=0; i< 26; i++) num[i]=0;
01025
01026 int kPiPlus = 0;
01027 int kPiMinus = 1;
01028 int kPi0 = 2;
01029 int kEta = 3;
01030 int kEtaPrime = 4;
01031 int kNeutralScalar = 5;
01032 int kUFVPlus = 6;
01033 int kUFVMinus = 7;
01034 int kRho0 = 8;
01035 int kOmega = 9;
01036 int kPhi = 10;
01037 int kKPlus = 11;
01038 int kKMinus = 12;
01039 int kK0S = 13;
01040 int kK0L = 14;
01041 int kFVPlus = 15;
01042 int kFVMinus = 16;
01043 int kKStar0 = 17;
01044 int kKStar0Bar = 18;
01045 int kK10 = 19;
01046 int kK10Bar = 20;
01047 int kLeptonPlus = 21;
01048 int kLeptonMinus = 22;
01049 int kNeutrino = 23;
01050 int kGamma = 24;
01051 int kUnknown = 25;
01052
01053
01054
01055 SmartDataPtr<Event::McParticleCol> mcParticles(eventSvc(), EventModel::MC::McParticleCol);
01056 int d0_pdg = charm == 1 ? kD0ID : kD0BarID;
01057
01058 HepLorentzVector piPlus, piMinus, k0;
01059
01060 for( Event::McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end(); it++)
01061 {
01062 int pdg_code = (*it)->particleProperty();
01063 if ((*it)->primaryParticle()) continue;
01064 Event::McParticle it2 = (*it)->mother();
01065 int mother_pdg = 0;
01066
01067
01068 mother_pdg = it2.particleProperty();
01069
01070
01071 if (pdg_code == kKStar0ID || pdg_code == kKStar0BarID || pdg_code == kK0ID || pdg_code == kK0BarID
01072 || pdg_code == kK10ID || pdg_code == kK10BarID
01073 || pdg_code == kK0Star0ID || pdg_code == kK0Star0BarID || pdg_code == kK0StarPlusID || pdg_code == kK0BarID ) continue;
01074
01075 if (mother_pdg == kK0ID || mother_pdg == kK0BarID) {
01076 it2 = it2.mother();
01077 mother_pdg = it2.particleProperty(); }
01078
01079 if (mother_pdg == kKStar0ID || mother_pdg == kKStar0BarID) {
01080
01081
01082 if (pdg_code == kPiPlusID || pdg_code == kPiMinusID) continue;
01083 if (mother_pdg == kKStar0ID && pdg_code == kKPlusID) pdg_code = kKStar0ID;
01084 if (mother_pdg == kKStar0BarID && pdg_code == kKMinusID) pdg_code = kKStar0BarID;
01085 it2 = it2.mother();
01086 mother_pdg = it2.particleProperty();
01087 }
01088
01089 if ( mother_pdg == kK0Star0ID || mother_pdg == kK0Star0BarID ) {
01090 it2 = it2.mother();
01091 mother_pdg = it2.particleProperty();
01092 }
01093
01094 if (mother_pdg == kK10ID || mother_pdg == kK10BarID) {
01095 it2 = it2.mother();
01096 mother_pdg = it2.particleProperty();
01097 }
01098
01099
01100
01101
01102 if (mother_pdg == d0_pdg)
01103 {
01104 if (pdg_code == kPiPlusID || pdg_code == kA0PlusID ) num[0]++;
01105 else if (pdg_code == kPiMinusID || pdg_code == kA0MinusID ) num[1]++;
01106 else if (pdg_code == kPi0ID ) num[2]++;
01107 else if (pdg_code == kEtaID ) num[3]++;
01108 else if (pdg_code == kEtaPrimeID ) num[4]++;
01109 else if (pdg_code == kF0ID || pdg_code == kA00ID
01110 || pdg_code == kFPrime0ID|| pdg_code == kF0m1500ID
01111 || pdg_code == kF2ID) num[5]++;
01112 else if (pdg_code == kRhoPlusID || pdg_code == kRho2SPlusID
01113 || pdg_code == kA1PlusID ) num[6]++;
01114 else if (pdg_code == kRhoMinusID || pdg_code == kRho2SMinusID
01115 || pdg_code == kA1MinusID) num[7]++;
01116 else if (pdg_code == kRho0ID || pdg_code == kRho2S0ID) num[8]++;
01117 else if (pdg_code == kOmegaID ) num[9]++;
01118 else if (pdg_code == kPhiID ) num[10]++;
01119 else if (pdg_code == kKPlusID ) num[11]++;
01120 else if (pdg_code == kKMinusID ) num[12]++;
01121 else if (pdg_code == kK0SID ) num[13]++;
01122 else if (pdg_code == kK0LID ) num[14]++;
01123 else if (pdg_code == kKStarPlusID || pdg_code == kK1PlusID
01124 || pdg_code == kK2StarPlusID || pdg_code == kK1PrimePlusID
01125 || pdg_code == kK0StarPlusID) num[15]++;
01126 else if (pdg_code == kKStarMinusID || pdg_code == kK1MinusID
01127 || pdg_code == kK2StarMinusID || pdg_code == kK1PrimeMinusID
01128 || pdg_code == kK0StarMinusID ) num[16]++;
01129 else if (pdg_code == kKStar0ID ) num[17]++;
01130 else if (pdg_code == kKStar0BarID ) num[18]++;
01131 else if (pdg_code == kK10ID || pdg_code == kK1Prime0ID) num[19]++;
01132 else if (pdg_code == kK10BarID || pdg_code == kK1Prime0BarID) num[20]++;
01133 else if (pdg_code == kEPlusID || pdg_code == kMuPlusID ) num[21]++;
01134 else if (pdg_code == kEMinusID || pdg_code == kMuMinusID ) num[22]++;
01135 else if (pdg_code == kNuEID || pdg_code == kNuEBarID
01136 || pdg_code == kNuMuID || pdg_code == kNuMuBarID ) num[23]++;
01137 else if (pdg_code == kGammaID ) num[24]++;
01138 else if (pdg_code == kGammaFSRID ) continue;
01139 else {
01140 num[25]++;
01141 cout <<"Unknown particle: "<< pdg_code << endl;
01142 }
01143
01144
01145
01146
01147
01148 if (pdg_code == kPiPlusID) piPlus.setVectM((*it)->initialFourMomentum().vect(), xmpion);
01149 if (pdg_code == kPiMinusID) piMinus.setVectM((*it)->initialFourMomentum().vect(), xmpion);
01150 if (pdg_code == kK0SID) k0.setVectM((*it)->initialFourMomentum().vect(), xmk0);
01151 if (pdg_code == kK0LID) k0.setVectM((*it)->initialFourMomentum().vect(), xmk0);
01152 }
01153 }
01154
01155
01156 int nDaughters = 0 ;
01157 for( int i = 0 ; i < 26 ; i++ )
01158 {
01159 nDaughters += num[ i ] ;
01160 }
01161
01162 int nUFP0 = num[ kPi0 ] + num[ kEta ] + num[ kEtaPrime ] ;
01163 int nUFV0 = num[ kRho0 ] + num[ kPhi ] + num[ kOmega ] + num[ kGamma ] ;
01164 int nFV0 = num[ kKStar0 ] + num[ kK10 ] ;
01165 int nFV0Bar = num[ kKStar0Bar ] + num[ kK10Bar ] ;
01166 int nStrange = num[ kKMinus ] + num[ kFVMinus ] + nFV0Bar ;
01167 int nAntiStrange = num[ kKPlus ] + num[ kFVPlus ] + nFV0 ;
01168 int nCPPlusEig = num[ kNeutralScalar ] + num[ kRho0 ] + num[ kOmega ] + num[ kPhi ] + num[ kK0S ] + num[ kGamma ] ;
01169 int nCPMinusEig = num[ kPi0 ] + num[ kEta ] + num[ kEtaPrime ] + num[ kK0L ] ;
01170 int nCPEig = nCPPlusEig + nCPMinusEig ;
01171 int nChargedPiK = num[ kPiPlus ] + num[ kPiMinus ] + num[ kKPlus ] + num[ kKMinus ] ;
01172 int nK0 = num[ kK0S ] + num[ kK0L ] ;
01173
01174
01175
01176 double mnm_gen = 0. ;
01177 double mpn_gen = 0. ;
01178 double mpm_gen = 0. ;
01179 bool isKsPiPi = false ;
01180
01181 if( nK0 == 1 && num[ kPiPlus ] == 1 && num[ kPiMinus ] && nDaughters == 3 )
01182 {
01183 decay_mode = kDalitz ;
01184 k0.boost(-0.011, 0, 0);
01185 piMinus.boost(-0.011, 0, 0);
01186 piPlus.boost(-0.011, 0, 0);
01187 mnm_gen = (k0 + piMinus).m2();
01188 mpn_gen = (k0 + piPlus).m2();
01189 mpm_gen = (piPlus + piMinus).m2();
01190 if( num[ kK0S ] == 1) isKsPiPi = true ;
01191 }
01192
01193
01194 if( decay_mode == kDalitz )
01195 {
01196 ++m_nDalitz ;
01197
01198 Dalitz t(8) ;
01199 TComplex D0, D0bar;
01200 if (m_dalitzModel == 0) {
01201 D0 = t.CLEO_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
01202 D0bar = t.CLEO_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
01203 if (m_dalitzModel == 1) {
01204 D0 = t.Babar_Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
01205 D0bar = t.Babar_Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
01206 if (m_dalitzModel == 2) {
01207 D0 = t.Amplitude( mpn_gen, mnm_gen, mpm_gen ) ;
01208 D0bar = t.Amplitude( mnm_gen, mpn_gen, mpm_gen ) ;}
01209
01210 if( charm == 1){
01211 r2D0 = norm( D0bar ) / norm( D0 ) ;
01212 deltaD0 = arg( D0 * conj( D0bar ) ) + ( isKsPiPi ? PI : 0. );
01213
01214 double cosd = cos( deltaD0 ) ;
01215 double sind = sin( deltaD0 ) ;
01216 if( mpn_gen - mnm_gen > 0. )
01217 {
01218 m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd ;
01219 m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x ;
01220 m_dalitzDenom += -2. * r2D0 * cosd * cosd ;
01221 }
01222 }
01223 else {
01224 r2D0 = norm( D0 ) / norm( D0bar ) ;
01225 deltaD0 = arg( D0bar * conj( D0 ) ) + ( isKsPiPi ? PI : 0. ) ;
01226
01227 double cosd = cos( deltaD0 ) ;
01228 double sind = sin( deltaD0 ) ;
01229 if( mpn_gen - mnm_gen < 0. )
01230 {
01231 m_dalitzNumer1 += -2. * sqrt( r2D0 ) * cosd ;
01232 m_dalitzNumer2 += 2. * r2D0 * cosd * sind * m_x ;
01233 m_dalitzDenom += -2. * r2D0 * cosd * cosd ;
01234 }
01235 }
01236 }
01237 else if( num[ kLeptonPlus ] == 1 )
01238 {
01239 decay_mode = kSLPlus ;
01240 ++m_nSL ;
01241
01242 r2D0 = 0. ;
01243 deltaD0 = 0. ;
01244 }
01245 else if( num[ kLeptonMinus ] == 1 )
01246 {
01247 decay_mode = kSLMinus ;
01248 ++m_nSL ;
01249
01250 r2D0 = 0. ;
01251 deltaD0 = 0. ;
01252 }
01253
01254 else if(
01255 ( nDaughters == 2 &&
01256 ( ( num[ kPiPlus ] == 1 && num[ kPiMinus ] == 1 ) ||
01257 nUFP0 == 2 ||
01258 num[ kNeutralScalar ] == 2 ||
01259 ( nUFP0 == 1 && nUFV0 == 1 ) ||
01260 ( num[ kNeutralScalar ] == 1 && nUFV0 == 1 ) ||
01261 ( num[ kKPlus ] == 1 && num[ kKMinus ] == 1 ) ||
01262 num[ kK0S ] == 2 ||
01263 num[ kK0L ] == 2 ||
01264 ( num[ kK0L ] == 1 && nUFP0 == 1 ) ||
01265 ( num[ kK0S ] == 1 && num[ kNeutralScalar ] == 1 ) ||
01266 ( num[ kK0L ] == 1 && nUFV0 == 1 ) ) ) ||
01267 ( nDaughters == 3 &&
01268 ( ( nCPPlusEig == 1 && num[ kPi0 ] == 2 ) ||
01269 ( num[ kK0S ] == 3 ) ) ) ||
01270 ( nDaughters == 4 &&
01271 num[ kK0L ] == 1 && num[ kPi0 ] == 3 )
01272 )
01273 {
01274 decay_mode = kCPPlus ;
01275 ++m_nCPPlus ;
01276
01277 r2D0 = 1. ;
01278 deltaD0 = PI;
01279 }
01280
01281 else if(
01282 ( nDaughters == 2 &&
01283 ( ( num[ kK0S ] == 1 && nUFP0 == 1 ) ||
01284 ( num[ kK0L ] == 1 && num[ kNeutralScalar ] == 1 ) ||
01285 ( num[ kK0S ] == 1 && nUFV0 == 1 ) ||
01286 ( nUFP0 == 1 && num[ kNeutralScalar ] == 1 ) ) ) ||
01287 ( nDaughters == 3 &&
01288 ( ( nCPMinusEig == 3 && num[ kPi0 ] == 2 ) ||
01289 ( num[ kPi0 ] == 3 ) ||
01290 ( num[ kK0L ] == 3 ) ) ) ||
01291 ( nDaughters == 4 &&
01292 num[ kK0S ] == 1 && num[ kPi0 ] == 3 )
01293 )
01294 {
01295 decay_mode = kCPMinus ;
01296 ++m_nCPMinus ;
01297
01298 r2D0 = 1. ;
01299 deltaD0 = 0. ;
01300 }
01301 else if( nStrange == nAntiStrange + 1 )
01302 {
01303 decay_mode = kFlavoredCF ;
01304
01305 if( charm == 1 )
01306 {
01307 ++m_nFlavoredCFD0 ;
01308 r2D0 = pow(m_rCF,2) ;
01309 deltaD0 = m_deltaCF ;
01310 }
01311 else
01312 {
01313 ++m_nFlavoredDCSD0 ;
01314 r2D0 = 1. / pow( m_rCF,2 ) ;
01315 deltaD0 = -m_deltaCF ;
01316 }
01317 }
01318 else if( nAntiStrange == nStrange + 1 )
01319 {
01320 decay_mode = kFlavoredCFBar ;
01321
01322 if( charm == -1 )
01323 {
01324 ++m_nFlavoredCFD0 ;
01325 r2D0 = pow( m_rCF, 2) ;
01326 deltaD0 = m_deltaCF ;
01327 }
01328 else
01329 {
01330 ++m_nFlavoredDCSD0 ;
01331 r2D0 = 1. / pow( m_rCF, 2) ;
01332 deltaD0 = -m_deltaCF ;
01333 }
01334 }
01335 else if( nStrange == nAntiStrange )
01336 {
01337 if( ( num[ kKPlus ] > 0 &&
01338 ( num[ kKPlus ] == num[ kFVMinus ] ||
01339 num[ kKPlus ] == nFV0Bar ) ) ||
01340 ( nK0 > 0 &&
01341 nFV0 != nFV0Bar &&
01342 nK0 == nFV0Bar ) ||
01343 ( num[ kPiPlus ] > 0 &&
01344 num[ kPiPlus ] == num[ kUFVMinus ] ) )
01345 {
01346 decay_mode = kFlavoredCS ;
01347 ++m_nFlavoredCSD0 ;
01348
01349 if( charm == 1 )
01350 {
01351 r2D0 = pow( m_rCS, 2 ) ;
01352 deltaD0 = m_deltaCS ;
01353 }
01354 else
01355 {
01356 r2D0 = 1. / pow( m_rCS, 2 ) ;
01357 deltaD0 = -m_deltaCS ;
01358 }
01359 }
01360 else if( ( num[ kKMinus ] > 0 && ( num[ kKMinus ] == num[ kFVPlus ] || num[ kKMinus ] == nFV0 ) ) ||
01361 ( nK0 > 0 && nFV0 != nFV0Bar && nK0 == nFV0 ) ||
01362 ( num[ kPiMinus ] > 0 && num[ kPiMinus ] == num[ kUFVPlus ] ) )
01363 {
01364 decay_mode = kFlavoredCSBar ;
01365 ++m_nFlavoredCSD0 ;
01366
01367 if( charm == -1 )
01368 {
01369 r2D0 = pow( m_rCS, 2 ) ;
01370 deltaD0 = m_deltaCS ;
01371 }
01372 else
01373 {
01374 r2D0 = 1. / pow( m_rCS, 2 ) ;
01375 deltaD0 = -m_deltaCS ;
01376 }
01377 }
01378
01379
01380
01381
01382 else if( ( nDaughters >= 3 && num[ kPiPlus ] == num[ kPiMinus ] &&
01383 num[ kKPlus ] == num[ kKMinus ] && nChargedPiK + nCPEig == nDaughters ) ||
01384 nUFV0 == nDaughters ||
01385 ( num[ kKStar0 ] > 0 && num[ kKStar0 ] == num[ kKStar0Bar ] ) ||
01386 ( num[ kK10 ] > 0 && num[ kK10 ] == num[ kK10Bar ] ) ||
01387 ( num[ kUFVPlus ] == 1 && num[ kUFVMinus ] == 1 )
01388
01389 )
01390 {
01391 log << MSG::DEBUG <<" [ Self-conjugate mixed-CP ]" << endmsg ;
01392
01393 if( RandFlat::shoot() > 0.5 )
01394 {
01395 if( RandFlat::shoot() > 0.5 )
01396 {
01397 decay_mode = kCPPlus ;
01398 ++m_nCPPlus ;
01399
01400 r2D0 = 1. ;
01401 deltaD0 = PI;
01402 }
01403 else
01404 {
01405 decay_mode = kCPMinus ;
01406 ++m_nCPMinus ;
01407
01408 r2D0 = 1. ;
01409 deltaD0 = 0. ;
01410 }
01411 }
01412 else
01413 {
01414
01415 if( nK0 % 2 )
01416 {
01417
01418 double cutoff = m_rwsCF / ( 1. + m_rwsCF ) ;
01419
01420 if( charm == -1 )
01421 {
01422 cutoff = 1. - cutoff ;
01423 }
01424
01425 log << MSG::DEBUG <<" [ cutoff = " << cutoff << " ]" << endmsg ;
01426
01427 decay_mode = ( RandFlat::shoot() > cutoff ) ? kFlavoredCF : kFlavoredCFBar ;
01428
01429 if( ( decay_mode == kFlavoredCF && charm == 1 ) ||
01430 ( decay_mode == kFlavoredCFBar && charm == -1 ) )
01431 {
01432 ++m_nFlavoredCFD0 ;
01433
01434 r2D0 = sqrt( m_rCF ) ;
01435 deltaD0 = m_deltaCF ;
01436 }
01437 else
01438 {
01439 ++m_nFlavoredDCSD0 ;
01440
01441 r2D0 = 1. / sqrt( m_rCF ) ;
01442 deltaD0 = -m_deltaCF ;
01443 }
01444 }
01445 else
01446 {
01447
01448 double cutoff = m_rwsCS / ( 1. + m_rwsCS ) ;
01449
01450 if( charm == -1 )
01451 {
01452 cutoff = 1. - cutoff ;
01453 }
01454
01455 log << MSG::DEBUG <<" [ cutoff = " << cutoff << " ]" << endmsg ;
01456
01457 decay_mode = ( RandFlat::shoot() > cutoff ) ?
01458 kFlavoredCS : kFlavoredCSBar ;
01459 ++m_nFlavoredCSD0 ;
01460
01461 if( ( decay_mode == kFlavoredCS && charm == 1 ) ||
01462 ( decay_mode == kFlavoredCSBar && charm == -1 ) )
01463 {
01464 r2D0 = sqrt( m_rCS ) ;
01465 deltaD0 = m_deltaCS ;
01466 }
01467 else
01468 {
01469 r2D0 = 1. / sqrt( m_rCS ) ;
01470 deltaD0 = -m_deltaCS ;
01471 }
01472 }
01473 }
01474 }
01475 }
01476
01477 if (num[kUnknown] >= 1)
01478 {
01479 cout << "decay mode " << decay_mode << endl ;
01480 cout << "D #" << charm << endl ;
01481 cout << "pi+: " << num[ kPiPlus ] << endl ;
01482 cout << "pi-: " << num[ kPiMinus ] << endl ;
01483 cout << "pi0: " << num[ kPi0 ] << endl ;
01484 cout << "eta: " << num[ kEta ] << endl ;
01485 cout << "eta': " << num[ kEtaPrime ] << endl ;
01486 cout << "f0/a0: " << num[ kNeutralScalar ] << endl ;
01487 cout << "UFV+: " << num[ kUFVPlus ] << endl ;
01488 cout << "UFV-: " << num[ kUFVMinus ] << endl ;
01489 cout << "rho0: " << num[ kRho0 ] << endl ;
01490 cout << "omega: " << num[ kOmega ] << endl ;
01491 cout << "phi: " << num[ kPhi ] << endl ;
01492 cout << "K+: " << num[ kKPlus ] << endl ;
01493 cout << "K-: " << num[ kKMinus ] << endl ;
01494 cout << "K0S: " << num[ kK0S ] << endl ;
01495 cout << "K0L: " << num[ kK0L ] << endl ;
01496 cout << "FV+: " << num[ kFVPlus ] << endl ;
01497 cout << "FV-: " << num[ kFVMinus ] << endl ;
01498 cout << "K*0: " << num[ kKStar0 ] << endl ;
01499 cout << "K*0b: " << num[ kKStar0Bar ] << endl ;
01500 cout << "K10: " << num[ kK10 ] << endl ;
01501 cout << "K10b: " << num[ kK10Bar ] << endl ;
01502 cout << "l+: " << num[ kLeptonPlus ] << endl ;
01503 cout << "l-: " << num[ kLeptonMinus ] << endl ;
01504 cout << "nu: " << num[ kNeutrino ] << endl ;
01505 cout << "gamma: " << num[ kGamma ] << endl ;
01506 cout << "Unknown: " << num[ 25 ] << endl ;
01507 }
01508
01509
01510 d0_info[0]=decay_mode;
01511 d0_info[1]=r2D0;
01512 d0_info[2]=deltaD0;
01513 d0_info[3]=mnm_gen;
01514 d0_info[4]=mpn_gen;
01515 d0_info[5]= double (isKsPiPi);
01516 d0_info[6]=mpm_gen;
01517 return d0_info;
01518 }