00001 #include <vector>
00002
00003 #include "GaudiKernel/MsgStream.h"
00004 #include "GaudiKernel/AlgFactory.h"
00005 #include "GaudiKernel/ISvcLocator.h"
00006 #include "GaudiKernel/SmartDataPtr.h"
00007 #include "GaudiKernel/IDataProviderSvc.h"
00008 #include "GaudiKernel/PropertyMgr.h"
00009 #include "VertexFit/IVertexDbSvc.h"
00010 #include "GaudiKernel/Bootstrap.h"
00011 #include "GaudiKernel/ISvcLocator.h"
00012
00013 #include "EventModel/EventModel.h"
00014 #include "EventModel/Event.h"
00015
00016 #include "EvtRecEvent/EvtRecEvent.h"
00017 #include "EvtRecEvent/EvtRecTrack.h"
00018 #include "DstEvent/TofHitStatus.h"
00019 #include "EventModel/EventHeader.h"
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
00030 using CLHEP::Hep3Vector;
00031 using CLHEP::Hep2Vector;
00032 using CLHEP::HepLorentzVector;
00033 #include "CLHEP/Geometry/Point3D.h"
00034
00035 #include "VertexFit/KinematicFit.h"
00036 #include "VertexFit/VertexFit.h"
00037 #include "ParticleID/ParticleID.h"
00038
00039
00040 #include "EvtSelExample/EvtSelExample.h"
00041
00042 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00043 typedef HepGeom::Point3D<double> HepPoint3D;
00044 #endif
00045 using CLHEP::HepLorentzVector;
00046
00047 const double mpi = 0.13957;
00048 const double mk = 0.493677;
00049 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00050 const double velc = 299.792458;
00051 typedef std::vector<int> Vint;
00052 typedef std::vector<HepLorentzVector> Vp4;
00053
00055
00056 EvtSelExample::EvtSelExample(const std::string& name, ISvcLocator* pSvcLocator) :
00057 Algorithm(name, pSvcLocator) {
00058
00059
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 }
00068
00069
00070
00071 StatusCode EvtSelExample::initialize(){
00072 MsgStream log(msgSvc(), name());
00073
00074 log << MSG::INFO << "in initialize()" << endmsg;
00075 StatusCode status;
00076
00077
00078
00079
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 }
00097
00098
00099 StatusCode EvtSelExample::execute() {
00100
00101 MsgStream log(msgSvc(), name());
00102 log << MSG::INFO << "in execute()" << endreq;
00103
00104
00105
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
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
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
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
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
00211
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
00226
00227
00228
00229 iGam.push_back((*itTrk)->trackId());
00230 }
00231
00232
00233
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
00243
00244
00245
00246
00247
00248 m_tuple->write();
00249
00251
00252
00253
00254
00255
00256 (*(evtRecTrkCol->begin()+ipip[0]))->setPartId(2);
00257 (*(evtRecTrkCol->begin()+ipim[0]))->setPartId(2);
00258
00259
00260
00261
00262
00263 (*(evtRecTrkCol->begin()+ipip[0]))->setQuality(0);
00264 (*(evtRecTrkCol->begin()+ipim[0]))->setQuality(0);
00265
00266
00267
00268
00269 setFilterPassed(true);
00271
00272 return StatusCode::SUCCESS;
00273
00274 }
00275
00276
00277 StatusCode EvtSelExample::finalize() {
00278
00279 MsgStream log(msgSvc(), name());
00280 log << MSG::INFO << "in finalize()" << endmsg;
00281 return StatusCode::SUCCESS;
00282 }
00283
00284