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 #include "GaudiKernel/Bootstrap.h"
00008
00009 #include "GaudiKernel/INTupleSvc.h"
00010 #include "GaudiKernel/NTuple.h"
00011 #include "GaudiKernel/ITHistSvc.h"
00012
00013
00014 #include "EventModel/EventModel.h"
00015 #include "EventModel/Event.h"
00016
00017 #include "EvtRecEvent/EvtRecEvent.h"
00018 #include "EvtRecEvent/EvtRecTrack.h"
00019 #include "DstEvent/TofHitStatus.h"
00020 #include "EventModel/EventHeader.h"
00021
00022 #include "TMath.h"
00023 #include "CLHEP/Vector/ThreeVector.h"
00024 #include "CLHEP/Vector/LorentzVector.h"
00025 #include "CLHEP/Vector/TwoVector.h"
00026
00027 using CLHEP::Hep3Vector;
00028 using CLHEP::Hep2Vector;
00029 using CLHEP::HepLorentzVector;
00030 #include "CLHEP/Geometry/Point3D.h"
00031 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00032 typedef HepGeom::Point3D<double> HepPoint3D;
00033 #endif
00034
00035 #include "VertexFit/IVertexDbSvc.h"
00036 #include "VertexFit/KinematicFit.h"
00037 #include "VertexFit/VertexFit.h"
00038 #include "ParticleID/ParticleID.h"
00039 #include "VertexFit/SecondVertexFit.h"
00040
00041
00042 #include "DQAinclKs/inclks.h"
00043
00044 #include <vector>
00045
00046 const double mk0 = 0.497648;
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 inclks::inclks(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=10.0);
00069 declareProperty("Vz0cut", m_vz0cut=50.0);
00070 declareProperty("CheckDedx", m_checkDedx = 1);
00071 declareProperty("CheckTof", m_checkTof = 1);
00072
00073 }
00074
00075 StatusCode inclks::initialize(){
00076 MsgStream log(msgSvc(), name());
00077
00078 log << MSG::INFO << "in initialize()" << endmsg;
00079
00080 StatusCode status;
00081
00082 if(service("THistSvc", m_thsvc).isFailure()) {
00083 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
00084 return StatusCode::FAILURE;
00085 }
00086
00087 TH1F* inclks_mass = new TH1F("InclKs_mass","INCLUSIVE_Ks_MASS",70,0.46,0.53);
00088 inclks_mass->GetXaxis()->SetTitle("M_{#pi#pi} (GeV)");
00089 inclks_mass->GetYaxis()->SetTitle("Nentries/0.1MeV/c^{2}");
00090 if(m_thsvc->regHist("/DQAHist/InclKs/InclKs_mass", inclks_mass).isFailure()) {
00091 log << MSG::ERROR << "Couldn't register inclks_mass" << endreq;
00092 }
00093
00094
00095
00096 NTuplePtr nt1(ntupleSvc(), "DQAFILE/InclKs");
00097 if ( nt1 ) m_tuple1 = nt1;
00098 else {
00099 m_tuple1 = ntupleSvc()->book ("DQAFILE/InclKs", CLID_ColumnWiseTuple, "ksklx Ntuple");
00100 if ( m_tuple1 ) {
00101 status = m_tuple1->addItem ("npip", m_npip);
00102 status = m_tuple1->addItem ("npim", m_npim);
00103 status = m_tuple1->addItem ("ctau", m_ctau);
00104 status = m_tuple1->addItem ("lxyz", m_lxyz);
00105 status = m_tuple1->addItem ("elxyz", m_elxyz);
00106 status = m_tuple1->addItem ("kschi", m_chis);
00107 status = m_tuple1->addItem ("mksraw", m_ksRawMass);
00108 status = m_tuple1->addItem ("pk0", m_pk0);
00109
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
00121 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00122 return StatusCode::SUCCESS;
00123
00124 }
00125
00126
00127 StatusCode inclks::execute() {
00128 StatusCode sc = StatusCode::SUCCESS;
00129
00130 MsgStream log(msgSvc(), name());
00131 log << MSG::INFO << "in execute()" << endreq;
00132
00133
00134
00135
00136 setFilterPassed(false);
00137
00138 m_pass[0] += 1;
00139
00140 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00141 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00142 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00143
00144 Vint iks, ipip, ipim, iGood;
00145 iGood.clear();
00146 iks.clear();
00147 ipip.clear();
00148 ipim.clear();
00149
00150 Vp4 ppip, ppim;
00151 ppip.clear();
00152 ppim.clear();
00153
00154 int TotCharge = 0;
00155 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00156 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00157 if(!(*itTrk)->isMdcTrackValid()) continue;
00158 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00159 if(fabs(mdcTrk->z()) >= m_vz0cut) continue;
00160 if(mdcTrk->r() >= m_vr0cut) continue;
00161 iGood.push_back(i);
00162 TotCharge += mdcTrk->charge();
00163 }
00164
00165
00166
00167 int nGood = iGood.size();
00168
00169
00170
00171
00172
00173 if((nGood < 2) || (TotCharge!=0)) return sc;
00174
00175 m_pass[1] += 1;
00176
00177
00178
00179 ParticleID *pid = ParticleID::instance();
00180 for(int i = 0; i < nGood; i++) {
00181 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00182
00183 pid->init();
00184 pid->setMethod(pid->methodProbability());
00185 pid->setChiMinCut(4);
00186 pid->setRecTrack(*itTrk);
00187 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
00188 pid->identify(pid->onlyPion() | pid->onlyKaon());
00189
00190
00191 pid->calculate();
00192 if(!(pid->IsPidInfoValid())) continue;
00193 if(!(*itTrk)->isMdcKalTrackValid()) continue ;
00194 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
00195
00196
00197 if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
00198
00199 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
00200 HepLorentzVector ptrk;
00201 ptrk.setPx(mdcKalTrk->px());
00202 ptrk.setPy(mdcKalTrk->py());
00203 ptrk.setPz(mdcKalTrk->pz());
00204 double p3 = ptrk.mag();
00205 ptrk.setE(sqrt(p3*p3+xmass[2]*xmass[2]));
00206
00207 if(mdcKalTrk->charge() >0) {
00208 ipip.push_back(iGood[i]);
00209 ppip.push_back(ptrk);
00210 } else {
00211 ipim.push_back(iGood[i]);
00212 ppim.push_back(ptrk);
00213 }
00214 }
00215
00216 m_pass[2] += 1;
00217 int npip = ipip.size();
00218 int npim = ipim.size();
00219 m_npip=npip;
00220 m_npim=npim;
00221
00222 if(npip < 1 || npim <1) return sc;
00223
00224 m_pass[3] += 1;
00225
00226
00227
00228
00229 double chi, delm;
00230 double chisq=999.;
00231
00232
00233 HepPoint3D vx(0., 0., 0.);
00234 HepSymMatrix Evx(3, 0);
00235 double bx = 1E+6;
00236 double by = 1E+6;
00237 double bz = 1E+6;
00238 Evx[0][0] = bx*bx;
00239 Evx[1][1] = by*by;
00240 Evx[2][2] = bz*bz;
00241 VertexParameter vxpar;
00242 vxpar.setVx(vx);
00243 vxpar.setEvx(Evx);
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254 int ip1=-1, ip2=-1;
00255
00256 VertexFit *vtxfit0 = VertexFit::instance();
00257 SecondVertexFit *vtxfit = SecondVertexFit::instance();
00258
00259 for(int i1 = 0; i1 < ipip.size(); i1++) {
00260 RecMdcKalTrack *pipTrk = (*(evtRecTrkCol->begin()+ipip[i1]))->mdcKalTrack();
00261 pipTrk->setPidType(RecMdcKalTrack::pion);
00262 WTrackParameter wpiptrk(xmass[2], pipTrk->helix(), pipTrk->err());
00263
00264 for(int i2 = 0; i2 < ipim.size(); i2++) {
00265 RecMdcKalTrack *pimTrk = (*(evtRecTrkCol->begin()+ipim[i2]))->mdcKalTrack();
00266 pimTrk->setPidType(RecMdcKalTrack::pion);
00267 WTrackParameter wpimtrk(xmass[2], pimTrk->helix(), pimTrk->err());
00268
00269 vtxfit0->init();
00270 vtxfit0->AddTrack(0, wpiptrk);
00271 vtxfit0->AddTrack(1, wpimtrk);
00272 vtxfit0->AddVertex(0, vxpar, 0, 1);
00273 if(!(vtxfit0->Fit(0))) continue;
00274 vtxfit0->BuildVirtualParticle(0);
00275
00276 vtxfit->init();
00277 vtxfit->AddTrack(0, vtxfit0->wVirtualTrack(0));
00278 vtxfit->setVpar(vtxfit0->vpar(0));
00279 if(!vtxfit->Fit()) continue;
00280 chi = vtxfit->chisq();
00281
00282
00283 if(chi > chisq) continue;
00284 delm = fabs((vtxfit->p4par().m())-mk0);
00285 chisq = chi;
00286 ip1=ipip[i1];
00287 ip2=ipim[i2];
00288 }
00289 }
00290
00291
00292
00293 if(ip1==-1 || ip2==-1) return sc;
00294 m_pass[4] += 1;
00295
00296 RecMdcKalTrack *pi1Trk = (*(evtRecTrkCol->begin()+ip1))->mdcKalTrack();
00297 pi1Trk->setPidType(RecMdcKalTrack::pion);
00298 WTrackParameter wpi1trk(xmass[2], pi1Trk->helix(), pi1Trk->err());
00299
00300 RecMdcKalTrack *pi2Trk = (*(evtRecTrkCol->begin()+ip2))->mdcKalTrack();
00301 pi2Trk->setPidType(RecMdcKalTrack::pion);
00302 WTrackParameter wpi2trk(xmass[2], pi2Trk->helix(), pi2Trk->err());
00303 vtxfit0->init();
00304 vtxfit0->AddTrack(0, wpi1trk);
00305 vtxfit0->AddTrack(1, wpi2trk);
00306 vtxfit0->AddVertex(0, vxpar, 0, 1);
00307 if(!(vtxfit0->Fit(0))) return sc;
00308 vtxfit0->BuildVirtualParticle(0);
00309
00310 vtxfit->init();
00311 vtxfit->AddTrack(0, vtxfit0->wVirtualTrack(0));
00312 vtxfit->setVpar(vtxfit0->vpar(0));
00313 if(!vtxfit->Fit()) return sc;
00314
00315 m_ksRawMass = vtxfit->p4par().m();
00316 m_ctau = vtxfit->ctau();
00317 m_lxyz = vtxfit->decayLength();
00318 m_elxyz = vtxfit->decayLengthError();
00319 m_chis = vtxfit->chisq();
00320 m_pk0 = vtxfit->p4par().rho();
00321
00322
00323 if(m_lxyz>0.4 && m_chis<20.){
00324
00325 TH1 *h(0);
00326 if (m_thsvc->getHist("/DQAHist/InclKs/InclKs_mass", h).isSuccess()) {
00327 h->Fill(m_ksRawMass);
00328 } else {
00329 log << MSG::ERROR << "Couldn't retrieve inclks_mass" << endreq;
00330 }
00331 }
00332
00333 m_tuple1->write();
00335
00336
00337
00338
00339
00340 (*(evtRecTrkCol->begin()+ip1))->tagPion();
00341 (*(evtRecTrkCol->begin()+ip2))->tagPion();
00342
00343
00344
00345
00346
00347 (*(evtRecTrkCol->begin()+ip1))->setQuality(2);
00348 (*(evtRecTrkCol->begin()+ip2))->setQuality(2);
00349
00350
00351
00352 setFilterPassed(true);
00353
00354 return StatusCode::SUCCESS;
00355 }
00356
00357
00358 StatusCode inclks::finalize() {
00359
00360 MsgStream log(msgSvc(), name());
00361 log << MSG::INFO << "in finalize()" << endmsg;
00362 log << MSG::INFO << "Total Entries : " << m_pass[0] << endreq;
00363 log << MSG::INFO << "Qtot Cut : " << m_pass[1] << endreq;
00364 log << MSG::INFO << "PID : " << m_pass[2] << endreq;
00365 log << MSG::INFO << "NPI : " << m_pass[3] << endreq;
00366 log << MSG::INFO << "Second Vertex Cut : " << m_pass[4] << endreq;
00367 return StatusCode::SUCCESS;
00368 }
00369