BOSS6.5.5: DQAJpsi2PPbarAlg Class Reference

00207                                      {
00208   
00209   MsgStream log(msgSvc(), name());
00210   log << MSG::INFO << "in execute()" << endreq;
00211 
00212   // DQA
00213   // Add the line below at the beginning of execute()
00214   //
00215   setFilterPassed(false);
00216 
00217   counter[0]++;
00218   log << MSG::INFO <<"counter[0]="<<counter[0]<<endreq;
00219 
00220   SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00221   m_runNo=eventHeader->runNumber();
00222   m_event=eventHeader->eventNumber();
00223   log << MSG::DEBUG <<"run, evtnum = "
00224       << m_runNo << " , "
00225       << m_event <<endreq;
00226 
00227 
00228   SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00229   log << MSG::DEBUG <<"ncharg, nneu, tottks = " 
00230       << evtRecEvent->totalCharged() << " , "
00231       << evtRecEvent->totalNeutral() << " , "
00232       << evtRecEvent->totalTracks() <<endreq;
00233 
00234   SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00235 
00236   if(evtRecEvent->totalNeutral()>100) {
00237       return StatusCode::SUCCESS;
00238   }
00239 
00241   // read beam(vertex) position from the database    
00243   Hep3Vector xorigin(0,0,0);
00244  
00245   IVertexDbSvc*  vtxsvc;
00246   Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00247   if(vtxsvc->isVertexValid()){ 
00248   double* dbv = vtxsvc->PrimaryVertex();
00249   double*  vv = vtxsvc->SigmaPrimaryVertex();
00250     xorigin.setX(dbv[0]);
00251     xorigin.setY(dbv[1]);
00252     xorigin.setZ(dbv[2]);
00253   }
00254 
00256   // select good charged tracks   
00258   Vint iGood, iptrk, imtrk;
00259   iGood.clear();
00260   iptrk.clear();
00261   imtrk.clear();
00262   int nCharge = 0;
00263   for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00264     EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00265     if(!(*itTrk)->isMdcTrackValid()) continue;
00266     if (!(*itTrk)->isMdcKalTrackValid()) continue;
00267     RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00268 
00269     double pch =mdcTrk->p();
00270     double x0  =mdcTrk->x(); 
00271     double y0  =mdcTrk->y(); 
00272     double z0  =mdcTrk->z(); 
00273     double phi0=mdcTrk->helix(1);
00274     double xv=xorigin.x();
00275     double yv=xorigin.y();
00276     double Rxy=fabs((x0-xv)*cos(phi0)+(y0-yv)*sin(phi0));
00277 
00278     if(fabs(z0) >= m_vz1cut) continue;
00279     if(Rxy >= m_vr1cut) continue;
00280 //    if(fabs(m_Vct)>=m_cthcut) continue;
00281 
00282 //    iGood.push_back((*itTrk)->trackId());
00283     iGood.push_back(i);
00284     nCharge += mdcTrk->charge();
00285     if (mdcTrk->charge() > 0) {
00286         iptrk.push_back(i);
00287     } else {
00288         imtrk.push_back(i);
00289     }
00290   }
00291   
00292   //
00293   // Finish Good Charged Track Selection
00294   //
00295   int m_ngch = iGood.size();
00296   log << MSG::DEBUG << "ngood, totcharge = " << m_ngch << " , " << nCharge << endreq;
00297   if((m_ngch != 2)||(nCharge!=0)){
00298       return StatusCode::SUCCESS;
00299   }
00300   counter[1]++;
00301   log << MSG::INFO <<"counter[1]="<<counter[1]<<endreq;
00302 
00304   // read information for good charged tracks   
00306   Vp4 p_pp, p_pm; 
00307   p_pp.clear() ; 
00308   p_pm.clear() ; 
00309 
00310   int ii; 
00311   for(int i = 0; i < m_ngch; i++ ){
00312 
00313     // MDC information
00314     log << MSG::INFO <<"MDC information"<<endreq;
00315 
00316     EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + iGood[i];
00317     if(!(*itTrk)->isMdcTrackValid()) continue;
00318     RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00319     if(!(*itTrk)->isMdcKalTrackValid()) continue;
00320     RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
00321     
00322     if (mdcTrk->charge() > 0 ) {
00323       ii = 0 ; 
00324     }else{
00325       ii = 1 ; 
00326     }
00327 
00328     m_charge[ii] = mdcTrk->charge(); 
00329 
00330     double x0=mdcTrk->x();
00331     double y0=mdcTrk->y();
00332     double z0=mdcTrk->z();
00333     double phi0=mdcTrk->helix(1);
00334     double xv=xorigin.x();
00335     double yv=xorigin.y();
00336     double zv=xorigin.z();
00337     double rv=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
00338 
00339     m_vx0[ii]    = x0-xv ;
00340     m_vy0[ii]    = y0-yv ;
00341     m_vz0[ii]    = z0-zv ;
00342     m_vr0[ii]    = rv ;
00343     
00344     m_px[ii]   = mdcKalTrk->px();
00345     m_py[ii]   = mdcKalTrk->py();
00346     m_pz[ii]   = mdcKalTrk->pz();
00347     m_p[ii]    = mdcKalTrk->p();
00348     m_cos[ii]  = mdcKalTrk->pz()/mdcKalTrk->p();
00349 
00350     double e = sqrt(mdcKalTrk->p()*mdcKalTrk->p()+xmass[4]*xmass[4]);
00351     HepLorentzVector ptemp(m_px[ii],m_py[ii],m_pz[ii],e);
00352     if ( ii ==0 ) {
00353       p_pp.push_back(ptemp); 
00354     }else {
00355       p_pm.push_back(ptemp); 
00356     }
00357 
00358     // DEDX information  
00359     log << MSG::INFO <<"DEDX information"<<endreq;
00360     if ((*itTrk)->isMdcDedxValid()) { //Dedx information 
00361       RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
00362       m_chie[ii]  = dedxTrk->chiE();
00363       m_chimu[ii] = dedxTrk->chiMu();
00364       m_chipi[ii] = dedxTrk->chiPi();
00365       m_chik[ii]  = dedxTrk->chiK();
00366       m_chip[ii]  = dedxTrk->chiP();
00367       m_ghit[ii]  = dedxTrk->numGoodHits();
00368       m_thit[ii]  = dedxTrk->numTotalHits();
00369       m_probPH[ii]= dedxTrk->probPH();
00370       m_normPH[ii]= dedxTrk->normPH();
00371     }
00372 
00373    
00374     // TOF information 
00375     log << MSG::INFO <<"TOF information"<<endreq;
00376     double ptrk = mdcKalTrk->p() ; 
00377     if((*itTrk)->isTofTrackValid()) {
00378 
00379       SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
00380       
00381       SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
00382       for(;iter_tof != tofTrkCol.end(); iter_tof++ ) { 
00383         TofHitStatus *status = new TofHitStatus; 
00384         status->setStatus((*iter_tof)->status());
00385 
00386         if((status->is_barrel())){//barrel
00387           if( !(status->is_counter()) ) continue; // ? 
00388           if(status->layer()==1){ //layer1
00389             double path=(*iter_tof)->path(); // ? 
00390             double tof  = (*iter_tof)->tof();
00391             double ph   = (*iter_tof)->ph();
00392             double rhit = (*iter_tof)->zrhit();
00393             double qual = 0.0 + (*iter_tof)->quality();
00394             double cntr = 0.0 + (*iter_tof)->tofID();
00395             double texp[5];
00396             for(int j = 0; j < 5; j++) {
00397               double gb = ptrk/xmass[j];
00398               double beta = gb/sqrt(1+gb*gb);
00399               texp[j] = path /beta/velc;
00400             }
00401             
00402             m_cntr_btof1[ii]  = cntr;
00403             m_ptot_btof1[ii] = ptrk;
00404             m_ph_btof1[ii]   = ph;
00405             m_zhit_btof1[ii]  = rhit;
00406             m_qual_btof1[ii]  = qual;
00407             m_tof_btof1[ii]   = tof ; 
00408             m_te_btof1[ii]    = tof - texp[0];
00409             m_tmu_btof1[ii]   = tof - texp[1];
00410             m_tpi_btof1[ii]   = tof - texp[2];
00411             m_tk_btof1[ii]    = tof - texp[3];
00412             m_tp_btof1[ii]    = tof - texp[4];
00413           }
00414 
00415           if(status->layer()==2){//layer2
00416             double path=(*iter_tof)->path(); // ? 
00417             double tof  = (*iter_tof)->tof();
00418             double ph   = (*iter_tof)->ph();
00419             double rhit = (*iter_tof)->zrhit();
00420             double qual = 0.0 + (*iter_tof)->quality();
00421             double cntr = 0.0 + (*iter_tof)->tofID();
00422             double texp[5];
00423             for(int j = 0; j < 5; j++) {
00424               double gb = ptrk/xmass[j];
00425               double beta = gb/sqrt(1+gb*gb);
00426               texp[j] = path /beta/velc;
00427             }
00428  
00429             m_cntr_btof2[ii]  = cntr;
00430             m_ptot_btof2[ii] = ptrk;
00431             m_ph_btof2[ii]   = ph;
00432             m_zhit_btof2[ii]  = rhit;
00433             m_qual_btof2[ii]  = qual;
00434             m_tof_btof2[ii]   = tof ; 
00435             m_te_btof2[ii]    = tof - texp[0];
00436             m_tmu_btof2[ii]   = tof - texp[1];
00437             m_tpi_btof2[ii]   = tof - texp[2];
00438             m_tk_btof2[ii]    = tof - texp[3];
00439             m_tp_btof2[ii]    = tof - texp[4];
00440           } 
00441         }
00442         delete status; 
00443       } 
00444     }
00445     
00446     // EMC information 
00447     log << MSG::INFO <<"EMC information"<<endreq;
00448     if((*itTrk)->isEmcShowerValid()) {
00449       RecEmcShower *emcTrk = (*itTrk)->emcShower();
00450       m_e_emc[ii] = emcTrk->energy();
00451     }
00452  
00453   }
00454   counter[2]++;
00455   log << MSG::INFO <<"counter[2]="<<counter[2]<<endreq;
00456 
00457   double deltatof = 0.0 ; 
00458   log << MSG::INFO <<"11111111111"<<endreq;
00459   if(m_tof_btof1[0]*m_tof_btof1[1] != 0.0) deltatof+=fabs(m_tof_btof1[0]-m_tof_btof1[1]);
00460   log << MSG::INFO <<"22222222222"<<endreq;
00461   if(m_tof_btof2[0]*m_tof_btof2[1] != 0.0) deltatof+=fabs(m_tof_btof2[0]-m_tof_btof2[1]);
00462   log << MSG::INFO <<"33333333333"<<endreq;
00463   m_deltatof = deltatof ;  
00464 
00465   //boosted   mdcTrk
00466   p_pp[0].boost(u_lab);
00467   p_pm[0].boost(u_lab);
00468   for (int i=0; i<m_ngch; i++ ) {
00469     HepLorentzVector p; 
00470     if (i==0) p = p_pp[0]; 
00471     if (i==1) p = p_pm[0]; 
00472     
00473     m_bst_p[i]  = p.rho(); 
00474     m_bst_cos[i]= p.cosTheta(); 
00475   }
00476   log << MSG::INFO <<"44444444444"<<endreq;
00477 
00478   m_mpp = (p_pp[0] + p_pm[0]).m();
00479   m_angle = (p_pp[0].vect()).angle((p_pm[0]).vect()) * 180.0/(CLHEP::pi) ;
00480 
00481   counter[3]++;
00482   log << MSG::INFO <<"counter[3]="<<counter[3]<<endreq;
00483 
00484 
00485   // Vertex Fit
00486     
00487   // Kinamatic Fit
00488       
00489   // finale selection
00490   if (m_e_emc[0]>0.7) {return StatusCode::SUCCESS ; }       
00491   if (m_angle<178.0)    {return StatusCode::SUCCESS ; }
00492   if (fabs(m_bst_p[0]-1.232)<0.02 ) {return StatusCode::SUCCESS ; }
00493   if (fabs(m_bst_p[1]-1.232)<0.02 ) {return StatusCode::SUCCESS ; }
00494   if (fabs(m_deltatof)>4.0) {return StatusCode::SUCCESS ; }
00495 
00496 
00497   // Vertex Fit
00498 
00499   // Kinamatic Fit
00500 
00501   // finale selection
00502   
00503   m_tuple->write();
00504 
00505 
00506   // DQA
00507   TH1 *h(0);
00508   if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hbst_p", h).isSuccess()) {
00509       h->Fill(m_bst_p[0]);
00510   } else {
00511       log << MSG::ERROR << "Couldn't retrieve hbst_p" << endreq;
00512   }
00513 
00514   if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hbst_cos", h).isSuccess()) {
00515       h->Fill(m_bst_cos[0]);
00516   } else {
00517       log << MSG::ERROR << "Couldn't retrieve hbst_cos" << endreq;
00518   }
00519 
00520   if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hmpp", h).isSuccess()) {
00521       h->Fill(m_mpp);
00522   } else {
00523       log << MSG::ERROR << "Couldn't retrieve hmpp" << endreq;
00524   }
00525 
00526   if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hangle", h).isSuccess()) {
00527       h->Fill(m_angle);
00528   } else {
00529       log << MSG::ERROR << "Couldn't retrieve hangle" << endreq;
00530   }
00531 
00532   if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hdeltatof", h).isSuccess()) {
00533       h->Fill(m_deltatof);
00534   } else {
00535       log << MSG::ERROR << "Couldn't retrieve hdeltatof" << endreq;
00536   }
00537 
00538   if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/he_emc1", h).isSuccess()) {
00539       h->Fill(m_e_emc[0]);
00540   } else {
00541       log << MSG::ERROR << "Couldn't retrieve he_emc1" << endreq;
00542   }
00543 
00544   if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/he_emc2", h).isSuccess()) {
00545       h->Fill(m_e_emc[0]);
00546   } else {
00547       log << MSG::ERROR << "Couldn't retrieve he_emc2" << endreq;
00548   }
00549 
00551   // DQA
00552   // set tag and quality
00553 
00554   // Pid: 1 - electron, 2 - muon, 3 - pion, 4 - kaon, 5 - proton
00555   
00556   (*(evtRecTrkCol->begin()+iptrk[0]))->setPartId(5);
00557   (*(evtRecTrkCol->begin()+imtrk[0]))->setPartId(5);
00558   // Quality: defined by whether dE/dx or TOF is used to identify particle
00559   // 0 - no dE/dx, no TOF (can be used for dE/dx and TOF calibration)
00560   // 1 - only dE/dx (can be used for TOF calibration)
00561   // 2 - only TOF (can be used for dE/dx calibration)
00562   // 3 - Both dE/dx and TOF
00563   (*(evtRecTrkCol->begin()+iptrk[0]))->setQuality(0);
00564   (*(evtRecTrkCol->begin()+imtrk[0]))->setQuality(0);
00565 
00566   // DQA
00567   // Add the line below at the end of execute(), (before return)
00568   //
00569   setFilterPassed(true);
00571   
00572   return StatusCode::SUCCESS;
00573 
00574 }


