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 "VertexFit/SecondVertexFit.h"
00037 #include "ParticleID/ParticleID.h"
00038
00039 #include "DQAinclLamda/incllambda.h"
00040
00041 #include <vector>
00042
00043 const double ecms = 3.097;
00044 const double mpi = 0.13957;
00045 const double mp = 0.93827203;
00046 const double mlam = 1.115683;
00047 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00048
00049
00050 typedef std::vector<int> Vint;
00051 typedef std::vector<HepLorentzVector> Vp4;
00052
00054
00055
00056
00057
00058
00060
00061 incllambda::incllambda(const std::string& name, ISvcLocator* pSvcLocator) :
00062 Algorithm(name, pSvcLocator) {
00063
00064 m_tuple1 = 0;
00065 for(int i = 0; i < 10; i++) m_pass[i] = 0;
00066
00067
00068 declareProperty("Vr0cut", m_vr0cut=20.0);
00069 declareProperty("Vz0cut", m_vz0cut=50.0);
00070 declareProperty("CheckDedx", m_checkDedx = 1);
00071 declareProperty("CheckTof", m_checkTof = 1);
00072
00073 }
00074
00075
00076 StatusCode incllambda::initialize(){
00077 MsgStream log(msgSvc(), name());
00078
00079 log << MSG::INFO << "in initialize()" << endmsg;
00080
00081 StatusCode status;
00082
00083 if(service("THistSvc", m_thsvc).isFailure()) {
00084 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
00085 return StatusCode::FAILURE;
00086 }
00087
00088 TH1F* incllambda_mass = new TH1F("InclLam_mass","INCLUSIVE_LAMBDA_MASS",150,1.11,1.125);
00089 incllambda_mass->GetXaxis()->SetTitle("M_{p#pi} (GeV)");
00090 incllambda_mass->GetYaxis()->SetTitle("Nentries/0.1MeV/c^{2}");
00091 if(m_thsvc->regHist("/DQAHist/InclLam/InclLam_mass", incllambda_mass).isFailure()) {
00092 log << MSG::ERROR << "Couldn't register incllambda_mass" << endreq;
00093 }
00094
00095
00096
00097 NTuplePtr nt1(ntupleSvc(), "DQAFILE/InclLam");
00098 if ( nt1 ) m_tuple1 = nt1;
00099 else {
00100 m_tuple1 = ntupleSvc()->book ("DQAFILE/InclLam", CLID_ColumnWiseTuple, "incllam Ntuple");
00101 if ( m_tuple1 ) {
00102 status = m_tuple1->addItem ("npi", m_npi);
00103 status = m_tuple1->addItem ("np", m_np);
00104 status = m_tuple1->addItem ("lxyz", m_lxyz);
00105 status = m_tuple1->addItem ("ctau", m_ctau);
00106 status = m_tuple1->addItem ("elxyz", m_elxyz);
00107 status = m_tuple1->addItem ("lamchi", m_chis);
00108 status = m_tuple1->addItem ("mlamraw", m_lamRawMass);
00109 status = m_tuple1->addItem ("plam", m_plam);
00110 }
00111 else {
00112 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00113 return StatusCode::FAILURE;
00114 }
00115 }
00116
00117
00118
00119
00120 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00121 return StatusCode::SUCCESS;
00122
00123 }
00124
00125
00126 StatusCode incllambda::execute() {
00127 StatusCode sc = StatusCode::SUCCESS;
00128
00129 MsgStream log(msgSvc(), name());
00130 log << MSG::INFO << "in execute()" << endreq;
00131
00132
00133
00134
00135 setFilterPassed(false);
00136
00137 m_pass[0] += 1;
00138
00139 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00140 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00141 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00142
00143 Vint ilam, ip, ipi, iGood;
00144 iGood.clear();
00145 ilam.clear();
00146 ip.clear();
00147 ipi.clear();
00148
00149 Vp4 ppi, pp;
00150 ppi.clear();
00151 pp.clear();
00152
00153 int TotCharge = 0;
00154 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00155 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00156 if(!(*itTrk)->isMdcTrackValid()) continue;
00157 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00158 if(fabs(mdcTrk->z()) >= m_vz0cut) continue;
00159 if(mdcTrk->r() >= m_vr0cut) continue;
00160 iGood.push_back(i);
00161 TotCharge += mdcTrk->charge();
00162 }
00163
00164
00165
00166 int nGood = iGood.size();
00167
00168
00169
00170
00171
00172 if((nGood < 2) || (TotCharge!=0)) return sc;
00173
00174 m_pass[1] += 1;
00175
00176
00177
00178 ParticleID *pid = ParticleID::instance();
00179 for(int i = 0; i < nGood; i++) {
00180 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00181
00182 pid->init();
00183 pid->setMethod(pid->methodProbability());
00184 pid->setChiMinCut(4);
00185 pid->setRecTrack(*itTrk);
00186 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
00187 pid->identify(pid->onlyPion() | pid->onlyProton());
00188 pid->calculate();
00189 if(!(pid->IsPidInfoValid())) continue;
00190
00191 if(!(*itTrk)->isMdcKalTrackValid()) continue ;
00192 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
00193
00194 if(mdcKalTrk->charge() == 0) continue;
00195
00196 if(pid->probPion() > 0.001 && (pid->probPion() > pid->probProton())){
00197 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
00198 ipi.push_back(iGood[i]);
00199 HepLorentzVector ptrk;
00200 ptrk.setPx(mdcKalTrk->px());
00201 ptrk.setPy(mdcKalTrk->py());
00202 ptrk.setPz(mdcKalTrk->pz());
00203 double p3 = ptrk.mag();
00204 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00205 ppi.push_back(ptrk);
00206 }
00207 if(pid->probProton() > 0.001 && (pid->probProton() > pid->probPion())){
00208 mdcKalTrk->setPidType(RecMdcKalTrack::proton);
00209 ip.push_back(iGood[i]);
00210 HepLorentzVector ptrk;
00211 ptrk.setPx(mdcKalTrk->px());
00212 ptrk.setPy(mdcKalTrk->py());
00213 ptrk.setPz(mdcKalTrk->pz());
00214 double p3 = ptrk.mag();
00215 ptrk.setE(sqrt(p3*p3+mp*mp));
00216 pp.push_back(ptrk);
00217 }
00218 }
00219
00220 m_pass[2] += 1;
00221
00222 int npi = ipi.size();
00223 int np = ip.size();
00224 m_npi = npi;
00225 m_np = np;
00226 if(npi < 1 || np <1) return sc;
00227
00228 m_pass[3] += 1;
00229
00230
00231
00232 double chi=999.;
00233 double delm;
00234
00235 HepPoint3D vx(0., 0., 0.);
00236 HepSymMatrix Evx(3, 0);
00237 double bx = 1E+6;
00238 double by = 1E+6;
00239 double bz = 1E+6;
00240 Evx[0][0] = bx*bx;
00241 Evx[1][1] = by*by;
00242 Evx[2][2] = bz*bz;
00243 VertexParameter vxpar;
00244 vxpar.setVx(vx);
00245 vxpar.setEvx(Evx);
00246
00247 int ipion=-1, iproton=-1;
00248
00249 VertexFit *vtxfit0 = VertexFit::instance();
00250 SecondVertexFit *vtxfit = SecondVertexFit::instance();
00251
00252 for(unsigned int i1 = 0; i1 < ipi.size(); i1++) {
00253 RecMdcKalTrack *piTrk = (*(evtRecTrkCol->begin()+ipi[i1]))->mdcKalTrack();
00254 piTrk->setPidType(RecMdcKalTrack::pion);
00255 WTrackParameter wpitrk(mpi, piTrk->helix(), piTrk->err());
00256
00257 for(unsigned int i2 = 0; i2 < ip.size(); i2++) {
00258 RecMdcKalTrack *protonTrk = (*(evtRecTrkCol->begin()+ip[i2]))->mdcKalTrack();
00259 protonTrk->setPidType(RecMdcKalTrack::proton);
00260 WTrackParameter wptrk(mp, protonTrk->helix(), protonTrk->err());
00261
00262 int NCharge = 0;
00263 NCharge =piTrk->charge() + protonTrk->charge();
00264
00265
00266
00267 if(NCharge != 0) continue;
00268
00269 vtxfit0->init();
00270 vtxfit0->AddTrack(0, wpitrk);
00271 vtxfit0->AddTrack(1, wptrk);
00272 vtxfit0->AddVertex(0, vxpar, 0, 1);
00273 if(!(vtxfit0->Fit(0))) continue;
00274 vtxfit0->BuildVirtualParticle(0);
00275
00276 vtxfit->init();
00277
00278 vtxfit->AddTrack(0, vtxfit0->wVirtualTrack(0));
00279 vtxfit->setVpar(vtxfit0->vpar(0));
00280 if(!vtxfit->Fit()) continue;
00281
00282
00283
00284
00285 if(vtxfit->chisq() > chi) continue;
00286
00287 delm = fabs((vtxfit->p4par()).m()-mlam);
00288 chi = vtxfit->chisq();
00289 ipion=ipi[i1];
00290 iproton=ip[i2];
00291 }
00292 }
00293
00294
00295
00296 if(ipion==-1 || iproton==-1) return sc;
00297 m_pass[4] += 1;
00298
00299 RecMdcKalTrack *pionTrk = (*(evtRecTrkCol->begin()+ipion))->mdcKalTrack();
00300 pionTrk->setPidType(RecMdcKalTrack::pion);
00301 WTrackParameter wpiontrk(mpi, pionTrk->helix(), pionTrk->err());
00302
00303 RecMdcKalTrack *protonTrk = (*(evtRecTrkCol->begin()+iproton))->mdcKalTrack();
00304 protonTrk->setPidType(RecMdcKalTrack::proton);
00305 WTrackParameter wprotontrk(mp, protonTrk->helix(), protonTrk->err());
00306 vtxfit0->init();
00307 vtxfit0->AddTrack(0, wpiontrk);
00308 vtxfit0->AddTrack(1, wprotontrk);
00309 vtxfit0->AddVertex(0, vxpar, 0, 1);
00310 if(!(vtxfit0->Fit(0))) return sc;
00311 vtxfit0->BuildVirtualParticle(0);
00312
00313 vtxfit->init();
00314
00315 vtxfit->AddTrack(0, vtxfit0->wVirtualTrack(0));
00316 vtxfit->setVpar(vtxfit0->vpar(0));
00317 if(!vtxfit->Fit()) return sc;
00318
00319 m_lamRawMass = vtxfit->p4par().m();
00320 m_ctau = vtxfit->ctau();
00321 m_lxyz = vtxfit->decayLength();
00322 m_elxyz = vtxfit->decayLengthError();
00323 m_chis = vtxfit->chisq();
00324 m_plam = vtxfit->p4par().rho();
00325
00326 if(m_lxyz>3.5){
00327
00328 TH1 *h(0);
00329 if (m_thsvc->getHist("/DQAHist/InclLam/InclLam_mass", h).isSuccess()) {
00330 h->Fill(m_lamRawMass);
00331 } else {
00332 log << MSG::ERROR << "Couldn't retrieve incllam_mass" << endreq;
00333 }
00334
00335 }
00336 m_tuple1->write();
00338
00339
00340
00341
00342
00343 (*(evtRecTrkCol->begin()+ipion))->tagPion();
00344 (*(evtRecTrkCol->begin()+iproton))->tagProton();
00345
00346
00347
00348
00349
00350 (*(evtRecTrkCol->begin()+ipion))->setQuality(1);
00351 (*(evtRecTrkCol->begin()+iproton))->setQuality(1);
00352
00353
00354
00355
00356 setFilterPassed(true);
00357
00358 return StatusCode::SUCCESS;
00359 }
00360
00361
00362 StatusCode incllambda::finalize() {
00363
00364 MsgStream log(msgSvc(), name());
00365 log << MSG::INFO << "in finalize()" << endmsg;
00366 log << MSG::INFO << "Total Entries : " << m_pass[0] << endreq;
00367 log << MSG::INFO << "Qtot Cut : " << m_pass[1] << endreq;
00368 log << MSG::INFO << "Nks : " << m_pass[2] << endreq;
00369 log << MSG::INFO << "Good Ks cut : " << m_pass[3] << endreq;
00370 return StatusCode::SUCCESS;
00371 }
00372