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 #include "VertexFit/IVertexDbSvc.h"
00008 #include "GaudiKernel/Bootstrap.h"
00009 #include "GaudiKernel/ISvcLocator.h"
00010
00011 #include "EventModel/EventModel.h"
00012 #include "EventModel/Event.h"
00013
00014 #include "EvtRecEvent/EvtRecEvent.h"
00015 #include "EvtRecEvent/EvtRecTrack.h"
00016 #include "DstEvent/TofHitStatus.h"
00017 #include "EventModel/EventHeader.h"
00018
00019
00020
00021 #include "TMath.h"
00022 #include "GaudiKernel/INTupleSvc.h"
00023 #include "GaudiKernel/NTuple.h"
00024 #include "GaudiKernel/Bootstrap.h"
00025 #include "GaudiKernel/IHistogramSvc.h"
00026 #include "CLHEP/Vector/ThreeVector.h"
00027 #include "CLHEP/Vector/LorentzVector.h"
00028 #include "CLHEP/Vector/TwoVector.h"
00029 using CLHEP::Hep3Vector;
00030 using CLHEP::Hep2Vector;
00031 using CLHEP::HepLorentzVector;
00032 #include "CLHEP/Geometry/Point3D.h"
00033 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00034 typedef HepGeom::Point3D<double> HepPoint3D;
00035 #endif
00036 #include "PpjrhopiAlg/Ppjrhopi.h"
00037
00038 #include "VertexFit/KinematicFit.h"
00039 #include "VertexFit/VertexFit.h"
00040 #include "ParticleID/ParticleID.h"
00041
00042 #include <vector>
00043
00044
00045 const double mpi = 0.13957;
00046 const double mk = 0.493677;
00047 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00048
00049 const double velc = 299.792458;
00050 typedef std::vector<int> Vint;
00051 typedef std::vector<HepLorentzVector> Vp4;
00052
00053 int Ncut0,Ncut1,Ncut2,Ncut3,Ncut4,Ncut5,Ncut6,Ncut7,Ncut8,Ncut9,Ncut10;
00054
00056
00057 Ppjrhopi::Ppjrhopi(const std::string& name, ISvcLocator* pSvcLocator) :
00058 Algorithm(name, pSvcLocator) {
00059
00060
00061 declareProperty("Vr0cut", m_vr0cut=5.0);
00062 declareProperty("Vz0cut", m_vz0cut=20.0);
00063 declareProperty("Vr1cut", m_vr1cut=1.0);
00064 declareProperty("Vz1cut", m_vz1cut=5.0);
00065 declareProperty("Vctcut", m_cthcut=0.93);
00066 declareProperty("EnergyThreshold", m_energyThreshold=0.04);
00067 declareProperty("GammaAngCut", m_gammaAngCut=20.0);
00068 declareProperty("Test4C", m_test4C = 1);
00069 declareProperty("Test5C", m_test5C = 1);
00070 declareProperty("CheckDedx", m_checkDedx = 1);
00071 declareProperty("CheckTof", m_checkTof = 1);
00072 }
00073
00074
00075 StatusCode Ppjrhopi::initialize(){
00076 MsgStream log(msgSvc(), name());
00077
00078 log << MSG::INFO << "in initialize()" << endmsg;
00079
00080 StatusCode status;
00081
00082 if(m_test4C==1) {
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286 }
00287
00288
00289
00290
00291
00292 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00293 return StatusCode::SUCCESS;
00294
00295 }
00296
00297
00298 StatusCode Ppjrhopi::execute() {
00299
00300
00301
00302 MsgStream log(msgSvc(), name());
00303 log << MSG::INFO << "in execute()" << endreq;
00304
00305 setFilterPassed(false);
00306
00307 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00308 int runNo=eventHeader->runNumber();
00309 int event=eventHeader->eventNumber();
00310 log << MSG::DEBUG <<"runNo, evtnum = "
00311 << runNo << " , "
00312 << event <<endreq;
00313
00314 Ncut0++;
00315
00316
00317 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00318 log << MSG::INFO << "get event tag OK" << endreq;
00319 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00320 << evtRecEvent->totalCharged() << " , "
00321 << evtRecEvent->totalNeutral() << " , "
00322 << evtRecEvent->totalTracks() <<endreq;
00323
00324
00325
00326 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00327
00328
00329
00330
00331 if(evtRecEvent->totalNeutral()>100) {
00332 return StatusCode::SUCCESS;
00333 }
00334
00335 Vint iGood, ipip, ipim;
00336 iGood.clear();
00337 ipip.clear();
00338 ipim.clear();
00339 Vp4 ppip, ppim;
00340 ppip.clear();
00341 ppim.clear();
00342
00343 Hep3Vector xorigin(0,0,0);
00344
00345
00346 IVertexDbSvc* vtxsvc;
00347 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00348 if(vtxsvc->isVertexValid()){
00349 double* dbv = vtxsvc->PrimaryVertex();
00350 double* vv = vtxsvc->SigmaPrimaryVertex();
00351
00352
00353 xorigin.setX(dbv[0]);
00354 xorigin.setY(dbv[1]);
00355 xorigin.setZ(dbv[2]);
00356 }
00357
00358 int nCharge = 0;
00359 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00360 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00361 if(!(*itTrk)->isMdcTrackValid()) continue;
00362 if (!(*itTrk)->isMdcKalTrackValid()) continue;
00363 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00364
00365 double pch =mdcTrk->p();
00366 double x0 =mdcTrk->x();
00367 double y0 =mdcTrk->y();
00368 double z0 =mdcTrk->z();
00369 double phi0=mdcTrk->helix(1);
00370 double xv=xorigin.x();
00371 double yv=xorigin.y();
00372 double Rxy=fabs((x0-xv)*cos(phi0)+(y0-yv)*sin(phi0));
00373
00374 double m_vx0 = x0;
00375 double m_vy0 = y0;
00376 double m_vz0 = z0-xorigin.z();
00377 double m_vr0 = Rxy;
00378 double m_Vctc=z0/sqrt(Rxy*Rxy+z0*z0);
00379 double m_Vct =cos(mdcTrk->theta());
00380
00381
00382
00383 if(fabs(m_vz0) >= m_vz0cut) continue;
00384 if(m_vr0 >= m_vr0cut) continue;
00385
00386 iGood.push_back((*itTrk)->trackId());
00387 nCharge += mdcTrk->charge();
00388 }
00389
00390
00391
00392
00393 int nGood = iGood.size();
00394
00395 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
00396 if((nGood != 4)||(nCharge!=0)){
00397 return StatusCode::SUCCESS;
00398 }
00399
00400 Ncut1++;
00401
00402 Vint iGam;
00403 iGam.clear();
00404 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
00405 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00406 if(!(*itTrk)->isEmcShowerValid()) continue;
00407 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00408 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
00409
00410 double dthe = 200.;
00411 double dphi = 200.;
00412 double dang = 200.;
00413 log << MSG::DEBUG << "liuf neu= " <<i <<endreq;
00414 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
00415 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
00416 if(!(*jtTrk)->isExtTrackValid()) continue;
00417 RecExtTrack *extTrk = (*jtTrk)->extTrack();
00418 if(extTrk->emcVolumeNumber() == -1) continue;
00419 Hep3Vector extpos = extTrk->emcPosition();
00420 log << MSG::DEBUG << "liuf charge= " <<j <<endreq;
00421
00422 double angd = extpos.angle(emcpos);
00423 double thed = extpos.theta() - emcpos.theta();
00424 double phid = extpos.deltaPhi(emcpos);
00425 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00426 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00427
00428
00429
00430 if(angd < dang) {
00431 dang = angd;
00432 dthe = thed;
00433 dphi = phid;
00434 }
00435 }
00436 if(dang>=200) continue;
00437 double eraw = emcTrk->energy();
00438 dthe = dthe * 180 / (CLHEP::pi);
00439 dphi = dphi * 180 / (CLHEP::pi);
00440 dang = dang * 180 / (CLHEP::pi);
00441
00442 double m_dthe = dthe;
00443 double m_dphi = dphi;
00444 double m_dang = dang;
00445 double m_eraw = eraw;
00446
00447
00448 log << MSG::DEBUG << "eraw dang= " << eraw << " , " <<dang <<"," <<i <<endreq;
00449 if(eraw < m_energyThreshold) continue;
00450 if(dang < m_gammaAngCut) continue;
00451
00452
00453
00454
00455 iGam.push_back((*itTrk)->trackId());
00456 }
00457
00458
00459
00460
00461 int nGam = iGam.size();
00462
00463 log << MSG::DEBUG << "num Good Photon " << nGam << " , " <<evtRecEvent->totalNeutral()<<endreq;
00464 if(nGam<2){
00465 return StatusCode::SUCCESS;
00466 }
00467
00468 Ncut2++;
00469
00470
00471
00472
00473
00474
00475
00476
00477 Vp4 pGam;
00478 pGam.clear();
00479 for(int i = 0; i < nGam; i++) {
00480 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
00481 RecEmcShower* emcTrk = (*itTrk)->emcShower();
00482 double eraw = emcTrk->energy();
00483 double phi = emcTrk->phi();
00484 double the = emcTrk->theta();
00485 HepLorentzVector ptrk;
00486 ptrk.setPx(eraw*sin(the)*cos(phi));
00487 ptrk.setPy(eraw*sin(the)*sin(phi));
00488 ptrk.setPz(eraw*cos(the));
00489 ptrk.setE(eraw);
00490
00491
00492
00493 pGam.push_back(ptrk);
00494 }
00495
00496
00497 for(int i = 0; i < nGood; i++) {
00498 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00499 if (!(*itTrk)->isMdcTrackValid()) continue;
00500 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00501 if (!(*itTrk)->isMdcKalTrackValid()) continue;
00502 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
00503 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
00504 if(mdcKalTrk->charge() >0 ) {
00505 ipip.push_back(iGood[i]);
00506 HepLorentzVector ptrk;
00507 ptrk.setPx(mdcKalTrk->px());
00508 ptrk.setPy(mdcKalTrk->py());
00509 ptrk.setPz(mdcKalTrk->pz());
00510 double p3 = ptrk.mag();
00511 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00512 ppip.push_back(ptrk);
00513 } else {
00514 ipim.push_back(iGood[i]);
00515 HepLorentzVector ptrk;
00516 ptrk.setPx(mdcKalTrk->px());
00517 ptrk.setPy(mdcKalTrk->py());
00518 ptrk.setPz(mdcKalTrk->pz());
00519 double p3 = ptrk.mag();
00520 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00521 ppim.push_back(ptrk);
00522 }
00523 }
00524
00525
00526 int npip = ipip.size();
00527 int npim = ipim.size();
00528 if(npip!=2||npim != 2) return SUCCESS;
00529
00530
00531
00532
00533 Ncut3++;
00534
00535
00536
00537
00538
00539
00540
00541
00542 HepLorentzVector pTot0(0.011*3.6862,0,0,3.6862);
00543 HepLorentzVector pTrec1,pTrec2,pTrec3,pTrec4;
00544 HepLorentzVector pTrecf;
00545 double m_recjpsi1,m_recjpsi2,m_recjpsi3,m_recjpsi4,m_recppf;
00546 double deljp1,deljp2,deljp3,deljp4;
00547 pTrec1 = pTot0 - ppip[0] - ppim[0];
00548 pTrec2 = pTot0 - ppip[0] - ppim[1];
00549 pTrec3 = pTot0 - ppip[1] - ppim[0];
00550 pTrec4 = pTot0 - ppip[1] - ppim[1];
00551 m_recjpsi1 = pTrec1.m();
00552 m_recjpsi2 = pTrec2.m();
00553 m_recjpsi3 = pTrec3.m();
00554 m_recjpsi4 = pTrec4.m();
00555 deljp1=fabs(m_recjpsi1-3.097);
00556 deljp2=fabs(m_recjpsi2-3.097);
00557 deljp3=fabs(m_recjpsi3-3.097);
00558 deljp4=fabs(m_recjpsi4-3.097);
00559
00560 int itmp,itmp1,itmp2;
00561 HepLorentzVector ptmp,ptmp1,ptmp2;
00562
00563 pTrecf =pTrec1;
00564 m_recppf=pTrec1.m();
00565
00566 if(deljp2<deljp1&&deljp2<deljp3&&deljp2<deljp4)
00567 { itmp= ipim[1];
00568 ipim[1]=ipim[0];
00569 ipim[0]=itmp;
00570
00571 ptmp =ppim[1];
00572 ppim[1]=ppim[0];
00573 ppim[0]=ptmp;
00574
00575 pTrecf =pTrec2;
00576 m_recppf=pTrec2.m();
00577 }
00578
00579 if(deljp3<deljp1&&deljp3<deljp2&&deljp3<deljp4)
00580 { itmp= ipip[1];
00581 ipip[1]=ipip[0];
00582 ipip[0]=itmp;
00583
00584 ptmp =ppip[1];
00585 ppip[1]=ppip[0];
00586 ppip[0]=ptmp;
00587
00588 pTrecf =pTrec3;
00589 m_recppf=pTrec3.m();
00590 }
00591
00592 if(deljp4<deljp1&&deljp4<deljp2&&deljp4<deljp3)
00593 { itmp1= ipip[1];
00594 ipip[1]=ipip[0];
00595 ipip[0]=itmp1;
00596 itmp2= ipim[1];
00597 ipim[1]=ipim[0];
00598 ipim[0]=itmp2;
00599
00600 ptmp1 =ppip[1];
00601 ppip[1]=ppip[0];
00602 ppip[0]=ptmp1;
00603 ptmp2 =ppim[1];
00604 ppim[1]=ppim[0];
00605 ppim[0]=ptmp2;
00606
00607 pTrecf =pTrec4;
00608 m_recppf=pTrec4.m();
00609 }
00610
00611 if(fabs(m_recppf-3.097)>0.2) return SUCCESS;
00612
00613 log << MSG::DEBUG << "ngood track ID after jpsi = " << ipip[0] << " , "
00614 << ipim[0]<< " , " << ipip[1] << " , " << ipim[1] << endreq;
00615 Ncut4++;
00616
00617 HepLorentzVector ppi2_no1 = ppip[0] + ppim[0];
00618 HepLorentzVector ppi2_no2 = ppip[1] + ppim[1];
00619 HepLorentzVector ppi2_no3 = ppip[0] + ppim[1];
00620 HepLorentzVector ppi2_no4 = ppip[1] + ppim[0];
00621 HepLorentzVector p4pi_no = ppi2_no1+ ppi2_no2;
00622
00623 double emcTg1=0.0;
00624 double emcTg2=0.0;
00625 double emcTg3=0.0;
00626 double emcTg4=0.0;
00627 double laypi1=-1.0;
00628 double laypi2=-1.0;
00629 double laypi3=-1.0;
00630 double laypi4=-1.0;
00631
00632 EvtRecTrackIterator itTrkp1=evtRecTrkCol->begin() + ipip[0];
00633 RecMdcTrack* mdcTrkp1 = (*itTrkp1)->mdcTrack();
00634 RecMdcKalTrack *mdcKalTrkp1 = (*itTrkp1)->mdcKalTrack();
00635 RecEmcShower* emcTrkp1 = (*itTrkp1)->emcShower();
00636 RecMucTrack *mucTrkp1=(*itTrkp1)->mucTrack();
00637
00638 double phi01=mdcTrkp1->helix(1);
00639 double m_p1vx = mdcTrkp1->x();
00640 double m_p1vy = mdcTrkp1->y();
00641 double m_p1vz = mdcTrkp1->z()-xorigin.z();
00642 double m_p1vr = fabs((mdcTrkp1->x()-xorigin.x())*cos(phi01)+(mdcTrkp1->y()-xorigin.y())*sin(phi01));
00643 double m_p1vct=cos(mdcTrkp1->theta());
00644 double m_p1ptot=mdcKalTrkp1->p();
00645 double m_p1pxy=sqrt(mdcKalTrkp1->px()*mdcKalTrkp1->px()+mdcKalTrkp1->py()*mdcKalTrkp1->py());
00646
00647 if((*itTrkp1)->isEmcShowerValid()){
00648 emcTg1=emcTrkp1->energy();
00649 }
00650 if((*itTrkp1)->isMucTrackValid()){
00651 laypi1=mucTrkp1->numLayers();
00652 }
00653 double m_laypip1=laypi1;
00654
00655 EvtRecTrackIterator itTrkm1=evtRecTrkCol->begin() + ipim[0];
00656 RecMdcTrack* mdcTrkm1 = (*itTrkm1)->mdcTrack();
00657 RecMdcKalTrack *mdcKalTrkm1 = (*itTrkm1)->mdcKalTrack();
00658 RecEmcShower* emcTrkm1 = (*itTrkm1)->emcShower();
00659 RecMucTrack *mucTrkm1=(*itTrkm1)->mucTrack();
00660
00661 double phi02=mdcTrkm1->helix(1);
00662 double m_m1vx = mdcTrkm1->x();
00663 double m_m1vy = mdcTrkm1->y();
00664 double m_m1vz = mdcTrkm1->z()-xorigin.z();
00665 double m_m1vr = fabs((mdcTrkm1->x()-xorigin.x())*cos(phi02)+(mdcTrkm1->y()-xorigin.y())*sin(phi02));
00666 double m_m1vct=cos(mdcTrkm1->theta());
00667 double m_m1ptot=mdcKalTrkm1->p();
00668 double m_m1pxy=sqrt(mdcKalTrkm1->px()*mdcKalTrkm1->px()+mdcKalTrkm1->py()*mdcKalTrkm1->py());
00669
00670 if((*itTrkm1)->isEmcShowerValid()){
00671 emcTg2= emcTrkm1->energy();
00672 }
00673 if((*itTrkm1)->isMucTrackValid()){
00674 laypi2=mucTrkm1->numLayers();
00675 }
00676 double m_laypim1=laypi2;
00677
00678 EvtRecTrackIterator itTrkp2=evtRecTrkCol->begin() + ipip[1];
00679 RecMdcTrack* mdcTrkp2 = (*itTrkp2)->mdcTrack();
00680 RecMdcKalTrack *mdcKalTrkp2 = (*itTrkp2)->mdcKalTrack();
00681 RecEmcShower* emcTrkp2 = (*itTrkp2)->emcShower();
00682 RecMucTrack *mucTrkp2=(*itTrkp2)->mucTrack();
00683
00684 double phi03=mdcTrkp2->helix(1);
00685 double m_p2vx = mdcTrkp2->x();
00686 double m_p2vy = mdcTrkp2->y();
00687 double m_p2vz = mdcTrkp2->z()-xorigin.z();
00688 double m_p2vr = fabs((mdcTrkp2->x()-xorigin.x())*cos(phi03)+(mdcTrkp2->y()-xorigin.y())*sin(phi03));
00689 double m_p2vct=cos(mdcTrkp2->theta());
00690 double m_p2ptot=mdcKalTrkp2->p();
00691 double m_p2pxy=sqrt(mdcKalTrkp2->px()*mdcKalTrkp2->px()+mdcKalTrkp2->py()*mdcKalTrkp2->py());
00692
00693 if((*itTrkp2)->isEmcShowerValid()){
00694 emcTg3= emcTrkp2->energy();
00695 }
00696 if((*itTrkp2)->isMucTrackValid()){
00697 laypi3=mucTrkp2->numLayers();
00698 }
00699 double m_laypip2=laypi3;
00700
00701 EvtRecTrackIterator itTrkm2=evtRecTrkCol->begin() + ipim[1];
00702 RecMdcTrack* mdcTrkm2 = (*itTrkm2)->mdcTrack();
00703 RecMdcKalTrack *mdcKalTrkm2 = (*itTrkm2)->mdcKalTrack();
00704 RecEmcShower* emcTrkm2 = (*itTrkm2)->emcShower();
00705 RecMucTrack *mucTrkm2=(*itTrkm2)->mucTrack();
00706
00707 double phi04=mdcTrkm2->helix(1);
00708 double m_m2vx = mdcTrkm2->x();
00709 double m_m2vy = mdcTrkm2->y();
00710 double m_m2vz = mdcTrkm2->z()-xorigin.z();
00711 double m_m2vr = fabs((mdcTrkm2->x()-xorigin.x())*cos(phi04)+(mdcTrkm2->y()-xorigin.y())*sin(phi04));
00712 double m_m2vct=cos(mdcTrkm2->theta());
00713 double m_m2ptot=mdcKalTrkm2->p();
00714 double m_m2pxy=sqrt(mdcKalTrkm2->px()*mdcKalTrkm2->px()+mdcKalTrkm2->py()*mdcKalTrkm2->py());
00715
00716 if((*itTrkm2)->isEmcShowerValid()){
00717 emcTg4= emcTrkm2->energy();
00718 }
00719 if((*itTrkm2)->isMucTrackValid()){
00720 laypi4=mucTrkm2->numLayers();
00721 }
00722 double m_laypim2=laypi4;
00723
00724 double m_emcTp1 =emcTg1;
00725 double m_emcTm1 =emcTg2;
00726 double m_emcTp2 =emcTg3;
00727 double m_emcTm2 =emcTg4;
00728
00729 if(fabs(m_p1vz) >= m_vz1cut) return SUCCESS;
00730 if(m_p1vr >= m_vr1cut) return SUCCESS;
00731 if(fabs(m_p1vct)>=m_cthcut) return SUCCESS;
00732
00733 if(fabs(m_m1vz) >= m_vz1cut) return SUCCESS;
00734 if(m_m1vr >= m_vr1cut) return SUCCESS;
00735 if(fabs(m_m1vct)>=m_cthcut) return SUCCESS;
00736 Ncut5++;
00737
00738 HepLorentzVector p4muonp = ppip[1];
00739 HepLorentzVector p4muonm = ppim[1];
00740 HepLorentzVector p4uu = pTrecf;
00741
00742
00743 Hep3Vector p3jpsiUnit = (p4uu.vect()).unit();
00744 double jBeta = p4uu.beta();
00745
00746
00747
00748
00749
00750
00751
00752
00753 HepLorentzVector pTot;
00754 double minpi0=999.0;
00755 for(int i = 0; i < nGam - 1; i++){
00756 for(int j = i+1; j < nGam; j++) {
00757 HepLorentzVector p2g = pGam[i] + pGam[j];
00758 pTot = ppip[0] + ppim[0] + ppip[1] + ppim[1];
00759 pTot += p2g;
00760 if(fabs(p2g.m()-0.135)<minpi0){
00761 minpi0 = fabs(p2g.m()-0.135);
00762
00763 double m_m2gg = p2g.m();
00764
00765 HepLorentzVector prho0_no = ppi2_no2;
00766 HepLorentzVector prhop_no = ppip[1] + p2g;
00767 HepLorentzVector prhom_no = ppim[1] + p2g;
00768 HepLorentzVector prho0pi0 = ppi2_no2 + p2g;
00769 HepLorentzVector frho1pi0 = ppi2_no1 + p2g;
00770 HepLorentzVector frho2pi0 = ppi2_no3 + p2g;
00771 HepLorentzVector frho3pi0 = ppi2_no4 + p2g;
00772 HepLorentzVector prho0g1 = ppi2_no2 + pGam[i];
00773 HepLorentzVector prho0g2 = ppi2_no2 + pGam[j];
00774 HepLorentzVector frho1g1 = ppi2_no1 + pGam[i];
00775 HepLorentzVector frho1g2 = ppi2_no1 + pGam[j];
00776 HepLorentzVector frho2g1 = ppi2_no3 + pGam[i];
00777 HepLorentzVector frho2g2 = ppi2_no3 + pGam[j];
00778 HepLorentzVector frho3g1 = ppi2_no4 + pGam[i];
00779 HepLorentzVector frho3g2 = ppi2_no4 + pGam[j];
00780 HepLorentzVector p5pi_no = p4pi_no + p2g;
00781
00782
00783 double m_prho0_no = prho0_no.m();
00784 double m_prhop_no = prhop_no.m();
00785 double m_prhom_no = prhom_no.m();
00786 double m_prho0pi0 = prho0pi0.m();
00787 double m_frho1pi0 = frho1pi0.m();
00788 double m_frho2pi0 = frho2pi0.m();
00789 double m_frho3pi0 = frho3pi0.m();
00790 double m_prho0g1 = prho0g1.m();
00791 double m_prho0g2 = prho0g2.m();
00792 double m_frho1g1 = frho1g1.m();
00793 double m_frho1g2 = frho1g2.m();
00794 double m_frho2g1 = frho2g1.m();
00795 double m_frho2g2 = frho2g2.m();
00796 double m_frho3g1 = frho3g1.m();
00797 double m_frho3g2 = frho3g2.m();
00798 double m_p4pi_no = p4pi_no.m();
00799 double m_p5pi_no = p5pi_no.m();
00800 double m_mdcpx1=ppip[0].px();
00801 double m_mdcpy1=ppip[0].py();
00802 double m_mdcpz1=ppip[0].pz();
00803 double m_mdcpe1=ppip[0].e();
00804 double m_mdcpx2=ppim[0].px();
00805 double m_mdcpy2=ppim[0].py();
00806 double m_mdcpz2=ppim[0].pz();
00807 double m_mdcpe2=ppim[0].e();
00808 double m_mdcpx3=ppip[1].px();
00809 double m_mdcpy3=ppip[1].py();
00810 double m_mdcpz3=ppip[1].pz();
00811 double m_mdcpe3=ppip[1].e();
00812 double m_mdcpx4=ppim[1].px();
00813 double m_mdcpy4=ppim[1].py();
00814 double m_mdcpz4=ppim[1].pz();
00815 double m_mdcpe4=ppim[1].e();
00816 double m_mdcpxg1=pGam[i].px();
00817 double m_mdcpyg1=pGam[i].py();
00818 double m_mdcpzg1=pGam[i].pz();
00819 double m_mdcpeg1=pGam[i].e();
00820 double m_mdcpxg2=pGam[j].px();
00821 double m_mdcpyg2=pGam[j].py();
00822 double m_mdcpzg2=pGam[j].pz();
00823 double m_mdcpeg2=pGam[j].e();
00824 double m_etot = pTot.e();
00825 double m_mrecjp1=m_recjpsi1;
00826 double m_mrecjp2=m_recjpsi2;
00827 double m_mrecjp3=m_recjpsi3;
00828 double m_mrecjp4=m_recjpsi4;
00829
00830
00831 }
00832 }
00833 }
00834 Ncut6++;
00835
00836
00837
00838
00839
00840
00841 HepPoint3D vx(0., 0., 0.);
00842 HepSymMatrix Evx(3, 0);
00843 double bx = 1E+6;
00844 double by = 1E+6;
00845 double bz = 1E+6;
00846 Evx[0][0] = bx*bx;
00847 Evx[1][1] = by*by;
00848 Evx[2][2] = bz*bz;
00849
00850 VertexParameter vxpar;
00851 vxpar.setVx(vx);
00852 vxpar.setEvx(Evx);
00853
00854 VertexFit* vtxfit = VertexFit::instance();
00855 vtxfit->init();
00856
00857
00858
00859
00860 RecMdcKalTrack *pipTrk1 = (*(evtRecTrkCol->begin()+ipip[0]))->mdcKalTrack();
00861 RecMdcKalTrack *pimTrk1 = (*(evtRecTrkCol->begin()+ipim[0]))->mdcKalTrack();
00862 RecMdcKalTrack *pipTrk2 = (*(evtRecTrkCol->begin()+ipip[1]))->mdcKalTrack();
00863 RecMdcKalTrack *pimTrk2 = (*(evtRecTrkCol->begin()+ipim[1]))->mdcKalTrack();
00864
00865 WTrackParameter wvpipTrk1, wvpimTrk1,wvpipTrk2, wvpimTrk2;
00866 wvpipTrk1 = WTrackParameter(mpi, pipTrk1->getZHelix(), pipTrk1->getZError());
00867 wvpimTrk1 = WTrackParameter(mpi, pimTrk1->getZHelix(), pimTrk1->getZError());
00868 wvpipTrk2 = WTrackParameter(mpi, pipTrk2->getZHelix(), pipTrk2->getZError());
00869 wvpimTrk2 = WTrackParameter(mpi, pimTrk2->getZHelix(), pimTrk2->getZError());
00870
00871 vtxfit->AddTrack(0, wvpipTrk1);
00872 vtxfit->AddTrack(1, wvpimTrk1);
00873 vtxfit->AddVertex(0, vxpar,0, 1);
00874 if(!vtxfit->Fit(0)) return SUCCESS;
00875 vtxfit->Swim(0);
00876
00877 Ncut7++;
00878
00879 WTrackParameter wpip1 = vtxfit->wtrk(0);
00880 WTrackParameter wpim1 = vtxfit->wtrk(1);
00881
00882 KinematicFit * kmfit = KinematicFit::instance();
00883
00884
00885
00886
00887 int igbf1 = -1;
00888 int igbf2 = -1;
00889 HepLorentzVector pTgam1(0,0,0,0);
00890 HepLorentzVector pTgam2(0,0,0,0);
00891
00892 if(m_test4C==1) {
00893
00894 HepLorentzVector ecms(0.011*3.6862,0,0,3.6862);
00895
00896
00897
00898
00899 WTrackParameter wvkipTrk2, wvkimTrk2;
00900 wvkipTrk2 = WTrackParameter(mk, pipTrk2->getZHelixK(), pipTrk2->getZErrorK());
00901 wvkimTrk2 = WTrackParameter(mk, pimTrk2->getZHelixK(), pimTrk2->getZErrorK());
00902 double chisq = 9999.;
00903 int ig11 = -1;
00904 int ig21 = -1;
00905 double chikk=9999.;
00906 for(int i = 0; i < nGam-1; i++) {
00907 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
00908 for(int j = i+1; j < nGam; j++) {
00909 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[j]))->emcShower();
00910 kmfit->init();
00911 kmfit->setDynamicerror(1);
00912 kmfit->AddTrack(0, wpip1);
00913 kmfit->AddTrack(1, wpim1);
00914 kmfit->AddTrack(2, wvkipTrk2);
00915 kmfit->AddTrack(3, wvkimTrk2);
00916 kmfit->AddTrack(4, 0.0, g1Trk);
00917 kmfit->AddTrack(5, 0.0, g2Trk);
00918 kmfit->AddFourMomentum(0, ecms);
00919 bool oksq = kmfit->Fit();
00920 if(oksq&&kmfit->chisq()<chikk) {
00921 chikk = kmfit->chisq();
00922 }
00923 }
00924 }
00925 Ncut8++;
00926
00927
00928
00929
00930
00931 chisq = 9999.;
00932 int ig1 = -1;
00933 int ig2 = -1;
00934 for(int i = 0; i < nGam-1; i++) {
00935 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
00936 for(int j = i+1; j < nGam; j++) {
00937 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[j]))->emcShower();
00938 kmfit->init();
00939 kmfit->setDynamicerror(1);
00940 kmfit->AddTrack(0, wpip1);
00941 kmfit->AddTrack(1, wpim1);
00942 kmfit->AddTrack(2, wvpipTrk2);
00943 kmfit->AddTrack(3, wvpimTrk2);
00944 kmfit->AddTrack(4, 0.0, g1Trk);
00945 kmfit->AddTrack(5, 0.0, g2Trk);
00946 kmfit->AddFourMomentum(0, ecms);
00947 bool oksq = kmfit->Fit();
00948 if(oksq) {
00949 double chi2 = kmfit->chisq();
00950 if(chi2 < chisq) {
00951 chisq = chi2;
00952 ig1 = iGam[i];
00953 ig2 = iGam[j];
00954 igbf1 = iGam[i];
00955 igbf2 = iGam[j];
00956 pTgam1=pGam[i];
00957 pTgam2=pGam[j];
00958 }
00959 }
00960 }
00961 }
00962
00963
00964 if(chisq > 200) return SUCCESS;
00965 Ncut9++;
00966
00967
00968 Vint jGood;
00969 jGood.clear();
00970 jGood.push_back(ipip[0]);
00971 jGood.push_back(ipim[0]);
00972 jGood.push_back(ipip[1]);
00973 jGood.push_back(ipim[1]);
00974
00975
00976
00977 Vint jGgam;
00978 jGgam.clear();
00979 jGgam.push_back(igbf1);
00980 jGgam.push_back(igbf2);
00981
00982 double chi1_pp=9999.0;
00983
00984 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+ig1))->emcShower();
00985 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+ig2))->emcShower();
00986 kmfit->init();
00987 kmfit->AddTrack(0, wpip1);
00988 kmfit->AddTrack(1, wpim1);
00989 kmfit->AddTrack(2, wvpipTrk2);
00990 kmfit->AddTrack(3, wvpimTrk2);
00991 kmfit->AddTrack(4, 0.0, g1Trk);
00992 kmfit->AddTrack(5, 0.0, g2Trk);
00993 kmfit->AddFourMomentum(0, ecms);
00994 bool oksq = kmfit->Fit();
00995 if(oksq) {
00996 chi1_pp = kmfit->chisq();
00997 HepLorentzVector ppi0 = kmfit->pfit(4) + kmfit->pfit(5);
00998 HepLorentzVector prho0= kmfit->pfit(2) + kmfit->pfit(3);
00999 HepLorentzVector prhop= kmfit->pfit(2) + ppi0;
01000 HepLorentzVector prhom= kmfit->pfit(3) + ppi0;
01001 HepLorentzVector pjjj = prho0 + ppi0;
01002
01003 HepLorentzVector p2pi1=kmfit->pfit(0) + kmfit->pfit(1);
01004 HepLorentzVector f2pi1g1= p2pi1 + kmfit->pfit(4);
01005 HepLorentzVector f2pi1g2= p2pi1 + kmfit->pfit(5);
01006 HepLorentzVector f2pi1pi0=p2pi1 + ppi0;
01007
01008 HepLorentzVector t2pi2g1= prho0 + kmfit->pfit(4);
01009 HepLorentzVector t2pi2g2= prho0 + kmfit->pfit(5);
01010
01011 HepLorentzVector p2pi3=kmfit->pfit(0) + kmfit->pfit(3);
01012 HepLorentzVector f2pi3g1= p2pi3 + kmfit->pfit(4);
01013 HepLorentzVector f2pi3g2= p2pi3 + kmfit->pfit(5);
01014 HepLorentzVector f2pi3pi0=p2pi3 + ppi0;
01015
01016 HepLorentzVector p2pi4=kmfit->pfit(1) + kmfit->pfit(2);
01017 HepLorentzVector f2pi4g1= p2pi4 + kmfit->pfit(4);
01018 HepLorentzVector f2pi4g2= p2pi4 + kmfit->pfit(5);
01019 HepLorentzVector f2pi4pi0=p2pi4 + ppi0;
01020
01021 HepLorentzVector p4pi= p2pi1 + prho0;
01022 HepLorentzVector p4pig1= p4pi + kmfit->pfit(4);
01023 HepLorentzVector p4pig2= p4pi + kmfit->pfit(5);
01024 HepLorentzVector ppptot= p4pi + ppi0;
01025
01026
01027 HepLorentzVector be4cpi0= pTgam1 + pTgam2;
01028 HepLorentzVector be4c_ppi1 = ppip[0] + ppim[0];
01029 HepLorentzVector be4c_ppi2 = ppip[1] + ppim[1];
01030 HepLorentzVector be4cjp= be4cpi0 + be4c_ppi2;
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425 Ncut10++;
01426
01427 }
01428 }
01429
01430
01431
01432
01433
01434
01435
01436 setFilterPassed(true);
01437
01438 (*(evtRecTrkCol->begin()+ipip[0]))->setPartId(2);
01439 (*(evtRecTrkCol->begin()+ipim[0]))->setPartId(2);
01440 (*(evtRecTrkCol->begin()+ipip[1]))->setPartId(2);
01441 (*(evtRecTrkCol->begin()+ipim[1]))->setPartId(2);
01442
01443
01444 return StatusCode::SUCCESS;
01445 }
01446
01447
01448
01449 StatusCode Ppjrhopi::finalize() {
01450 cout<<"total number: "<<Ncut0<<endl;
01451 cout<<"nGood==4, nCharge==0: "<<Ncut1<<endl;
01452 cout<<"nGam>=2: "<<Ncut2<<endl;
01453 cout<<"Pass no Pid: "<<Ncut3<<endl;
01454 cout<<"ChangeID recfrom psp: "<<Ncut4<<endl;
01455 cout<<"vetex position: "<<Ncut5<<endl;
01456 cout<<"Mass from MDC: "<<Ncut6<<endl;
01457 cout<<"primary vetex fit: "<<Ncut7<<endl;
01458 cout<<"Pass 4C for ppkkp0: "<<Ncut8<<endl;
01459 cout<<"Pass 4C for 4pi+pi0: "<<Ncut9<<endl;
01460 cout<<"Pass 4C <200: "<<Ncut10<<endl;
01461 MsgStream log(msgSvc(), name());
01462 log << MSG::INFO << "in finalize()" << endmsg;
01463 return StatusCode::SUCCESS;
01464 }
01465