#include <EvtSelExample.h>
Public Member Functions | |
EvtSelExample (const std::string &name, ISvcLocator *pSvcLocator) | |
StatusCode | initialize () |
StatusCode | execute () |
StatusCode | finalize () |
Private Attributes | |
double | m_vr0cut |
double | m_vz0cut |
double | m_vr1cut |
double | m_vz1cut |
double | m_cthcut |
double | m_energyThreshold |
double | m_gammaAngCut |
NTuple::Tuple * | m_tuple |
NTuple::Item< long > | m_runNo |
NTuple::Item< long > | m_event |
Definition at line 22 of file EvtSelExample.h.
EvtSelExample::EvtSelExample | ( | const std::string & | name, | |
ISvcLocator * | pSvcLocator | |||
) |
Definition at line 56 of file EvtSelExample.cxx.
References m_cthcut, m_energyThreshold, m_gammaAngCut, m_vr0cut, m_vr1cut, m_vz0cut, and m_vz1cut.
00056 : 00057 Algorithm(name, pSvcLocator) { 00058 00059 //Declare the properties 00060 declareProperty("Vr0cut", m_vr0cut=5.0); 00061 declareProperty("Vz0cut", m_vz0cut=20.0); 00062 declareProperty("Vr1cut", m_vr1cut=1.0); 00063 declareProperty("Vz1cut", m_vz1cut=5.0); 00064 declareProperty("Vctcut", m_cthcut=0.93); 00065 declareProperty("EnergyThreshold", m_energyThreshold=0.04); 00066 declareProperty("GammaAngCut", m_gammaAngCut=20.0); 00067 }
StatusCode EvtSelExample::execute | ( | ) |
Definition at line 99 of file EvtSelExample.cxx.
References DstMdcTrack::charge(), cos(), Bes_Common::DEBUG, DstExtTrack::emcPosition(), DstExtTrack::emcVolumeNumber(), DstEmcShower::energy(), EventModel::EvtRec::EvtRecEvent, EventModel::EvtRec::EvtRecTrackCol, DstMdcTrack::helix(), genRecEmupikp::i, Bes_Common::INFO, IVertexDbSvc::isVertexValid(), ganga-rec::j, m_energyThreshold, m_event, m_gammaAngCut, m_runNo, m_tuple, m_vr0cut, m_vz0cut, msgSvc(), DstMdcTrack::p(), phi0, pi, IVertexDbSvc::PrimaryVertex(), IVertexDbSvc::SigmaPrimaryVertex(), sin(), DstEmcShower::x(), DstMdcTrack::x(), DstEmcShower::y(), DstMdcTrack::y(), DstEmcShower::z(), and DstMdcTrack::z().
00099 { 00100 00101 MsgStream log(msgSvc(), name()); 00102 log << MSG::INFO << "in execute()" << endreq; 00103 00104 // DQA 00105 // Add the line below at the beginning of execute() 00106 // 00107 setFilterPassed(false); 00108 00109 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader"); 00110 m_runNo=eventHeader->runNumber(); 00111 m_event=eventHeader->eventNumber(); 00112 log << MSG::DEBUG <<"run, evtnum = " 00113 << m_runNo << " , " 00114 << m_event <<endreq; 00115 00116 00117 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent); 00118 log << MSG::DEBUG <<"ncharg, nneu, tottks = " 00119 << evtRecEvent->totalCharged() << " , " 00120 << evtRecEvent->totalNeutral() << " , " 00121 << evtRecEvent->totalTracks() <<endreq; 00122 00123 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol); 00124 00125 if(evtRecEvent->totalNeutral()>100) { 00126 return StatusCode::SUCCESS; 00127 } 00128 00129 Vint iGood, ipip, ipim; 00130 iGood.clear(); 00131 ipip.clear(); 00132 ipim.clear(); 00133 Vp4 ppip, ppim; 00134 ppip.clear(); 00135 ppim.clear(); 00136 00137 Hep3Vector xorigin(0,0,0); 00138 00139 IVertexDbSvc* vtxsvc; 00140 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc); 00141 if(vtxsvc->isVertexValid()){ 00142 double* dbv = vtxsvc->PrimaryVertex(); 00143 double* vv = vtxsvc->SigmaPrimaryVertex(); 00144 xorigin.setX(dbv[0]); 00145 xorigin.setY(dbv[1]); 00146 xorigin.setZ(dbv[2]); 00147 } 00148 00149 int nCharge = 0; 00150 for(int i = 0; i < evtRecEvent->totalCharged(); i++){ 00151 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i; 00152 if(!(*itTrk)->isMdcTrackValid()) continue; 00153 if (!(*itTrk)->isMdcKalTrackValid()) continue; 00154 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack(); 00155 00156 double pch =mdcTrk->p(); 00157 double x0 =mdcTrk->x(); 00158 double y0 =mdcTrk->y(); 00159 double z0 =mdcTrk->z(); 00160 double phi0=mdcTrk->helix(1); 00161 double xv=xorigin.x(); 00162 double yv=xorigin.y(); 00163 double Rxy=fabs((x0-xv)*cos(phi0)+(y0-yv)*sin(phi0)); 00164 00165 if(fabs(z0) >= m_vz0cut) continue; 00166 if(Rxy >= m_vr0cut) continue; 00167 // if(fabs(m_Vct)>=m_cthcut) continue; 00168 iGood.push_back((*itTrk)->trackId()); 00169 nCharge += mdcTrk->charge(); 00170 if (mdcTrk->charge() > 0) { 00171 ipip.push_back((*itTrk)->trackId()); 00172 } else { 00173 ipim.push_back((*itTrk)->trackId()); 00174 } 00175 } 00176 00177 // 00178 // Finish Good Charged Track Selection 00179 // 00180 int nGood = iGood.size(); 00181 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq; 00182 if((nGood != 2)||(nCharge!=0)){ 00183 return StatusCode::SUCCESS; 00184 } 00185 00186 Vint iGam; 00187 iGam.clear(); 00188 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) { 00189 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i; 00190 if(!(*itTrk)->isEmcShowerValid()) continue; 00191 RecEmcShower *emcTrk = (*itTrk)->emcShower(); 00192 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z()); 00193 // find the nearest charged track 00194 double dthe = 200.; 00195 double dphi = 200.; 00196 double dang = 200.; 00197 for(int j = 0; j < evtRecEvent->totalCharged(); j++) { 00198 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j; 00199 if(!(*jtTrk)->isExtTrackValid()) continue; 00200 RecExtTrack *extTrk = (*jtTrk)->extTrack(); 00201 if(extTrk->emcVolumeNumber() == -1) continue; 00202 Hep3Vector extpos = extTrk->emcPosition(); 00203 // double ctht = extpos.cosTheta(emcpos); 00204 double angd = extpos.angle(emcpos); 00205 double thed = extpos.theta() - emcpos.theta(); 00206 double phid = extpos.deltaPhi(emcpos); 00207 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi; 00208 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi; 00209 00210 // if(fabs(thed) < fabs(dthe)) dthe = thed; 00211 // if(fabs(phid) < fabs(dphi)) dphi = phid; 00212 if(angd < dang) { 00213 dang = angd; 00214 dthe = thed; 00215 dphi = phid; 00216 } 00217 } 00218 if(dang>=200) continue; 00219 double eraw = emcTrk->energy(); 00220 dthe = dthe * 180 / (CLHEP::pi); 00221 dphi = dphi * 180 / (CLHEP::pi); 00222 dang = dang * 180 / (CLHEP::pi); 00223 if(eraw < m_energyThreshold) continue; 00224 if(dang < m_gammaAngCut) continue; 00225 // if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue; 00226 // 00227 // good photon cut will be set here 00228 // 00229 iGam.push_back((*itTrk)->trackId()); 00230 } 00231 00232 // 00233 // Finish Good Photon Selection 00234 // 00235 int nGam = iGam.size(); 00236 00237 log << MSG::DEBUG << "num Good Photon " << nGam << " , " <<evtRecEvent->totalNeutral()<<endreq; 00238 if(nGam<2){ 00239 return StatusCode::SUCCESS; 00240 } 00241 00242 // Vertex Fit 00243 00244 // Kinamatic Fit 00245 00246 // finale selection 00247 00248 m_tuple->write(); 00249 00251 // DQA 00252 // set tag and quality 00253 00254 // Pid: 1 - electron, 2 - muon, 3 - pion, 4 - kaon, 5 - proton 00255 00256 (*(evtRecTrkCol->begin()+ipip[0]))->setPartId(2); 00257 (*(evtRecTrkCol->begin()+ipim[0]))->setPartId(2); 00258 // Quality: defined by whether dE/dx or TOF is used to identify particle 00259 // 0 - no dE/dx, no TOF (can be used for dE/dx and TOF calibration) 00260 // 1 - only dE/dx (can be used for TOF calibration) 00261 // 2 - only TOF (can be used for dE/dx calibration) 00262 // 3 - Both dE/dx and TOF 00263 (*(evtRecTrkCol->begin()+ipip[0]))->setQuality(0); 00264 (*(evtRecTrkCol->begin()+ipim[0]))->setQuality(0); 00265 00266 // DQA 00267 // Add the line below at the end of execute(), (before return) 00268 // 00269 setFilterPassed(true); 00271 00272 return StatusCode::SUCCESS; 00273 00274 }
StatusCode EvtSelExample::finalize | ( | ) |
Definition at line 277 of file EvtSelExample.cxx.
References Bes_Common::INFO, and msgSvc().
00277 { 00278 00279 MsgStream log(msgSvc(), name()); 00280 log << MSG::INFO << "in finalize()" << endmsg; 00281 return StatusCode::SUCCESS; 00282 }
StatusCode EvtSelExample::initialize | ( | ) |
Definition at line 71 of file EvtSelExample.cxx.
References calibUtil::ERROR, Bes_Common::INFO, m_event, m_runNo, m_tuple, msgSvc(), and ntupleSvc().
00071 { 00072 MsgStream log(msgSvc(), name()); 00073 00074 log << MSG::INFO << "in initialize()" << endmsg; 00075 StatusCode status; 00076 00077 // DQA 00078 // The first directory specifier must be "DQAFILE" 00079 // The second is the control sample name, the first letter is in upper format. eg. "Rhopi" 00080 NTuplePtr nt(ntupleSvc(), "DQAFILE/Rhopi"); 00081 if ( nt ) m_tuple = nt; 00082 else { 00083 m_tuple = ntupleSvc()->book("DQAFILE/Rhopi", CLID_ColumnWiseTuple, "Rhopi ntuple"); 00084 if( m_tuple ) { 00085 status = m_tuple->addItem("runNo", m_runNo); 00086 status = m_tuple->addItem("event", m_event); 00087 } else { 00088 log << MSG::ERROR << "Can not book N-tuple:" << long(m_tuple) << endreq; 00089 } 00090 } 00091 00092 log << MSG::INFO << "successfully return from initialize()" <<endmsg; 00093 return StatusCode::SUCCESS; 00094 00095 00096 }
double EvtSelExample::m_cthcut [private] |
double EvtSelExample::m_energyThreshold [private] |
NTuple::Item<long> EvtSelExample::m_event [private] |
double EvtSelExample::m_gammaAngCut [private] |
NTuple::Item<long> EvtSelExample::m_runNo [private] |
NTuple::Tuple* EvtSelExample::m_tuple [private] |
double EvtSelExample::m_vr0cut [private] |
double EvtSelExample::m_vr1cut [private] |
double EvtSelExample::m_vz0cut [private] |
double EvtSelExample::m_vz1cut [private] |