00001 #include "GaudiKernel/MsgStream.h"
00002 #include "GaudiKernel/AlgFactory.h"
00003 #include "GaudiKernel/ISvcLocator.h"
00004 #include "GaudiKernel/SmartDataPtr.h"
00005 #include "GaudiKernel/IDataProviderSvc.h"
00006 #include "GaudiKernel/PropertyMgr.h"
00007
00008 #include "EventModel/EventModel.h"
00009 #include "EventModel/Event.h"
00010 #include "EventModel/EventHeader.h"
00011
00012 #include "EvtRecEvent/EvtRecEvent.h"
00013 #include "EvtRecEvent/EvtRecTrack.h"
00014 #include "EvtRecEvent/EvtRecDTag.h"
00015 #include "EvtRecEvent/EvtRecVeeVertex.h"
00016 #include "EvtRecEvent/EvtRecPi0.h"
00017 #include "DstEvent/TofHitStatus.h"
00018
00019 #include "TMath.h"
00020 #include "GaudiKernel/INTupleSvc.h"
00021 #include "GaudiKernel/NTuple.h"
00022 #include "GaudiKernel/Bootstrap.h"
00023 #include "GaudiKernel/IHistogramSvc.h"
00024 #include "CLHEP/Vector/ThreeVector.h"
00025 #include "CLHEP/Vector/LorentzVector.h"
00026 #include "CLHEP/Vector/TwoVector.h"
00027 #include "CLHEP/Geometry/Point3D.h"
00028
00029 using CLHEP::Hep3Vector;
00030 using CLHEP::Hep2Vector;
00031 using CLHEP::HepLorentzVector;
00032 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00033 typedef HepGeom::Point3D<double> HepPoint3D;
00034 #endif
00035
00036 #include "DTagTool/DTagTool.h"
00037 #include "DDecayAlg/DDecay.h"
00038
00039 #include "VertexFit/KinematicFit.h"
00040 #include "VertexFit/VertexFit.h"
00041 #include "ParticleID/ParticleID.h"
00042
00043 #include <vector>
00044 typedef std::vector<int> Vint;
00045 typedef std::vector<HepLorentzVector> Vp4;
00046
00047
00049
00050 DDecay::DDecay(const std::string& name, ISvcLocator* pSvcLocator) :
00051 Algorithm(name, pSvcLocator){
00052
00053 }
00054
00055 StatusCode DDecay::initialize(){
00056
00057 MsgStream log(msgSvc(), name());
00058
00059 log << MSG::INFO << "in initialize()" << endmsg;
00060
00061 StatusCode status;
00062 NTuplePtr nt1(ntupleSvc(), "FILE1/vxyz");
00063 if ( nt1 ) m_tuple1 = nt1;
00064 else {
00065 m_tuple1 = ntupleSvc()->book ("FILE1/vxyz", CLID_ColumnWiseTuple, "track N-Tuple example");
00066 if ( m_tuple1 ) {
00067 status = m_tuple1->addItem ("vx0", m_vx0);
00068 status = m_tuple1->addItem ("vy0", m_vy0);
00069 status = m_tuple1->addItem ("vz0", m_vz0);
00070 status = m_tuple1->addItem ("vr0", m_vr0);
00071 }
00072 else {
00073 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00074 return StatusCode::FAILURE;
00075 }
00076 }
00077
00078 NTuplePtr nt2(ntupleSvc(), "FILE1/ks");
00079 if ( nt2 ) m_tuple2 = nt2;
00080 else {
00081 m_tuple2 = ntupleSvc()->book ("FILE1/ks", CLID_ColumnWiseTuple, "ks N-Tuple example");
00082 if ( m_tuple2 ) {
00083 status = m_tuple2->addItem ("ksmass", m_ksmass);
00084 status = m_tuple2->addItem ("ksd", m_ksd);
00085 status = m_tuple2->addItem ("ksmode", m_ksmode);
00086 }
00087 else {
00088 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
00089 return StatusCode::FAILURE;
00090 }
00091 }
00092
00093 NTuplePtr nt3(ntupleSvc(), "FILE1/pi0");
00094 if ( nt3 ) m_tuple3 = nt3;
00095 else {
00096 m_tuple3 = ntupleSvc()->book ("FILE1/pi0", CLID_ColumnWiseTuple, "pi0 N-Tuple example");
00097 if ( m_tuple3 ) {
00098 status = m_tuple3->addItem ("pi0mass", m_pi0mass);
00099 status = m_tuple3->addItem ("pi0mode", m_pi0mode);
00100 }
00101 else {
00102 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple3) << endmsg;
00103 return StatusCode::FAILURE;
00104 }
00105 }
00106
00107 NTuplePtr nt4(ntupleSvc(), "FILE1/tagD");
00108 if ( nt4 ) m_tuple4 = nt4;
00109 else {
00110 m_tuple4 = ntupleSvc()->book ("FILE1/tagD", CLID_ColumnWiseTuple, "DTag N-Tuple example");
00111 if ( m_tuple4 ) {
00112 status = m_tuple4->addItem ("mode", m_mode);
00113 status = m_tuple4->addItem ("type", m_type);
00114 status = m_tuple4->addItem ("charge", m_charge);
00115 status = m_tuple4->addItem ("charm", m_charm);
00116 status = m_tuple4->addItem ("numofchildren", m_numofchildren);
00117 status = m_tuple4->addItem ("mass", m_mass);
00118 status = m_tuple4->addItem ("mBC", m_mBC);
00119 status = m_tuple4->addItem ("deltaE", m_deltae);
00120 status = m_tuple4->addItem ("E", m_e);
00121 status = m_tuple4->addItem ("notherTrk", m_ntrk);
00122 }
00123 else {
00124 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple4) << endmsg;
00125 return StatusCode::FAILURE;
00126 }
00127 }
00128
00129
00130
00131
00132
00133
00134 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00135 return StatusCode::SUCCESS;
00136
00137 }
00138
00139
00140 StatusCode DDecay::execute() {
00141
00142 MsgStream log(msgSvc(), name());
00143 log << MSG::INFO << "in execute()" << endreq;
00144
00145 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00146 int runNo=eventHeader->runNumber();
00147 int eventNo=eventHeader->eventNumber();
00148
00149
00150
00151
00152
00154 SmartDataPtr<EvtRecVeeVertexCol> evtRecVeeVertexCol(eventSvc(), "/Event/EvtRec/EvtRecVeeVertexCol");
00155 if ( ! evtRecVeeVertexCol ) {
00156 log << MSG::FATAL << "Could not find EvtRecVeeVertexCol" << endreq;
00157 return StatusCode::FAILURE;
00158 }
00159
00161 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), "/Event/EvtRec/EvtRecPi0Col");
00162 if ( ! recPi0Col ) {
00163 log << MSG::FATAL << "Could not find EvtRecPi0Col" << endreq;
00164 return StatusCode::FAILURE;
00165 }
00166
00167
00168
00169 Hep3Vector xorigin(0,0,0);
00170 IVertexDbSvc* vtxsvc;
00171 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00172 if (vtxsvc->isVertexValid()) {
00173
00174
00175 double* vertex = vtxsvc->PrimaryVertex();
00176 xorigin.setX(vertex[0]);
00177 xorigin.setY(vertex[1]);
00178 xorigin.setZ(vertex[2]);
00179 }
00180
00181
00182
00183
00184
00185
00186 DTagTool dtagTool;
00187 if( dtagTool.isDTagListEmpty() ){
00188
00189 return StatusCode::SUCCESS;
00190 }
00191
00192 DTagToolIterator iter_begin = dtagTool.modes_begin();
00193 DTagToolIterator iter_end =dtagTool.modes_end();
00194
00195 cout<<"size of dtag:******:"<<iter_end-iter_begin<<endl;
00196
00197 int nCharge = 0;
00198
00199
00200
00201 cout<<"test single mode search :"<< EvtRecDTag::kD0toKPi<<endl;
00202
00203
00204 vector<int> mode = dtagTool.mode( EvtRecDTag::kD0toKPi );
00205 cout<<" there are "<< mode.size() <<" candidates for this mode" <<endl;
00206 for( int i=0; i < mode.size(); i++){
00207
00208 DTagToolIterator iter= dtagTool.modes_begin()+ mode[i];
00209 cout<<"No."<<i+1<<" candidate deltaE is : "<< (*iter)->deltaE()<<endl;
00210
00211 }
00212
00213
00214
00215
00216 for (DTagToolIterator iter_dtag=iter_begin; iter_dtag != iter_end; iter_dtag++){
00217
00218
00219 cout<<"***********"<<endl;
00220 cout<<"***********"<<endl;
00221 dtagTool<< iter_dtag;
00222
00223
00224
00225
00226
00227
00228 if((*iter_dtag)->decayMode()==EvtRecDTag::kD0toKPi) {
00229
00230 HepLorentzVector p4=(*iter_dtag)->p4();
00231 p4.boost(-0.011,0,0);
00232
00233 Hep3Vector p3=p4.v();
00234
00235 m_mode=(*iter_dtag)->decayMode();
00236 m_type=(*iter_dtag)->type();
00237 m_charge=(*iter_dtag)->charge();
00238 m_charm=(*iter_dtag)->charm();
00239 m_numofchildren=(*iter_dtag)->numOfChildren();
00240 m_mass=(*iter_dtag)->mass();
00241 m_mBC=(*iter_dtag)->mBC();
00242 m_e=(*iter_dtag)->beamE();
00243 m_deltae=(*iter_dtag)->deltaE();
00244
00245 SmartRefVector<EvtRecTrack> tracks=(*iter_dtag)->tracks();
00246 SmartRefVector<EvtRecTrack> othertracks=(*iter_dtag)->otherTracks();
00247 SmartRefVector<EvtRecTrack> othershowers=(*iter_dtag)->otherShowers();
00248 m_ntrk=othertracks.size();
00249
00250 m_tuple4->write();
00251
00252
00253 RecMdcKalTrack *mdcKalTrk1 = tracks[0]->mdcKalTrack();
00254 RecMdcKalTrack *mdcKalTrk2 = tracks[1]->mdcKalTrack();
00255 cout<<"same side track 1 charge is:"<<mdcKalTrk1->charge()<<endl;
00256 cout<<"same side track 2 charge is:"<<mdcKalTrk2->charge()<<endl;
00257
00258 for(int tk=0; tk<othertracks.size(); tk++){
00259 RecMdcTrack *mdcTrk = othertracks[tk]->mdcTrack();
00260 double pch=mdcTrk->p();
00261 double x0=mdcTrk->x();
00262 double y0=mdcTrk->y();
00263 double z0=mdcTrk->z();
00264 double phi0=mdcTrk->helix(1);
00265 double xp=xorigin.x();
00266 double yp=xorigin.y();
00267 double Rxy=(x0-xp)*cos(phi0)+(y0-yp)*sin(phi0);
00268
00269 m_vx0=x0;
00270 m_vy0=y0;
00271 m_vz0=z0;
00272 m_vr0=Rxy;
00273 m_tuple1->write();
00274 nCharge += mdcTrk->charge();
00275
00276 std::cout<<"other side track ID is: "<<othertracks[tk]->trackId()<<std::endl;
00277
00278 if(dtagTool.isPion(othertracks[tk]) )
00279 cout<<"it is pion"<<endl;
00280 if(dtagTool.isKaon(othertracks[tk]) )
00281 cout<<"it is kaon"<<endl;
00282
00283
00284
00285 }
00286
00287 for(int i=0; i<othershowers.size(); i++){
00288
00289
00290 }
00291
00292 }
00293
00294
00295
00296
00297
00298
00299 if((*iter_dtag)->decayMode()==EvtRecDTag::kD0toKPiPi0) {
00300
00301 m_mode=(*iter_dtag)->decayMode();
00302 m_type=(*iter_dtag)->type();
00303 m_charge=(*iter_dtag)->charge();
00304 m_charm=(*iter_dtag)->charm();
00305 m_numofchildren=(*iter_dtag)->numOfChildren();
00306 m_mass=(*iter_dtag)->mass();
00307 m_mBC=(*iter_dtag)->mBC();
00308 m_e=(*iter_dtag)->beamE();
00309 m_deltae=(*iter_dtag)->deltaE();
00310
00311 SmartRefVector<EvtRecTrack> othertracks=(*iter_dtag)->otherTracks();
00312 m_ntrk=othertracks.size();
00313
00314 m_tuple4->write();
00315
00316
00318
00319
00320 vector<int> pi0id= dtagTool.pi0Id(iter_dtag);
00321 cout<<"xxxxxxxxxxxxxxxxxxxxxxxxx"<<"num of pi0 is:"<<pi0id.size()<<endl;
00322
00323 for(int i=0; i<pi0id.size(); i++){
00324 pi0Iterator pi0Itr= dtagTool.pi0_begin()+pi0id[i];
00325 cout<<"pi0Mass: " << (*pi0Itr)->unconMass() << endl;
00326
00327 }
00328
00330 SmartRefVector<EvtRecTrack> showers=(*iter_dtag)->showers();
00331
00332 for(EvtRecPi0Col::iterator pi0Itr = recPi0Col->begin();
00333 pi0Itr < recPi0Col->end(); pi0Itr++){
00334
00336 EvtRecTrack* heGammaTrk = const_cast<EvtRecTrack*>((*pi0Itr)->hiEnGamma());
00337 EvtRecTrack* leGammaTrk = const_cast<EvtRecTrack*>((*pi0Itr)->loEnGamma());
00338
00339 int heGammaTrkId = heGammaTrk->trackId();
00340 int leGammaTrkId = leGammaTrk->trackId();
00341
00342 if((heGammaTrkId != showers[0]->trackId())&&
00343 (heGammaTrkId != showers[1]->trackId())) continue;
00344 if((leGammaTrkId != showers[0]->trackId())&&
00345 (leGammaTrkId != showers[1]->trackId())) continue;
00346
00347 const RecEmcShower* heGamma = heGammaTrk->emcShower();
00348 const RecEmcShower* leGamma = leGammaTrk->emcShower();
00349
00350 cout<<"pi0Mass: " << (*pi0Itr)->unconMass() << endl;
00351 cout<<" E(high): " << heGamma->energy() << endl;
00352 cout<<" E(low) : " << leGamma->energy() << endl;
00353
00354 m_pi0mass = (*pi0Itr)->unconMass();
00355 m_pi0mode = (*iter_dtag)->decayMode();
00356 m_tuple3->write();
00357
00359 break;
00360
00361 }
00362
00363
00364 }
00365
00366
00367
00368
00369
00370 if((*iter_dtag)->decayMode()==EvtRecDTag::kD0toKsPiPi) {
00371
00372 m_mode=(*iter_dtag)->decayMode();
00373 m_type=(*iter_dtag)->type();
00374 m_charge=(*iter_dtag)->charge();
00375 m_charm=(*iter_dtag)->charm();
00376 m_numofchildren=(*iter_dtag)->numOfChildren();
00377 m_mass=(*iter_dtag)->mass();
00378 m_mBC=(*iter_dtag)->mBC();
00379 m_e=(*iter_dtag)->beamE();
00380 m_deltae=(*iter_dtag)->deltaE();
00381
00382 SmartRefVector<EvtRecTrack> othertracks=(*iter_dtag)->otherTracks();
00383 m_ntrk=othertracks.size();
00384
00385 m_tuple4->write();
00386
00388
00389
00390 vector<int> ksid= dtagTool.ksId(iter_dtag);
00391 cout<<"xxxxxxxxxxxxxxxxxxxxxxxxx"<<"num of ks is:"<<ksid.size()<<endl;
00392
00393 for(int i=0; i<ksid.size(); i++){
00394 ksIterator ksItr= dtagTool.ks_begin()+ksid[i];
00395 cout<<"ksMass: " << (*ksItr)->mass() << endl;
00396
00397 }
00398
00399
00401 SmartRefVector<EvtRecTrack> tracks=(*iter_dtag)->tracks();
00402
00403 for(EvtRecVeeVertexCol::iterator ksItr = evtRecVeeVertexCol->begin();
00404 ksItr < evtRecVeeVertexCol->end(); ksItr++){
00405
00407 if((*ksItr)->vertexId() != 310) continue;
00408
00409 EvtRecTrack* aKsChild1Trk = (*ksItr)->daughter(0);
00410 EvtRecTrack* aKsChild2Trk = (*ksItr)->daughter(1);
00411
00412 int ksChild1TrkId = aKsChild1Trk->trackId();
00413 int ksChild2TrkId = aKsChild2Trk->trackId();
00414
00415 if((ksChild1TrkId != tracks[0]->trackId())&&
00416 (ksChild1TrkId != tracks[1]->trackId())) continue;
00417 if((ksChild2TrkId != tracks[0]->trackId())&&
00418 (ksChild2TrkId != tracks[1]->trackId())) continue;
00419
00420 cout<<"ksMass: " << (*ksItr)->mass() << endl;
00421
00422 Hep3Vector ks_D3(0,0,0);
00423 ks_D3.set((*ksItr)->w()[4],(*ksItr)->w()[5],(*ksItr)->w()[6]);
00424
00425 m_ksmass = (*ksItr)->mass();
00426 m_ksd = ks_D3.mag();
00427 m_ksmode = (*iter_dtag)->decayMode();
00428 m_tuple2->write();
00429
00431 break;
00432
00433 }
00434
00435 }
00436
00437
00438
00439
00440
00441 if((*iter_dtag)->decayMode()==EvtRecDTag::kDptoKPiPi) {
00442
00443 m_mode=(*iter_dtag)->decayMode();
00444 m_type=(*iter_dtag)->type();
00445 m_charge=(*iter_dtag)->charge();
00446 m_charm=(*iter_dtag)->charm();
00447 m_numofchildren=(*iter_dtag)->numOfChildren();
00448 m_mass=(*iter_dtag)->mass();
00449 m_mBC=(*iter_dtag)->mBC();
00450 m_e=(*iter_dtag)->beamE();
00451 m_deltae=(*iter_dtag)->deltaE();
00452
00453 SmartRefVector<EvtRecTrack> othertracks=(*iter_dtag)->otherTracks();
00454 m_ntrk=othertracks.size();
00455
00456 m_tuple4->write();
00457 }
00458
00459
00460
00461
00462 if((*iter_dtag)->decayMode()==EvtRecDTag::kDstoKsPi) {
00463
00464 m_mode=(*iter_dtag)->decayMode();
00465 m_type=(*iter_dtag)->type();
00466 m_charge=(*iter_dtag)->charge();
00467 m_charm=(*iter_dtag)->charm();
00468 m_numofchildren=(*iter_dtag)->numOfChildren();
00469 m_mass=(*iter_dtag)->mass();
00470 m_mBC=(*iter_dtag)->mBC();
00471 m_e=(*iter_dtag)->beamE();
00472 m_deltae=(*iter_dtag)->deltaE();
00473
00474 SmartRefVector<EvtRecTrack> othertracks=(*iter_dtag)->otherTracks();
00475 m_ntrk=othertracks.size();
00476
00477 m_tuple4->write();
00478 }
00479
00480
00481 }
00482
00483
00484
00485 cout<<"**************"<<endl;
00486 cout<<"**************"<<endl;
00487 cout<<"test print only D0/Dp/Ds modes"<<endl;
00488
00489
00490 vector<int> d0itindex= dtagTool.D0modes();
00491 for( int i=0; i< d0itindex.size(); i++){
00492 DTagToolIterator iter= dtagTool.modes_begin()+d0itindex[i];
00493 cout<<"No."<<i+1<<" D0 mode is :"<< (*iter)->decayMode()<<endl;
00494 }
00495
00496
00497 vector<int> dpitindex= dtagTool.Dpmodes();
00498 for( int i=0; i< dpitindex.size(); i++){
00499 DTagToolIterator iter= dtagTool.modes_begin()+dpitindex[i];
00500 cout<<"No."<<i+1<<" Dp mode is :"<< (*iter)->decayMode()<<endl;
00501 }
00502
00503
00504 vector<int> dsitindex= dtagTool.Dsmodes();
00505 for( int i=0; i< dsitindex.size(); i++){
00506 DTagToolIterator iter= dtagTool.modes_begin()+dsitindex[i];
00507 cout<<"No."<<i+1<<" Ds mode is :"<< (*iter)->decayMode()<<endl;
00508 }
00509
00510
00511
00512
00513 cout<<"**************"<<endl;
00514 cout<<"**************"<<endl;
00515 cout<<"test single tag "<<endl;
00516
00517 if( dtagTool.findSTag( EvtRecDTag::kD0toKPi , +1 ) )
00518 cout<<" find single tag mode: "<< (*dtagTool.stag())->decayMode() <<endl;
00519
00520
00521 if( dtagTool.findSTag( EvtRecDTag::kD0toKPi , +1 ) && dtagTool.cosmicandleptonVeto() ){
00522 cout<<"cosmic and lepton backgaround veto"<<endl;
00523 cout<<" find single tag mode: "<< (*dtagTool.stag())->decayMode() <<endl;
00524 }
00525
00526
00527
00528 cout<<"**************"<<endl;
00529 cout<<"**************"<<endl;
00530 cout<<"test double tag "<<endl;
00531
00532
00533
00534
00535 bool dtagflag=dtagTool.findDTag( EvtRecDTag::kD0toKPi , EvtRecDTag::kD0toKPiPi0);
00536 if( dtagflag){
00537 cout<<" find double tag mode 1 :"<<endl;
00538 dtagTool<<dtagTool.dtag1();
00539 cout<<" find double tag mode 2: "<< endl;
00540 dtagTool<<dtagTool.dtag2();
00541 }
00542
00543
00544 dtagflag=dtagTool.findDTag( EvtRecDTag::kD0toKPi , +1, EvtRecDTag::kD0toKPiPi0, -1 );
00545 if( dtagflag){
00546 cout<<" find double tag mode 1 :"<<endl;
00547 dtagTool<<dtagTool.dtag1();
00548 cout<<" find double tag mode 2: "<< endl;
00549 dtagTool<<dtagTool.dtag2();
00550 }
00551
00552
00553
00554 dtagflag=dtagTool.findDTag( EvtRecDTag::kD0toKPi , +1, EvtRecDTag::kD0toKPiPi0, -1 ,"inv");
00555 if( dtagflag){
00556 cout<<" double tag by invariant mass:"<<endl;
00557 cout<<" find double tag mode 1 :"<<endl;
00558 dtagTool<<dtagTool.dtag1();
00559 cout<<" find double tag mode 2: "<< endl;
00560 dtagTool<<dtagTool.dtag2();
00561 }
00562
00563
00564 dtagflag=dtagTool.findADTag( EvtRecDTag::kD0toKPi , +1, EvtRecDTag::kD0toKPiPi0, -1);
00565 if( dtagflag){
00566
00567 vector<DTagToolIterator> vdtag1=dtagTool.vdtag1();
00568 vector<DTagToolIterator> vdtag2=dtagTool.vdtag2();
00569
00570 cout<<" list of all doule tags:"<<endl;
00571 for(int i=0;i<vdtag1.size();i++){
00572 cout<<" find double tag mode 1 :"<<endl;
00573 dtagTool<<vdtag1[i];
00574 cout<<" find double tag mode 2: "<< endl;
00575 dtagTool<<vdtag2[i];
00576 }
00577
00578 }
00579
00580 dtagTool.clear();
00581
00582 }
00583
00584
00585
00586
00587 StatusCode DDecay::finalize() {
00588 MsgStream log(msgSvc(), name());
00589 log << MSG::INFO << "in finalize()" << endmsg;
00590 return StatusCode::SUCCESS;
00591 }
00592