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/ITHistSvc.h"
00020 #include "GaudiKernel/Bootstrap.h"
00021
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 "DQAinclKstar/inclkstar.h"
00039
00040 #include <vector>
00041
00042 const double ecms = 3.097;
00043 const double mpi = 0.13957;
00044 const double mk = 0.493677;
00045 const double mkstar0 = 0.896;
00046 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00047
00048 const double velc = 299.792458;
00049
00050 typedef std::vector<int> Vint;
00051 typedef std::vector<HepLorentzVector> Vp4;
00052
00054
00055
00056
00057
00058
00060
00061 inclkstar::inclkstar(const std::string& name, ISvcLocator* pSvcLocator) :
00062 Algorithm(name, pSvcLocator) {
00063
00064 m_tuple2 = 0;
00065 for(int i = 0; i < 10; i++) m_pass[i] = 0;
00066
00067
00068
00069 declareProperty("Vr0cut", m_vr0cut=5.0);
00070 declareProperty("Vz0cut", m_vz0cut=10.0);
00071 declareProperty("CheckDedx", m_checkDedx = 1);
00072 declareProperty("CheckTof", m_checkTof = 1);
00073 }
00074
00075
00076 StatusCode inclkstar::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* inclkstar_mass = new TH1F("InclKstar_mass","INCLUSIVE_Kstar_MASS",65,0.65,1.3);
00089 inclkstar_mass->GetXaxis()->SetTitle("M_{K#pi} (GeV)");
00090 inclkstar_mass->GetYaxis()->SetTitle("Nentries/10MeV/c^{2}");
00091
00092 if(m_thsvc->regHist("/DQAHist/InclKstar/InclKstar_mass", inclkstar_mass).isFailure()) {
00093 log << MSG::ERROR << "Couldn't register InclKstar_mass" << endreq;
00094 }
00095
00096
00097 NTuplePtr nt2(ntupleSvc(), "DQAFILE/InclKstar");
00098 if ( nt2 ) m_tuple2 = nt2;
00099 else {
00100 m_tuple2 = ntupleSvc()->book ("DQAFILE/InclKstar", CLID_ColumnWiseTuple, "inclkstar Ntuple");
00101 if ( m_tuple2 ) {
00102 status = m_tuple2->addItem ("nkaonm", m_nkm);
00103 status = m_tuple2->addItem ("nkaonp", m_nkp);
00104 status = m_tuple2->addItem ("npionp", m_npip);
00105 status = m_tuple2->addItem ("npionm", m_npim);
00106 status = m_tuple2->addItem ("ncharge", m_ncharge);
00107 status = m_tuple2->addItem ("difchikp", m_difchikp);
00108 status = m_tuple2->addItem ("difchikm", m_difchikm);
00109 status = m_tuple2->addItem ("mkstarkp", m_kstarkp);
00110 status = m_tuple2->addItem ("mkstarkm", m_kstarkm);
00111 status = m_tuple2->addItem ("mkstar", m_mkstar);
00112 status = m_tuple2->addItem ("pkstar", m_pkstar);
00113 }
00114 else {
00115 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
00116 return StatusCode::FAILURE;
00117 }
00118 }
00119
00120
00121
00122
00123 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00124 return StatusCode::SUCCESS;
00125
00126 }
00127
00128
00129 StatusCode inclkstar::execute() {
00130 StatusCode sc = StatusCode::SUCCESS;
00131
00132 MsgStream log(msgSvc(), name());
00133 log << MSG::INFO << "in execute()" << endreq;
00134
00135
00136
00137
00138 setFilterPassed(false);
00139
00140 m_pass[0] += 1;
00141
00142 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00143 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00144 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00145
00146 Vint iGood, ikm, ikp, ipip, ipim, iGam;
00147 iGood.clear();
00148 ikm.clear();
00149 ikp.clear();
00150 ipip.clear();
00151 ipim.clear();
00152 iGam.clear();
00153
00154 Vp4 ppip, ppim, pkm, pkp;
00155 ppip.clear();
00156 ppim.clear();
00157 pkm.clear();
00158 pkp.clear();
00159
00160 int TotCharge = 0;
00161 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00162 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00163 if(!(*itTrk)->isMdcTrackValid()) continue;
00164 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00165 if(fabs(mdcTrk->z()) >= m_vz0cut) continue;
00166 if(mdcTrk->r() >= m_vr0cut) continue;
00167 iGood.push_back(i);
00168 TotCharge += mdcTrk->charge();
00169 }
00170
00171
00172
00173 int nGood = iGood.size();
00174
00175
00176
00177
00178
00179 if((nGood < 2) || (TotCharge!=0)) return sc;
00180
00181 m_pass[1] += 1;
00182
00183
00184
00185
00186 ParticleID *pid = ParticleID::instance();
00187 for(int i = 0; i < nGood; i++) {
00188 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00189
00190 pid->init();
00191 pid->setMethod(pid->methodProbability());
00192 pid->setChiMinCut(4);
00193 pid->setRecTrack(*itTrk);
00194 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
00195 pid->identify(pid->onlyPion() | pid->onlyKaon());
00196
00197
00198 pid->calculate();
00199 if(!(pid->IsPidInfoValid())) continue;
00200
00201 if(!(*itTrk)->isMdcKalTrackValid()) continue ;
00202 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
00203 if(pid->probKaon() > 0.001 && (pid->probKaon() > pid->probPion())){
00204 mdcKalTrk->setPidType(RecMdcKalTrack::kaon);
00205 HepLorentzVector ptrk;
00206 ptrk.setPx(mdcKalTrk->px());
00207 ptrk.setPy(mdcKalTrk->py());
00208 ptrk.setPz(mdcKalTrk->pz());
00209 double p3 = ptrk.mag();
00210 ptrk.setE(sqrt(p3*p3+mk*mk));
00211 if(mdcKalTrk->charge() >0) {
00212 ikp.push_back(iGood[i]);
00213 pkp.push_back(ptrk);
00214 } else {
00215 ikm.push_back(iGood[i]);
00216 pkm.push_back(ptrk);
00217 }
00218 }
00219
00220 if(pid->probPion() > 0.001 && (pid->probPion() > pid->probKaon())){
00221 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
00222 HepLorentzVector ptrk;
00223 ptrk.setPx(mdcKalTrk->px());
00224 ptrk.setPy(mdcKalTrk->py());
00225 ptrk.setPz(mdcKalTrk->pz());
00226 double p3 = ptrk.mag();
00227 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00228 if(mdcKalTrk->charge() >0) {
00229 ipip.push_back(iGood[i]);
00230 ppip.push_back(ptrk);
00231 } else {
00232 ipim.push_back(iGood[i]);
00233 ppim.push_back(ptrk);
00234 }
00235 }
00236 }
00237 int npip = ipip.size();
00238 int npim = ipim.size();
00239 int nkm = ikm.size();
00240 int nkp = ikp.size();
00241
00242 m_nkm=nkm;
00243 m_nkp=nkp;
00244 m_npip=npip;
00245 m_npim=npim;
00246 m_ncharge=nGood;
00247 m_pass[3] += 1;
00248
00249 if(npip < 1 && npim < 1) return sc;
00250 if(nkp < 1 && nkm < 1) return sc;
00251
00252 m_pass[4] += 1;
00253
00254
00255
00256
00257
00258 HepLorentzVector pkstar, pkstar1, ppi1, ppi2, pTot;
00259 Vint ikstar;
00260 ikstar.clear();
00261
00262 double difchi0=99999.0;
00263 int ixpim = -1;
00264 int ixkp = -1;
00265
00266 for (int i = 0; i < npim; i++) {
00267 for (int j = 0; j < nkp; j++) {
00268 pkstar = ppim[i] + pkp[j];
00269 double difchi =fabs(pkstar.m()-mkstar0);
00270 if(difchi < difchi0){
00271 difchi0 = difchi;
00272 ixpim = i;
00273 ixkp = j;
00274 }
00275 }
00276 }
00277
00278 m_difchikp = difchi0;
00279
00280 if(ixpim != -1) m_kstarkp=(ppim[ixpim] + pkp[ixkp]).m();
00281
00282 double difchi1=99999.0;
00283 int ixpip = -1;
00284 int ixkm = -1;
00285
00286 for (int ii = 0; ii < npip; ii++) {
00287 for (int jj = 0; jj < nkm; jj++) {
00288 pkstar1 = ppip[ii] + pkm[jj];
00289 double difchi2 =fabs(pkstar1.m()-mkstar0);
00290 if(difchi2 < difchi1){
00291 difchi1 = difchi2;
00292 ixpip = ii;
00293 ixkm = jj;
00294 }
00295 }
00296 }
00297
00298 m_difchikm = difchi1;
00299
00300 if(ixpip != -1) m_kstarkm=(ppip[ixpip] + pkm[ixkm]).m();
00301
00302
00303 if(ixpip == -1 && ixpim == -1) return sc;
00304
00305 if(difchi0 < difchi1){
00306 pTot = ppim[ixpim] + pkp[ixkp];
00307
00308 } else{
00309 pTot = ppip[ixpip] + pkm[ixkm];
00310
00311 }
00312 m_mkstar = pTot.m();
00313 m_pkstar = pTot.rho();
00314
00315
00316 TH1 *h(0);
00317 if (m_thsvc->getHist("/DQAHist/InclKstar/InclKstar_mass", h).isSuccess()) {
00318 h->Fill(m_mkstar);
00319 } else {
00320 log << MSG::ERROR << "Couldn't retrieve InclKstar_mass" << endreq;
00321 }
00322
00323 m_tuple2->write();
00325
00326
00327 if(ixpim != -1){
00328
00329
00330 (*(evtRecTrkCol->begin()+ipim[ixpim]))->tagPion();
00331 (*(evtRecTrkCol->begin()+ikp[ixkp]))->tagKaon();
00332 (*(evtRecTrkCol->begin()+ipim[ixpim]))->setQuality(3);
00333 (*(evtRecTrkCol->begin()+ikp[ixkp]))->setQuality(3);
00334 } else {
00335
00336
00337 (*(evtRecTrkCol->begin()+ipip[ixpip]))->tagPion();
00338 (*(evtRecTrkCol->begin()+ikm[ixkm]))->tagKaon();
00339 (*(evtRecTrkCol->begin()+ipip[ixpip]))->setQuality(3);
00340 (*(evtRecTrkCol->begin()+ikm[ixkm]))->setQuality(3);
00341 }
00342
00343
00344
00345 setFilterPassed(true);
00346
00347
00348 return StatusCode::SUCCESS;
00349 }
00350
00351
00352 StatusCode inclkstar::finalize() {
00353
00354 MsgStream log(msgSvc(), name());
00355 log << MSG::INFO << "in finalize()" << endmsg;
00356 log << MSG::INFO << "Total Entries : " << m_pass[0] << endreq;
00357 log << MSG::INFO << "TOF,dEdx: " << m_pass[1] << endreq;
00358 log << MSG::INFO << "Ncharge Cut : " << m_pass[2] << endreq;
00359 log << MSG::INFO << "PID : " << m_pass[3] << endreq;
00360 log << MSG::INFO << "Npi and Nk Cut : " << m_pass[4] << endreq;
00361 return StatusCode::SUCCESS;
00362 }
00363