Public Member Functions
	DQAJpsi2PPbarAlg (const std::string &name, ISvcLocator *pSvcLocator)
	DQAJpsi2PPbarAlg (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	execute ()
StatusCode	execute ()
StatusCode	finalize ()
StatusCode	finalize ()
StatusCode	initialize ()
StatusCode	initialize ()
Private Attributes
NTuple::Item< double >	m_angle
NTuple::Item< double >	m_angle
NTuple::Array< double >	m_bst_cos
NTuple::Array< double >	m_bst_cos
NTuple::Array< double >	m_bst_p
NTuple::Array< double >	m_bst_p
NTuple::Array< long >	m_charge
NTuple::Array< long >	m_charge
NTuple::Array< double >	m_chie
NTuple::Array< double >	m_chie
NTuple::Array< double >	m_chik
NTuple::Array< double >	m_chik
NTuple::Array< double >	m_chimu
NTuple::Array< double >	m_chimu
NTuple::Array< double >	m_chip
NTuple::Array< double >	m_chip
NTuple::Array< double >	m_chipi
NTuple::Array< double >	m_chipi
NTuple::Array< double >	m_cntr_btof1
NTuple::Array< double >	m_cntr_btof1
NTuple::Array< double >	m_cntr_btof2
NTuple::Array< double >	m_cntr_btof2
NTuple::Array< double >	m_cos
NTuple::Array< double >	m_cos
double	m_cthcut
NTuple::Item< double >	m_deltatof
NTuple::Item< double >	m_deltatof
NTuple::Array< double >	m_e_emc
NTuple::Array< double >	m_e_emc
double	m_energyThreshold
NTuple::Item< long >	m_event
NTuple::Item< long >	m_event
double	m_gammaAngCut
NTuple::Array< double >	m_ghit
NTuple::Array< double >	m_ghit
NTuple::Item< double >	m_mpp
NTuple::Item< double >	m_mpp
NTuple::Item< long >	m_ngch
NTuple::Item< long >	m_ngch
NTuple::Array< double >	m_normPH
NTuple::Array< double >	m_normPH
NTuple::Item< long >	m_np
NTuple::Item< long >	m_np
NTuple::Item< long >	m_npb
NTuple::Item< long >	m_npb
NTuple::Array< double >	m_p
NTuple::Array< double >	m_p
NTuple::Array< double >	m_ph_btof1
NTuple::Array< double >	m_ph_btof1
NTuple::Array< double >	m_ph_btof2
NTuple::Array< double >	m_ph_btof2
NTuple::Array< double >	m_probPH
NTuple::Array< double >	m_probPH
NTuple::Array< double >	m_ptot_btof1
NTuple::Array< double >	m_ptot_btof1
NTuple::Array< double >	m_ptot_btof2
NTuple::Array< double >	m_ptot_btof2
NTuple::Array< double >	m_px
NTuple::Array< double >	m_px
NTuple::Array< double >	m_py
NTuple::Array< double >	m_py
NTuple::Array< double >	m_pz
NTuple::Array< double >	m_pz
NTuple::Array< double >	m_qual_btof1
NTuple::Array< double >	m_qual_btof1
NTuple::Array< double >	m_qual_btof2
NTuple::Array< double >	m_qual_btof2
NTuple::Item< long >	m_runNo
NTuple::Item< long >	m_runNo
NTuple::Array< double >	m_te_btof1
NTuple::Array< double >	m_te_btof1
NTuple::Array< double >	m_te_btof2
NTuple::Array< double >	m_te_btof2
NTuple::Array< double >	m_thit
NTuple::Array< double >	m_thit
ITHistSvc *	m_thsvc
ITHistSvc *	m_thsvc
NTuple::Array< double >	m_tk_btof1
NTuple::Array< double >	m_tk_btof1
NTuple::Array< double >	m_tk_btof2
NTuple::Array< double >	m_tk_btof2
NTuple::Array< double >	m_tmu_btof1
NTuple::Array< double >	m_tmu_btof1
NTuple::Array< double >	m_tmu_btof2
NTuple::Array< double >	m_tmu_btof2
NTuple::Array< double >	m_tof_btof1
NTuple::Array< double >	m_tof_btof1
NTuple::Array< double >	m_tof_btof2
NTuple::Array< double >	m_tof_btof2
NTuple::Array< double >	m_tp_btof1
NTuple::Array< double >	m_tp_btof1
NTuple::Array< double >	m_tp_btof2
NTuple::Array< double >	m_tp_btof2
NTuple::Array< double >	m_tpi_btof1
NTuple::Array< double >	m_tpi_btof1
NTuple::Array< double >	m_tpi_btof2
NTuple::Array< double >	m_tpi_btof2
NTuple::Tuple *	m_tuple
NTuple::Tuple *	m_tuple
NTuple::Array< double >	m_vr0
NTuple::Array< double >	m_vr0
double	m_vr0cut
double	m_vr1cut
NTuple::Array< double >	m_vx0
NTuple::Array< double >	m_vx0
NTuple::Array< double >	m_vy0
NTuple::Array< double >	m_vy0
NTuple::Array< double >	m_vz0
NTuple::Array< double >	m_vz0
double	m_vz0cut
double	m_vz1cut
NTuple::Array< double >	m_zhit_btof1
NTuple::Array< double >	m_zhit_btof1
NTuple::Array< double >	m_zhit_btof2
NTuple::Array< double >	m_zhit_btof2

DQAJpsi2PPbarAlg Class Reference

Public Member Functions

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation