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
00008 #include "EventModel/EventModel.h"
00009 #include "EventModel/Event.h"
00010
00011 #include "EvtRecEvent/EvtRecEvent.h"
00012 #include "EvtRecEvent/EvtRecTrack.h"
00013 #include "DstEvent/TofHitStatus.h"
00014 #include "EventModel/EventHeader.h"
00015
00016 #include "TMath.h"
00017 #include "GaudiKernel/INTupleSvc.h"
00018 #include "GaudiKernel/NTuple.h"
00019 #include "GaudiKernel/Bootstrap.h"
00020
00021 #include "GaudiKernel/ITHistSvc.h"
00022
00023 #include "CLHEP/Vector/ThreeVector.h"
00024 #include "CLHEP/Vector/LorentzVector.h"
00025 #include "CLHEP/Vector/TwoVector.h"
00026 using CLHEP::Hep3Vector;
00027 using CLHEP::Hep2Vector;
00028 using CLHEP::HepLorentzVector;
00029 #include "CLHEP/Geometry/Point3D.h"
00030 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00031 typedef HepGeom::Point3D<double> HepPoint3D;
00032 #endif
00033
00034 #include "VertexFit/KinematicFit.h"
00035 #include "VertexFit/VertexFit.h"
00036 #include "ParticleID/ParticleID.h"
00037
00038 #include "DQAinclPhi/inclphi.h"
00039
00040 #include <vector>
00041
00042 const double ecms = 3.097;
00043 const double mk = 0.493677;
00044 const double mphi = 1.01946;
00045 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00046 const double velc = 299.792458;
00047
00048
00049 typedef std::vector<int> Vint;
00050 typedef std::vector<HepLorentzVector> Vp4;
00051
00053
00054
00055
00056
00057
00059 inclphi::inclphi(const std::string& name, ISvcLocator* pSvcLocator) :
00060 Algorithm(name, pSvcLocator) {
00061
00062 m_tuple1 = 0;
00063 m_tuple2 = 0;
00064 for(int i = 0; i < 10; i++) m_pass[i] = 0;
00065
00066
00067 declareProperty("Vr0cut", m_vr0cut=5.0);
00068 declareProperty("Vz0cut", m_vz0cut=10.0);
00069 declareProperty("CheckDedx", m_checkDedx = 1);
00070 declareProperty("CheckTof", m_checkTof = 1);
00071 }
00072
00073
00074 StatusCode inclphi::initialize(){
00075 MsgStream log(msgSvc(), name());
00076
00077 log << MSG::INFO << "in initialize()" << endmsg;
00078
00079 StatusCode status;
00080
00081 if(service("THistSvc", m_thsvc).isFailure()) {
00082 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
00083 return StatusCode::FAILURE;
00084 }
00085
00086 TH1F* inclphi_mass = new TH1F("InclPhi_mass","INCLUSIVE_PHI_MASS",80,1.01,1.05);
00087 inclphi_mass->GetXaxis()->SetTitle("M_{KK} (GeV)");
00088 inclphi_mass->GetYaxis()->SetTitle("Nentries/0.5MeV/c^{2}");
00089 if(m_thsvc->regHist("/DQAHist/InclPhi/InclPhi_mass", inclphi_mass).isFailure()) {
00090 log << MSG::ERROR << "Couldn't register inclphi_mass" << endreq;
00091 }
00092
00093
00094
00095 NTuplePtr nt1(ntupleSvc(), "DQAFILE/InclPhi");
00096 if ( nt1 ) m_tuple1 = nt1;
00097 else {
00098 m_tuple1 = ntupleSvc()->book ("DQAFILE/InclPhi1", CLID_ColumnWiseTuple, "inclphi Ntuple");
00099 if ( m_tuple1 ) {
00100 status = m_tuple1->addItem ("mphiall", m_mphiall);
00101 status = m_tuple1->addItem ("mpphiall", m_pphiall);
00102 }
00103 else {
00104 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00105 return StatusCode::FAILURE;
00106 }
00107 }
00108 NTuplePtr nt2(ntupleSvc(), "DQAFILE/InclPhi2");
00109 if ( nt2 ) m_tuple2 = nt2;
00110 else {
00111 m_tuple2 = ntupleSvc()->book ("DQAFILE/InclPhi2", CLID_ColumnWiseTuple, "inclphi Ntuple");
00112 if ( m_tuple2 ) {
00113 status = m_tuple2->addItem ("nkp", m_nkp);
00114 status = m_tuple2->addItem ("nkm", m_nkm);
00115 status = m_tuple2->addItem ("ncharge", m_ncharge);
00116 status = m_tuple2->addItem ("difchi0", m_difchi);
00117 status = m_tuple2->addItem ("mmphi", m_mphi);
00118 status = m_tuple2->addItem ("pphi", m_pphi);
00119 status = m_tuple2->addItem ("pk1", m_pk1);
00120 }
00121 else {
00122 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
00123 return StatusCode::FAILURE;
00124 }
00125 }
00126
00127
00128
00129
00130 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00131 return StatusCode::SUCCESS;
00132
00133 }
00134
00135
00136 StatusCode inclphi::execute() {
00137 StatusCode sc = StatusCode::SUCCESS;
00138
00139 MsgStream log(msgSvc(), name());
00140 log << MSG::INFO << "in execute()" << endreq;
00141
00142
00143
00144
00145 setFilterPassed(false);
00146
00147 m_pass[0] += 1;
00148
00149 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00150 int run=eventHeader->runNumber();
00151 int event=eventHeader->eventNumber();
00152
00153 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00154 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00155
00156 m_pass[1] += 1;
00157
00158 Vint ikp, ikm, iGood;
00159 iGood.clear();
00160 ikp.clear();
00161 ikm.clear();
00162
00163 Vp4 pkp, pkm;
00164 pkp.clear();
00165 pkm.clear();
00166
00167 int TotCharge = 0;
00168 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00169 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00170 if(!(*itTrk)->isMdcTrackValid()) continue;
00171 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00172 if(fabs(mdcTrk->z()) >= m_vz0cut) continue;
00173 if(mdcTrk->r() >= m_vr0cut) continue;
00174 iGood.push_back(i);
00175 TotCharge += mdcTrk->charge();
00176 }
00177
00178
00179
00180 int nGood = iGood.size();
00181
00182
00183
00184
00185
00186 if((nGood < 2) || (TotCharge!=0)) return sc;
00187
00188 m_pass[2] += 1;
00189
00190
00191
00192
00193 ParticleID *pid = ParticleID::instance();
00194 for(int i = 0; i < nGood; i++) {
00195 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00196
00197 pid->init();
00198 pid->setMethod(pid->methodProbability());
00199 pid->setChiMinCut(4);
00200 pid->setRecTrack(*itTrk);
00201 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
00202 pid->identify(pid->onlyPion() | pid->onlyKaon());
00203 pid->calculate();
00204 if(!(pid->IsPidInfoValid())) continue;
00205 if(!(*itTrk)->isMdcKalTrackValid()) continue ;
00206 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
00207
00208
00209 if(pid->probKaon() < 0.001 || (pid->probKaon() < pid->probPion())) continue;
00210
00211 mdcKalTrk->setPidType(RecMdcKalTrack::kaon);
00212 HepLorentzVector ptrk;
00213 ptrk.setPx(mdcKalTrk->px());
00214 ptrk.setPy(mdcKalTrk->py());
00215 ptrk.setPz(mdcKalTrk->pz());
00216 double p3 = ptrk.mag();
00217 ptrk.setE(sqrt(p3*p3+mk*mk));
00218 if(mdcKalTrk->charge() >0) {
00219 ikp.push_back(iGood[i]);
00220 pkp.push_back(ptrk);
00221 } else {
00222 ikm.push_back(iGood[i]);
00223 pkm.push_back(ptrk);
00224 }
00225 }
00226
00227 m_pass[4] += 1;
00228 int nkp = ikp.size();
00229 int nkm = ikm.size();
00230
00231 m_nkp=nkp;
00232 m_nkm=nkm;
00233 m_ncharge=nGood;
00234
00235
00236 if(nkp < 1 || nkm <1) return sc;
00237
00238 m_pass[5] += 1;
00239
00240
00241
00242
00243
00244
00245 HepLorentzVector pphi, pTot;
00246 Vint iphi;
00247 iphi.clear();
00248
00249 double difchi0=99999.0;
00250 int ixk1 = -1;
00251 int ixk2 = -1;
00252
00253 for (int i = 0; i < nkm; i++) {
00254 for (int j = 0; j < nkp; j++) {
00255
00256 pphi = pkm[i] + pkp[j];
00257 m_mphiall = pphi.m();
00258 m_pphiall = pphi.rho();
00259 m_tuple1->write();
00260
00261 double difchi =fabs(pphi.m()-mphi);
00262 if(difchi < difchi0){
00263 difchi0 = difchi;
00264 ixk1 = i;
00265 ixk2 = j;
00266 }
00267 }
00268 }
00269
00270 m_difchi = difchi0;
00271
00272 if (ixk1 == -1) return sc;
00273
00274 m_pass[6] += 1;
00275
00276 pTot = pkm[ixk1] + pkp[ixk2];
00277
00278 m_pk1 = pkm[ixk1].m();
00279 m_mphi = pTot.m();
00280 m_pphi = pTot.rho();
00281
00282 TH1 *h(0);
00283 if (m_thsvc->getHist("/DQAHist/InclPhi/InclPhi_mass", h).isSuccess()) {
00284 h->Fill(pTot.m());
00285 } else {
00286 log << MSG::ERROR << "Couldn't retrieve inclphi_mass" << endreq;
00287 }
00288 m_tuple2->write();
00290
00291
00292
00293
00294
00295 (*(evtRecTrkCol->begin()+ikm[ixk1]))->tagKaon();
00296 (*(evtRecTrkCol->begin()+ikp[ixk2]))->tagKaon();
00297
00298
00299
00300
00301
00302 (*(evtRecTrkCol->begin()+ikm[ixk1]))->setQuality(3);
00303 (*(evtRecTrkCol->begin()+ikp[ixk2]))->setQuality(3);
00304
00305
00306
00307 setFilterPassed(true);
00308
00309 return StatusCode::SUCCESS;
00310 }
00311
00312
00313 StatusCode inclphi::finalize() {
00314
00315 MsgStream log(msgSvc(), name());
00316 log << MSG::INFO << "in finalize()" << endmsg;
00317 log << MSG::INFO << "Total Entries : " << m_pass[0] << endreq;
00318 log << MSG::INFO << "dstCol, trkCol: " << m_pass[1] << endreq;
00319 log << MSG::INFO << "Ncharge Cut : " << m_pass[2] << endreq;
00320 log << MSG::INFO << "TOF dEdx : " << m_pass[3] << endreq;
00321 log << MSG::INFO << "PID : " << m_pass[4] << endreq;
00322 log << MSG::INFO << "NK Cut : " << m_pass[5] << endreq;
00323 log << MSG::INFO << "Nphi select : " << m_pass[6] << endreq;
00324 return StatusCode::SUCCESS;
00325 }
00326