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 "RhopiAlg/Rhopi.h"
00037
00038
00039 #include "VertexFit/KalmanKinematicFit.h"
00040 #include "VertexFit/VertexFit.h"
00041 #include "VertexFit/Helix.h"
00042 #include "ParticleID/ParticleID.h"
00043
00044 #include <vector>
00045
00046
00047 const double mpi = 0.13957;
00048 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00049
00050 const double velc = 299.792458;
00051 typedef std::vector<int> Vint;
00052 typedef std::vector<HepLorentzVector> Vp4;
00053
00054 int Ncut0,Ncut1,Ncut2,Ncut3,Ncut4,Ncut5,Ncut6;
00055
00057
00058 Rhopi::Rhopi(const std::string& name, ISvcLocator* pSvcLocator) :
00059 Algorithm(name, pSvcLocator) {
00060
00061
00062 declareProperty("Vr0cut", m_vr0cut=1.0);
00063 declareProperty("Vz0cut", m_vz0cut=5.0);
00064 declareProperty("EnergyThreshold", m_energyThreshold=0.04);
00065 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
00066 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
00067 declareProperty("GammaAngleCut", m_gammaAngleCut=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 Rhopi::initialize(){
00076 MsgStream log(msgSvc(), name());
00077
00078 log << MSG::INFO << "in initialize()" << endmsg;
00079
00080 StatusCode status;
00081 NTuplePtr nt1(ntupleSvc(), "FILE1/vxyz");
00082 if ( nt1 ) m_tuple1 = nt1;
00083 else {
00084 m_tuple1 = ntupleSvc()->book ("FILE1/vxyz", CLID_ColumnWiseTuple, "ks N-Tuple example");
00085 if ( m_tuple1 ) {
00086 status = m_tuple1->addItem ("vx0", m_vx0);
00087 status = m_tuple1->addItem ("vy0", m_vy0);
00088 status = m_tuple1->addItem ("vz0", m_vz0);
00089 status = m_tuple1->addItem ("vr0", m_vr0);
00090 status = m_tuple1->addItem ("rvxy0", m_rvxy0);
00091 status = m_tuple1->addItem ("rvz0", m_rvz0);
00092 status = m_tuple1->addItem ("rvphi0", m_rvphi0);
00093 }
00094 else {
00095 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00096 return StatusCode::FAILURE;
00097 }
00098 }
00099
00100 NTuplePtr nt2(ntupleSvc(), "FILE1/photon");
00101 if ( nt2 ) m_tuple2 = nt2;
00102 else {
00103 m_tuple2 = ntupleSvc()->book ("FILE1/photon", CLID_ColumnWiseTuple, "ks N-Tuple example");
00104 if ( m_tuple2 ) {
00105 status = m_tuple2->addItem ("dthe", m_dthe);
00106 status = m_tuple2->addItem ("dphi", m_dphi);
00107 status = m_tuple2->addItem ("dang", m_dang);
00108 status = m_tuple2->addItem ("eraw", m_eraw);
00109 }
00110 else {
00111 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
00112 return StatusCode::FAILURE;
00113 }
00114 }
00115
00116
00117 NTuplePtr nt3(ntupleSvc(), "FILE1/etot");
00118 if ( nt3 ) m_tuple3 = nt3;
00119 else {
00120 m_tuple3 = ntupleSvc()->book ("FILE1/etot", CLID_ColumnWiseTuple, "ks N-Tuple example");
00121 if ( m_tuple3 ) {
00122 status = m_tuple3->addItem ("m2gg", m_m2gg);
00123 status = m_tuple3->addItem ("etot", m_etot);
00124 }
00125 else {
00126 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple3) << endmsg;
00127 return StatusCode::FAILURE;
00128 }
00129 }
00130 if(m_test4C==1) {
00131 NTuplePtr nt4(ntupleSvc(), "FILE1/fit4c");
00132 if ( nt4 ) m_tuple4 = nt4;
00133 else {
00134 m_tuple4 = ntupleSvc()->book ("FILE1/fit4c", CLID_ColumnWiseTuple, "ks N-Tuple example");
00135 if ( m_tuple4 ) {
00136 status = m_tuple4->addItem ("chi2", m_chi1);
00137 status = m_tuple4->addItem ("mpi0", m_mpi0);
00138 }
00139 else {
00140 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple4) << endmsg;
00141 return StatusCode::FAILURE;
00142 }
00143 }
00144 }
00145
00146
00147 if(m_test5C==1) {
00148 NTuplePtr nt5(ntupleSvc(), "FILE1/fit5c");
00149 if ( nt5 ) m_tuple5 = nt5;
00150 else {
00151 m_tuple5 = ntupleSvc()->book ("FILE1/fit5c", CLID_ColumnWiseTuple, "ks N-Tuple example");
00152 if ( m_tuple5 ) {
00153 status = m_tuple5->addItem ("chi2", m_chi2);
00154 status = m_tuple5->addItem ("mrh0", m_mrh0);
00155 status = m_tuple5->addItem ("mrhp", m_mrhp);
00156 status = m_tuple5->addItem ("mrhm", m_mrhm);
00157 }
00158 else {
00159 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple5) << endmsg;
00160 return StatusCode::FAILURE;
00161 }
00162 }
00163
00164 NTuplePtr nt6(ntupleSvc(), "FILE1/geff");
00165 if ( nt6 ) m_tuple6 = nt6;
00166 else {
00167 m_tuple6 = ntupleSvc()->book ("FILE1/geff", CLID_ColumnWiseTuple, "ks N-Tuple example");
00168 if ( m_tuple6 ) {
00169 status = m_tuple6->addItem ("fcos", m_fcos);
00170 status = m_tuple6->addItem ("elow", m_elow);
00171 }
00172 else {
00173 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple6) << endmsg;
00174 return StatusCode::FAILURE;
00175 }
00176 }
00177 }
00178
00179 if(m_checkDedx == 1) {
00180 NTuplePtr nt7(ntupleSvc(), "FILE1/dedx");
00181 if ( nt7 ) m_tuple7 = nt7;
00182 else {
00183 m_tuple7 = ntupleSvc()->book ("FILE1/dedx", CLID_ColumnWiseTuple, "ks N-Tuple example");
00184 if ( m_tuple7 ) {
00185 status = m_tuple7->addItem ("ptrk", m_ptrk);
00186 status = m_tuple7->addItem ("chie", m_chie);
00187 status = m_tuple7->addItem ("chimu", m_chimu);
00188 status = m_tuple7->addItem ("chipi", m_chipi);
00189 status = m_tuple7->addItem ("chik", m_chik);
00190 status = m_tuple7->addItem ("chip", m_chip);
00191 status = m_tuple7->addItem ("probPH", m_probPH);
00192 status = m_tuple7->addItem ("normPH", m_normPH);
00193 status = m_tuple7->addItem ("ghit", m_ghit);
00194 status = m_tuple7->addItem ("thit", m_thit);
00195 }
00196 else {
00197 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple7) << endmsg;
00198 return StatusCode::FAILURE;
00199 }
00200 }
00201 }
00202
00203 if(m_checkTof == 1) {
00204 NTuplePtr nt8(ntupleSvc(), "FILE1/tofe");
00205 if ( nt8 ) m_tuple8 = nt8;
00206 else {
00207 m_tuple8 = ntupleSvc()->book ("FILE1/tofe",CLID_ColumnWiseTuple, "ks N-Tuple example");
00208 if ( m_tuple8 ) {
00209 status = m_tuple8->addItem ("ptrk", m_ptot_etof);
00210 status = m_tuple8->addItem ("cntr", m_cntr_etof);
00211 status = m_tuple8->addItem ("ph", m_ph_etof);
00212 status = m_tuple8->addItem ("rhit", m_rhit_etof);
00213 status = m_tuple8->addItem ("qual", m_qual_etof);
00214 status = m_tuple8->addItem ("te", m_te_etof);
00215 status = m_tuple8->addItem ("tmu", m_tmu_etof);
00216 status = m_tuple8->addItem ("tpi", m_tpi_etof);
00217 status = m_tuple8->addItem ("tk", m_tk_etof);
00218 status = m_tuple8->addItem ("tp", m_tp_etof);
00219 }
00220 else {
00221 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple8) << endmsg;
00222 return StatusCode::FAILURE;
00223 }
00224 }
00225 }
00226
00227
00228
00229 if(m_checkTof == 1) {
00230 NTuplePtr nt9(ntupleSvc(), "FILE1/tof1");
00231 if ( nt9 ) m_tuple9 = nt9;
00232 else {
00233 m_tuple9 = ntupleSvc()->book ("FILE1/tof1", CLID_ColumnWiseTuple, "ks N-Tuple example");
00234 if ( m_tuple9 ) {
00235 status = m_tuple9->addItem ("ptrk", m_ptot_btof1);
00236 status = m_tuple9->addItem ("cntr", m_cntr_btof1);
00237 status = m_tuple9->addItem ("ph", m_ph_btof1);
00238 status = m_tuple9->addItem ("zhit", m_zhit_btof1);
00239 status = m_tuple9->addItem ("qual", m_qual_btof1);
00240 status = m_tuple9->addItem ("te", m_te_btof1);
00241 status = m_tuple9->addItem ("tmu", m_tmu_btof1);
00242 status = m_tuple9->addItem ("tpi", m_tpi_btof1);
00243 status = m_tuple9->addItem ("tk", m_tk_btof1);
00244 status = m_tuple9->addItem ("tp", m_tp_btof1);
00245 }
00246 else {
00247 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple9) << endmsg;
00248 return StatusCode::FAILURE;
00249 }
00250 }
00251 }
00252
00253
00254 if(m_checkTof == 1) {
00255 NTuplePtr nt10(ntupleSvc(), "FILE1/tof2");
00256 if ( nt10 ) m_tuple10 = nt10;
00257 else {
00258 m_tuple10 = ntupleSvc()->book ("FILE1/tof2", CLID_ColumnWiseTuple, "ks N-Tuple example");
00259 if ( m_tuple10 ) {
00260 status = m_tuple10->addItem ("ptrk", m_ptot_btof2);
00261 status = m_tuple10->addItem ("cntr", m_cntr_btof2);
00262 status = m_tuple10->addItem ("ph", m_ph_btof2);
00263 status = m_tuple10->addItem ("zhit", m_zhit_btof2);
00264 status = m_tuple10->addItem ("qual", m_qual_btof2);
00265 status = m_tuple10->addItem ("te", m_te_btof2);
00266 status = m_tuple10->addItem ("tmu", m_tmu_btof2);
00267 status = m_tuple10->addItem ("tpi", m_tpi_btof2);
00268 status = m_tuple10->addItem ("tk", m_tk_btof2);
00269 status = m_tuple10->addItem ("tp", m_tp_btof2);
00270 }
00271 else {
00272 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple10) << endmsg;
00273 return StatusCode::FAILURE;
00274 }
00275 }
00276 }
00277
00278
00279 NTuplePtr nt11(ntupleSvc(), "FILE1/pid");
00280 if ( nt11 ) m_tuple11 = nt11;
00281 else {
00282 m_tuple11 = ntupleSvc()->book ("FILE1/pid", CLID_ColumnWiseTuple, "ks N-Tuple example");
00283 if ( m_tuple11 ) {
00284 status = m_tuple11->addItem ("ptrk", m_ptrk_pid);
00285 status = m_tuple11->addItem ("cost", m_cost_pid);
00286 status = m_tuple11->addItem ("dedx", m_dedx_pid);
00287 status = m_tuple11->addItem ("tof1", m_tof1_pid);
00288 status = m_tuple11->addItem ("tof2", m_tof2_pid);
00289 status = m_tuple11->addItem ("prob", m_prob_pid);
00290 }
00291 else {
00292 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple11) << endmsg;
00293 return StatusCode::FAILURE;
00294 }
00295 }
00296
00297
00298
00299
00300
00301
00302 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00303 return StatusCode::SUCCESS;
00304
00305 }
00306
00307
00308 StatusCode Rhopi::execute() {
00309
00310 std::cout << "execute()" << std::endl;
00311
00312 MsgStream log(msgSvc(), name());
00313 log << MSG::INFO << "in execute()" << endreq;
00314
00315 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00316 int runNo=eventHeader->runNumber();
00317 int event=eventHeader->eventNumber();
00318 log << MSG::DEBUG <<"run, evtnum = "
00319 << runNo << " , "
00320 << event <<endreq;
00321 cout<<"event "<<event<<endl;
00322 Ncut0++;
00323
00324 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00325
00326 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00327 << evtRecEvent->totalCharged() << " , "
00328 << evtRecEvent->totalNeutral() << " , "
00329 << evtRecEvent->totalTracks() <<endreq;
00330
00331 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00332
00333
00334
00335
00336 Vint iGood, ipip, ipim;
00337 iGood.clear();
00338 ipip.clear();
00339 ipim.clear();
00340 Vp4 ppip, ppim;
00341 ppip.clear();
00342 ppim.clear();
00343
00344 int nCharge = 0;
00345
00346 Hep3Vector xorigin(0,0,0);
00347
00348
00349 IVertexDbSvc* vtxsvc;
00350 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00351 if(vtxsvc->isVertexValid()){
00352 double* dbv = vtxsvc->PrimaryVertex();
00353 double* vv = vtxsvc->SigmaPrimaryVertex();
00354
00355
00356 xorigin.setX(dbv[0]);
00357 xorigin.setY(dbv[1]);
00358 xorigin.setZ(dbv[2]);
00359 }
00360
00361 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00362 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00363 if(!(*itTrk)->isMdcTrackValid()) continue;
00364 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
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=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
00373 m_vx0 = x0;
00374 m_vy0 = y0;
00375 m_vz0 = z0;
00376 m_vr0 = Rxy;
00377
00378 HepVector a = mdcTrk->helix();
00379 HepSymMatrix Ea = mdcTrk->err();
00380 HepPoint3D point0(0.,0.,0.);
00381 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
00382 VFHelix helixip(point0,a,Ea);
00383 helixip.pivot(IP);
00384 HepVector vecipa = helixip.a();
00385 double Rvxy0=fabs(vecipa[0]);
00386 double Rvz0=vecipa[3];
00387 double Rvphi0=vecipa[1];
00388 m_rvxy0=Rvxy0;
00389 m_rvz0=Rvz0;
00390 m_rvphi0=Rvphi0;
00391
00392 m_tuple1->write();
00393
00394
00395
00396 if(fabs(Rvz0) >= 10.0) continue;
00397 if(fabs(Rvxy0) >= 1.0) continue;
00398
00399 iGood.push_back(i);
00400 nCharge += mdcTrk->charge();
00401 }
00402
00403
00404
00405
00406 int nGood = iGood.size();
00407 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
00408 if((nGood != 2)||(nCharge!=0)){
00409 return StatusCode::SUCCESS;
00410 }
00411 Ncut1++;
00412
00413 Vint iGam;
00414 iGam.clear();
00415 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
00416 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00417 if(!(*itTrk)->isEmcShowerValid()) continue;
00418 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00419 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
00420
00421 double dthe = 200.;
00422 double dphi = 200.;
00423 double dang = 200.;
00424 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
00425 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
00426 if(!(*jtTrk)->isExtTrackValid()) continue;
00427 RecExtTrack *extTrk = (*jtTrk)->extTrack();
00428 if(extTrk->emcVolumeNumber() == -1) continue;
00429 Hep3Vector extpos = extTrk->emcPosition();
00430
00431 double angd = extpos.angle(emcpos);
00432 double thed = extpos.theta() - emcpos.theta();
00433 double phid = extpos.deltaPhi(emcpos);
00434 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00435 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00436 if(angd < dang){
00437 dang = angd;
00438 dthe = thed;
00439 dphi = phid;
00440 }
00441 }
00442 if(dang>=200) continue;
00443 double eraw = emcTrk->energy();
00444 dthe = dthe * 180 / (CLHEP::pi);
00445 dphi = dphi * 180 / (CLHEP::pi);
00446 dang = dang * 180 / (CLHEP::pi);
00447 m_dthe = dthe;
00448 m_dphi = dphi;
00449 m_dang = dang;
00450 m_eraw = eraw;
00451 m_tuple2->write();
00452 if(eraw < m_energyThreshold) continue;
00453
00454 if(fabs(dang) < m_gammaAngleCut) continue;
00455
00456
00457
00458 iGam.push_back(i);
00459 }
00460
00461
00462
00463
00464 int nGam = iGam.size();
00465
00466 log << MSG::DEBUG << "num Good Photon " << nGam << " , " <<evtRecEvent->totalNeutral()<<endreq;
00467 if(nGam<2){
00468 return StatusCode::SUCCESS;
00469 }
00470 Ncut2++;
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480 if(m_checkDedx == 1) {
00481 for(int i = 0; i < nGood; i++) {
00482 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00483 if(!(*itTrk)->isMdcTrackValid()) continue;
00484 if(!(*itTrk)->isMdcDedxValid())continue;
00485 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00486 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
00487 m_ptrk = mdcTrk->p();
00488
00489 m_chie = dedxTrk->chiE();
00490 m_chimu = dedxTrk->chiMu();
00491 m_chipi = dedxTrk->chiPi();
00492 m_chik = dedxTrk->chiK();
00493 m_chip = dedxTrk->chiP();
00494 m_ghit = dedxTrk->numGoodHits();
00495 m_thit = dedxTrk->numTotalHits();
00496 m_probPH = dedxTrk->probPH();
00497 m_normPH = dedxTrk->normPH();
00498 m_tuple7->write();
00499 }
00500 }
00501
00502
00503
00504
00505
00506
00507 if(m_checkTof == 1) {
00508 for(int i = 0; i < nGood; i++) {
00509 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00510 if(!(*itTrk)->isMdcTrackValid()) continue;
00511 if(!(*itTrk)->isTofTrackValid()) continue;
00512
00513 RecMdcTrack * mdcTrk = (*itTrk)->mdcTrack();
00514 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
00515
00516 double ptrk = mdcTrk->p();
00517
00518 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
00519 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
00520 TofHitStatus *status = new TofHitStatus;
00521 status->setStatus((*iter_tof)->status());
00522 if(!(status->is_barrel())){
00523 if( !(status->is_counter()) ) continue;
00524 if( status->layer()!=0 ) continue;
00525 double path=(*iter_tof)->path();
00526 double tof = (*iter_tof)->tof();
00527 double ph = (*iter_tof)->ph();
00528 double rhit = (*iter_tof)->zrhit();
00529 double qual = 0.0 + (*iter_tof)->quality();
00530 double cntr = 0.0 + (*iter_tof)->tofID();
00531 double texp[5];
00532 for(int j = 0; j < 5; j++) {
00533 double gb = ptrk/xmass[j];
00534 double beta = gb/sqrt(1+gb*gb);
00535 texp[j] = 10 * path /beta/velc;
00536 }
00537 m_cntr_etof = cntr;
00538 m_ptot_etof = ptrk;
00539 m_ph_etof = ph;
00540 m_rhit_etof = rhit;
00541 m_qual_etof = qual;
00542 m_te_etof = tof - texp[0];
00543 m_tmu_etof = tof - texp[1];
00544 m_tpi_etof = tof - texp[2];
00545 m_tk_etof = tof - texp[3];
00546 m_tp_etof = tof - texp[4];
00547 m_tuple8->write();
00548 }
00549 else {
00550 if( !(status->is_counter()) ) continue;
00551 if(status->layer()==1){
00552 double path=(*iter_tof)->path();
00553 double tof = (*iter_tof)->tof();
00554 double ph = (*iter_tof)->ph();
00555 double rhit = (*iter_tof)->zrhit();
00556 double qual = 0.0 + (*iter_tof)->quality();
00557 double cntr = 0.0 + (*iter_tof)->tofID();
00558 double texp[5];
00559 for(int j = 0; j < 5; j++) {
00560 double gb = ptrk/xmass[j];
00561 double beta = gb/sqrt(1+gb*gb);
00562 texp[j] = 10 * path /beta/velc;
00563 }
00564
00565 m_cntr_btof1 = cntr;
00566 m_ptot_btof1 = ptrk;
00567 m_ph_btof1 = ph;
00568 m_zhit_btof1 = rhit;
00569 m_qual_btof1 = qual;
00570 m_te_btof1 = tof - texp[0];
00571 m_tmu_btof1 = tof - texp[1];
00572 m_tpi_btof1 = tof - texp[2];
00573 m_tk_btof1 = tof - texp[3];
00574 m_tp_btof1 = tof - texp[4];
00575 m_tuple9->write();
00576 }
00577
00578 if(status->layer()==2){
00579 double path=(*iter_tof)->path();
00580 double tof = (*iter_tof)->tof();
00581 double ph = (*iter_tof)->ph();
00582 double rhit = (*iter_tof)->zrhit();
00583 double qual = 0.0 + (*iter_tof)->quality();
00584 double cntr = 0.0 + (*iter_tof)->tofID();
00585 double texp[5];
00586 for(int j = 0; j < 5; j++) {
00587 double gb = ptrk/xmass[j];
00588 double beta = gb/sqrt(1+gb*gb);
00589 texp[j] = 10 * path /beta/velc;
00590 }
00591
00592 m_cntr_btof2 = cntr;
00593 m_ptot_btof2 = ptrk;
00594 m_ph_btof2 = ph;
00595 m_zhit_btof2 = rhit;
00596 m_qual_btof2 = qual;
00597 m_te_btof2 = tof - texp[0];
00598 m_tmu_btof2 = tof - texp[1];
00599 m_tpi_btof2 = tof - texp[2];
00600 m_tk_btof2 = tof - texp[3];
00601 m_tp_btof2 = tof - texp[4];
00602 m_tuple10->write();
00603 }
00604 }
00605
00606 delete status;
00607 }
00608 }
00609 }
00610
00611
00612
00613
00614
00615
00616 Vp4 pGam;
00617 pGam.clear();
00618 for(int i = 0; i < nGam; i++) {
00619 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
00620 RecEmcShower* emcTrk = (*itTrk)->emcShower();
00621 double eraw = emcTrk->energy();
00622 double phi = emcTrk->phi();
00623 double the = emcTrk->theta();
00624 HepLorentzVector ptrk;
00625 ptrk.setPx(eraw*sin(the)*cos(phi));
00626 ptrk.setPy(eraw*sin(the)*sin(phi));
00627 ptrk.setPz(eraw*cos(the));
00628 ptrk.setE(eraw);
00629
00630
00631
00632 pGam.push_back(ptrk);
00633 }
00634 cout<<"before pid"<<endl;
00635
00636
00637
00638 ParticleID *pid = ParticleID::instance();
00639 for(int i = 0; i < nGood; i++) {
00640 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00641
00642 pid->init();
00643 pid->setMethod(pid->methodProbability());
00644
00645
00646 pid->setChiMinCut(4);
00647 pid->setRecTrack(*itTrk);
00648 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2() | pid->useTofE());
00649 pid->identify(pid->onlyPion() | pid->onlyKaon());
00650
00651
00652 pid->calculate();
00653 if(!(pid->IsPidInfoValid())) continue;
00654 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00655 m_ptrk_pid = mdcTrk->p();
00656 m_cost_pid = cos(mdcTrk->theta());
00657 m_dedx_pid = pid->chiDedx(2);
00658 m_tof1_pid = pid->chiTof1(2);
00659 m_tof2_pid = pid->chiTof2(2);
00660 m_prob_pid = pid->probPion();
00661 m_tuple11->write();
00662
00663
00664 if(pid->probPion() < 0.001) continue;
00665
00666
00667 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
00668 RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);
00669
00670 if(mdcKalTrk->charge() >0) {
00671 ipip.push_back(iGood[i]);
00672 HepLorentzVector ptrk;
00673 ptrk.setPx(mdcKalTrk->px());
00674 ptrk.setPy(mdcKalTrk->py());
00675 ptrk.setPz(mdcKalTrk->pz());
00676 double p3 = ptrk.mag();
00677 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00678
00679
00680
00681 ppip.push_back(ptrk);
00682 } else {
00683 ipim.push_back(iGood[i]);
00684 HepLorentzVector ptrk;
00685 ptrk.setPx(mdcKalTrk->px());
00686 ptrk.setPy(mdcKalTrk->py());
00687 ptrk.setPz(mdcKalTrk->pz());
00688 double p3 = ptrk.mag();
00689 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00690
00691
00692
00693 ppim.push_back(ptrk);
00694 }
00695 }
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723 int npip = ipip.size();
00724 int npim = ipim.size();
00725 if(npip*npim != 1) return SUCCESS;
00726
00727 Ncut3++;
00728
00729
00730
00731
00732
00733
00734 HepLorentzVector pTot;
00735 for(int i = 0; i < nGam - 1; i++){
00736 for(int j = i+1; j < nGam; j++) {
00737 HepLorentzVector p2g = pGam[i] + pGam[j];
00738 pTot = ppip[0] + ppim[0];
00739 pTot += p2g;
00740 m_m2gg = p2g.m();
00741 m_etot = pTot.e();
00742 m_tuple3 -> write();
00743
00744 }
00745 }
00746
00747
00748 RecMdcKalTrack *pipTrk = (*(evtRecTrkCol->begin()+ipip[0]))->mdcKalTrack();
00749 RecMdcKalTrack *pimTrk = (*(evtRecTrkCol->begin()+ipim[0]))->mdcKalTrack();
00750
00751 WTrackParameter wvpipTrk, wvpimTrk;
00752 wvpipTrk = WTrackParameter(mpi, pipTrk->getZHelix(), pipTrk->getZError());
00753 wvpimTrk = WTrackParameter(mpi, pimTrk->getZHelix(), pimTrk->getZError());
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765 HepPoint3D vx(0., 0., 0.);
00766 HepSymMatrix Evx(3, 0);
00767 double bx = 1E+6;
00768 double by = 1E+6;
00769 double bz = 1E+6;
00770 Evx[0][0] = bx*bx;
00771 Evx[1][1] = by*by;
00772 Evx[2][2] = bz*bz;
00773
00774 VertexParameter vxpar;
00775 vxpar.setVx(vx);
00776 vxpar.setEvx(Evx);
00777
00778 VertexFit* vtxfit = VertexFit::instance();
00779 vtxfit->init();
00780 vtxfit->AddTrack(0, wvpipTrk);
00781 vtxfit->AddTrack(1, wvpimTrk);
00782 vtxfit->AddVertex(0, vxpar,0, 1);
00783 if(!vtxfit->Fit(0)) return SUCCESS;
00784 vtxfit->Swim(0);
00785
00786 WTrackParameter wpip = vtxfit->wtrk(0);
00787 WTrackParameter wpim = vtxfit->wtrk(1);
00788
00789
00790 KalmanKinematicFit * kmfit = KalmanKinematicFit::instance();
00791
00792
00793
00794
00795 cout<<"before 4c"<<endl;
00796 if(m_test4C==1) {
00797
00798 HepLorentzVector ecms(0.034,0,0,3.097);
00799
00800 double chisq = 9999.;
00801 int ig1 = -1;
00802 int ig2 = -1;
00803 for(int i = 0; i < nGam-1; i++) {
00804 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
00805 for(int j = i+1; j < nGam; j++) {
00806 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[j]))->emcShower();
00807 kmfit->init();
00808 kmfit->AddTrack(0, wpip);
00809 kmfit->AddTrack(1, wpim);
00810 kmfit->AddTrack(2, 0.0, g1Trk);
00811 kmfit->AddTrack(3, 0.0, g2Trk);
00812 kmfit->AddFourMomentum(0, ecms);
00813 bool oksq = kmfit->Fit();
00814 if(oksq) {
00815 double chi2 = kmfit->chisq();
00816 if(chi2 < chisq) {
00817 chisq = chi2;
00818 ig1 = iGam[i];
00819 ig2 = iGam[j];
00820 }
00821 }
00822 }
00823 }
00824
00825 if(chisq < 200) {
00826
00827 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+ig1))->emcShower();
00828 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+ig2))->emcShower();
00829 kmfit->init();
00830 kmfit->AddTrack(0, wpip);
00831 kmfit->AddTrack(1, wpim);
00832 kmfit->AddTrack(2, 0.0, g1Trk);
00833 kmfit->AddTrack(3, 0.0, g2Trk);
00834 kmfit->AddFourMomentum(0, ecms);
00835 bool oksq = kmfit->Fit();
00836 if(oksq) {
00837 HepLorentzVector ppi0 = kmfit->pfit(2) + kmfit->pfit(3);
00838 m_mpi0 = ppi0.m();
00839 m_chi1 = kmfit->chisq();
00840 m_tuple4->write();
00841 Ncut4++;
00842 }
00843 }
00844 }
00845
00846
00847
00848
00849
00850
00851 if(m_test5C==1) {
00852
00853 HepLorentzVector ecms(0.034,0,0,3.097);
00854 double chisq = 9999.;
00855 int ig1 = -1;
00856 int ig2 = -1;
00857 for(int i = 0; i < nGam-1; i++) {
00858 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
00859 for(int j = i+1; j < nGam; j++) {
00860 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[j]))->emcShower();
00861 kmfit->init();
00862 kmfit->AddTrack(0, wpip);
00863 kmfit->AddTrack(1, wpim);
00864 kmfit->AddTrack(2, 0.0, g1Trk);
00865 kmfit->AddTrack(3, 0.0, g2Trk);
00866 kmfit->AddResonance(0, 0.135, 2, 3);
00867 kmfit->AddFourMomentum(1, ecms);
00868 if(!kmfit->Fit(0)) continue;
00869 if(!kmfit->Fit(1)) continue;
00870 bool oksq = kmfit->Fit();
00871 if(oksq) {
00872 double chi2 = kmfit->chisq();
00873 if(chi2 < chisq) {
00874 chisq = chi2;
00875 ig1 = iGam[i];
00876 ig2 = iGam[j];
00877 }
00878 }
00879 }
00880 }
00881
00882
00883 log << MSG::INFO << " chisq = " << chisq <<endreq;
00884
00885 if(chisq < 200) {
00886 RecEmcShower* g1Trk = (*(evtRecTrkCol->begin()+ig1))->emcShower();
00887 RecEmcShower* g2Trk = (*(evtRecTrkCol->begin()+ig2))->emcShower();
00888
00889 kmfit->init();
00890 kmfit->AddTrack(0, wpip);
00891 kmfit->AddTrack(1, wpim);
00892 kmfit->AddTrack(2, 0.0, g1Trk);
00893 kmfit->AddTrack(3, 0.0, g2Trk);
00894 kmfit->AddResonance(0, 0.135, 2, 3);
00895 kmfit->AddFourMomentum(1, ecms);
00896 bool oksq = kmfit->Fit();
00897 if(oksq){
00898 HepLorentzVector ppi0 = kmfit->pfit(2) + kmfit->pfit(3);
00899 HepLorentzVector prho0 = kmfit->pfit(0) + kmfit->pfit(1);
00900 HepLorentzVector prhop = ppi0 + kmfit->pfit(0);
00901 HepLorentzVector prhom = ppi0 + kmfit->pfit(1);
00902
00903 m_chi2 = kmfit->chisq();
00904 m_mrh0 = prho0.m();
00905 m_mrhp = prhop.m();
00906 m_mrhm = prhom.m();
00907 double eg1 = (kmfit->pfit(2)).e();
00908 double eg2 = (kmfit->pfit(3)).e();
00909 double fcos = abs(eg1-eg2)/ppi0.rho();
00910 m_tuple5->write();
00911 Ncut5++;
00912
00913
00914
00915
00916 if(fabs(prho0.m()-0.770)<0.150) {
00917 if(fabs(fcos)<0.99) {
00918 m_fcos = (eg1-eg2)/ppi0.rho();
00919 m_elow = (eg1 < eg2) ? eg1 : eg2;
00920 m_tuple6->write();
00921 Ncut6++;
00922 }
00923 }
00924 }
00925 }
00926 }
00927 return StatusCode::SUCCESS;
00928 }
00929
00930
00931
00932 StatusCode Rhopi::finalize() {
00933 cout<<"total number: "<<Ncut0<<endl;
00934 cout<<"nGood==2, nCharge==0: "<<Ncut1<<endl;
00935 cout<<"nGam>=2: "<<Ncut2<<endl;
00936 cout<<"Pass Pid: "<<Ncut3<<endl;
00937 cout<<"Pass 4C: "<<Ncut4<<endl;
00938 cout<<"Pass 5C: "<<Ncut5<<endl;
00939 cout<<"J/psi->rho0 pi0: "<<Ncut6<<endl;
00940 MsgStream log(msgSvc(), name());
00941 log << MSG::INFO << "in finalize()" << endmsg;
00942 return StatusCode::SUCCESS;
00943 }
00944