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 "VertexFit/IVertexDbSvc.h"
00008 #include "GaudiKernel/Bootstrap.h"
00009
00010 #include "EventModel/EventModel.h"
00011 #include "EventModel/Event.h"
00012 #include "EventModel/EventHeader.h"
00013
00014 #include "EvtRecEvent/EvtRecEvent.h"
00015 #include "EvtRecEvent/EvtRecTrack.h"
00016 #include "DstEvent/TofHitStatus.h"
00017
00018 #include "TMath.h"
00019 #include "CLHEP/Vector/ThreeVector.h"
00020 #include "CLHEP/Vector/LorentzVector.h"
00021 #include "CLHEP/Vector/TwoVector.h"
00022
00023 #include "VertexFit/KinematicFit.h"
00024 #include "VertexFit/VertexFit.h"
00025 #include "VertexFit/SecondVertexFit.h"
00026 #include "VertexFit/IVertexDbSvc.h"
00027 #include "VertexFit/Helix.h"
00028
00029 #include "ParticleID/ParticleID.h"
00030
00031 #include "TMath.h"
00032 #include "TH1F.h"
00033 #include "GaudiKernel/INTupleSvc.h"
00034 #include "GaudiKernel/NTuple.h"
00035 #include "GaudiKernel/Bootstrap.h"
00036 #include "GaudiKernel/ISvcLocator.h"
00037 #include "GaudiKernel/ITHistSvc.h"
00038
00039
00040 using CLHEP::Hep3Vector;
00041 using CLHEP::Hep2Vector;
00042 using CLHEP::HepLorentzVector;
00043 #include "CLHEP/Geometry/Point3D.h"
00044 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00045 typedef HepGeom::Point3D<double> HepPoint3D;
00046 #endif
00047
00048 #include "VertexFit/KinematicFit.h"
00049 #include "VertexFit/VertexFit.h"
00050 #include "VertexFit/SecondVertexFit.h"
00051 #include "VertexFit/IVertexDbSvc.h"
00052 #include "ParticleID/ParticleID.h"
00053 #include <vector>
00054
00055 #include "DQAKsKpiAlg/DQAKsKpi.h"
00056
00057 const double mpi = 0.13957;
00058 const double mk = 0.493677;
00059 const double mks0 = 0.497648;
00060 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00061 const double velc = 299.792458;
00062 typedef std::vector<int> Vint;
00063 typedef std::vector<HepLorentzVector> Vp4;
00064
00065
00066 const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097);
00067 const Hep3Vector u_cms = -p_cms.boostVector();
00068
00069 static int counter[10]={0,0,0,0,0,0,0,0,0,0};
00070
00072
00073 DQAKsKpi::DQAKsKpi(const std::string& name, ISvcLocator* pSvcLocator) :
00074 Algorithm(name, pSvcLocator) {
00075
00076
00077 declareProperty("Vr0cut", m_vr0cut=5.0);
00078 declareProperty("Vz0cut", m_vz0cut=20.0);
00079 declareProperty("Vr1cut", m_vr1cut=1.0);
00080 declareProperty("Vz1cut", m_vz1cut=5.0);
00081 declareProperty("Vctcut", m_cthcut=0.93);
00082 declareProperty("EnergyThreshold", m_energyThreshold=0.04);
00083 declareProperty("GammaAngCut", m_gammaAngCut=20.0);
00084 }
00085
00086
00087
00088 StatusCode DQAKsKpi::initialize(){
00089 MsgStream log(msgSvc(), name());
00090
00091 log << MSG::INFO << "in initialize()" << endmsg;
00092 StatusCode status;
00093
00094 if(service("THistSvc", m_thsvc).isFailure()) {
00095 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
00096 return StatusCode::FAILURE;
00097 }
00098
00099
00100
00101 TH1F* hks_dl = new TH1F("ks_dl", "ks_dl", 300, -5.0, 25.0);
00102 if(m_thsvc->regHist("/DQAHist/DQAKsKpi/hks_dl", hks_dl).isFailure()) {
00103 log << MSG::ERROR << "Couldn't register ks_dl" << endreq;
00104 }
00105
00106 TH1F* hks_m = new TH1F("ks_m", "ks_m", 200,0.4, 0.6);
00107 if(m_thsvc->regHist("/DQAHist/DQAKsKpi/hks_m", hks_m).isFailure()) {
00108 log << MSG::ERROR << "Couldn't register ks_m" << endreq;
00109 }
00110
00111 TH1F* hkspi_m = new TH1F("kspi_m", "kspi_m", 200,0.6, 2.6);
00112 if(m_thsvc->regHist("/DQAHist/DQAKsKpi/hkspi_m", hkspi_m).isFailure()) {
00113 log << MSG::ERROR << "Couldn't register kspi_m" << endreq;
00114 }
00115
00116 TH1F* hks_p = new TH1F("ks_p", "ks_p", 100,0.0, 1.5);
00117 if(m_thsvc->regHist("/DQAHist/DQAKsKpi/hks_p", hks_p).isFailure()) {
00118 log << MSG::ERROR << "Couldn't register ks_p" << endreq;
00119 }
00120
00121 TH1F* hkpi_m = new TH1F("kpi_m", "kpi_m", 200,0.6, 2.6);
00122 if(m_thsvc->regHist("/DQAHist/DQAKsKpi/hkpi_m", hkpi_m).isFailure()) {
00123 log << MSG::ERROR << "Couldn't register kpi_m" << endreq;
00124 }
00125
00126
00127
00128
00129 NTuplePtr nt(ntupleSvc(), "DQAFILE/KsKpi");
00130 if ( nt ) m_tuple = nt;
00131 else {
00132 m_tuple = ntupleSvc()->book("DQAFILE/KsKpi", CLID_ColumnWiseTuple, "KsKpi ntuple");
00133 if( m_tuple ) {
00134 status = m_tuple->addItem("runNo", m_runNo);
00135 status = m_tuple->addItem("event", m_event);
00136
00137
00138
00139 status = m_tuple->addItem("npip", m_npip);
00140 status = m_tuple->addItem("npim", m_npim);
00141 status = m_tuple->addItem("nkp", m_nkp);
00142 status = m_tuple->addItem("nkm", m_nkm);
00143 status = m_tuple->addItem("np", m_np);
00144 status = m_tuple->addItem("npb", m_npb);
00145
00146 status = m_tuple->addItem("vfits_chi", m_vfits_chi);
00147 status = m_tuple->addItem("vfits_vx", m_vfits_vx);
00148 status = m_tuple->addItem("vfits_vy", m_vfits_vy);
00149 status = m_tuple->addItem("vfits_vz", m_vfits_vz);
00150 status = m_tuple->addItem("vfits_vr", m_vfits_vr);
00151
00152 status = m_tuple->addItem("vfitp_chi", m_vfitp_chi);
00153 status = m_tuple->addItem("vfitp_vx", m_vfitp_vx);
00154 status = m_tuple->addItem("vfitp_vy", m_vfitp_vy);
00155 status = m_tuple->addItem("vfitp_vz", m_vfitp_vz);
00156 status = m_tuple->addItem("vfitp_vr", m_vfitp_vr);
00157
00158 status = m_tuple->addItem("vfit2_chi", m_vfit2_chi);
00159 status = m_tuple->addItem("vfit2_mks", m_vfit2_mks);
00160 status = m_tuple->addItem("vfit2_ct", m_vfit2_ct);
00161 status = m_tuple->addItem("vfit2_dl", m_vfit2_dl);
00162 status = m_tuple->addItem("vfit2_dle", m_vfit2_dle);
00163
00164 status = m_tuple->addItem("chi2_fs4c", m_chi2_fs4c);
00165 status = m_tuple->addItem("mks_fs4c", m_mks_fs4c);
00166 status = m_tuple->addItem("mkspi_fs4c",m_mkspi_fs4c);
00167 status = m_tuple->addItem("mksk_fs4c", m_mksk_fs4c);
00168 status = m_tuple->addItem("mkpi_fs4c", m_mkpi_fs4c);
00169
00170 status = m_tuple->addItem("4c_chi2", m_4c_chi2);
00171 status = m_tuple->addItem("4c_mks", m_4c_mks);
00172 status = m_tuple->addItem("4c_mkspi", m_4c_mkspi);
00173 status = m_tuple->addItem("4c_mksk", m_4c_mksk);
00174 status = m_tuple->addItem("4c_mkpi", m_4c_mkpi);
00175 status = m_tuple->addItem("4c_ks_px", m_4c_ks_px);
00176 status = m_tuple->addItem("4c_ks_py", m_4c_ks_py);
00177 status = m_tuple->addItem("4c_ks_pz", m_4c_ks_pz);
00178 status = m_tuple->addItem("4c_ks_p", m_4c_ks_p);
00179 status = m_tuple->addItem("4c_ks_cos", m_4c_ks_cos);
00180
00181 status = m_tuple->addItem("NGch", m_ngch, 0, 10);
00182 status = m_tuple->addIndexedItem("pidcode", m_ngch, m_pidcode);
00183 status = m_tuple->addIndexedItem("pidprob", m_ngch, m_pidprob);
00184 status = m_tuple->addIndexedItem("pidchiDedx", m_ngch, m_pidchiDedx);
00185 status = m_tuple->addIndexedItem("pidchiTof1", m_ngch, m_pidchiTof1);
00186 status = m_tuple->addIndexedItem("pidchiTof2", m_ngch, m_pidchiTof2);
00187
00188 status = m_tuple->addIndexedItem("charge",m_ngch, m_charge);
00189 status = m_tuple->addIndexedItem("vx0", m_ngch, m_vx0);
00190 status = m_tuple->addIndexedItem("vy0", m_ngch, m_vy0);
00191 status = m_tuple->addIndexedItem("vz0", m_ngch, m_vz0);
00192 status = m_tuple->addIndexedItem("vr0", m_ngch, m_vr0);
00193
00194 status = m_tuple->addIndexedItem("vx", m_ngch, m_vx);
00195 status = m_tuple->addIndexedItem("vy", m_ngch, m_vy);
00196 status = m_tuple->addIndexedItem("vz", m_ngch, m_vz);
00197 status = m_tuple->addIndexedItem("vr", m_ngch, m_vr);
00198
00199 status = m_tuple->addIndexedItem("px", m_ngch, m_px);
00200 status = m_tuple->addIndexedItem("py", m_ngch, m_py);
00201 status = m_tuple->addIndexedItem("pz", m_ngch, m_pz);
00202 status = m_tuple->addIndexedItem("p", m_ngch, m_p);
00203 status = m_tuple->addIndexedItem("cost", m_ngch, m_cost);
00204
00205 status = m_tuple->addIndexedItem("probPH", m_ngch, m_probPH);
00206 status = m_tuple->addIndexedItem("normPH", m_ngch, m_normPH);
00207 status = m_tuple->addIndexedItem("chie", m_ngch, m_chie);
00208 status = m_tuple->addIndexedItem("chimu", m_ngch, m_chimu);
00209 status = m_tuple->addIndexedItem("chipi", m_ngch, m_chipi);
00210 status = m_tuple->addIndexedItem("chik", m_ngch, m_chik);
00211 status = m_tuple->addIndexedItem("chip", m_ngch, m_chip);
00212 status = m_tuple->addIndexedItem("ghit", m_ngch, m_ghit);
00213 status = m_tuple->addIndexedItem("thit", m_ngch, m_thit);
00214
00215 status = m_tuple->addIndexedItem("e_emc", m_ngch, m_e_emc);
00216
00217 status = m_tuple->addIndexedItem("qual_etof", m_ngch, m_qual_etof);
00218 status = m_tuple->addIndexedItem("tof_etof", m_ngch, m_tof_etof);
00219 status = m_tuple->addIndexedItem("te_etof", m_ngch, m_te_etof);
00220 status = m_tuple->addIndexedItem("tmu_etof", m_ngch, m_tmu_etof);
00221 status = m_tuple->addIndexedItem("tpi_etof", m_ngch, m_tpi_etof);
00222 status = m_tuple->addIndexedItem("tk_etof", m_ngch, m_tk_etof);
00223 status = m_tuple->addIndexedItem("tp_etof", m_ngch, m_tp_etof);
00224
00225 status = m_tuple->addIndexedItem("qual_btof1", m_ngch, m_qual_btof1);
00226 status = m_tuple->addIndexedItem("tof_btof1", m_ngch, m_tof_btof1);
00227 status = m_tuple->addIndexedItem("te_btof1", m_ngch, m_te_btof1);
00228 status = m_tuple->addIndexedItem("tmu_btof1", m_ngch, m_tmu_btof1);
00229 status = m_tuple->addIndexedItem("tpi_btof1", m_ngch, m_tpi_btof1);
00230 status = m_tuple->addIndexedItem("tk_btof1", m_ngch, m_tk_btof1);
00231 status = m_tuple->addIndexedItem("tp_btof1", m_ngch, m_tp_btof1);
00232
00233 status = m_tuple->addIndexedItem("qual_btof2", m_ngch, m_qual_btof2);
00234 status = m_tuple->addIndexedItem("tof_btof2", m_ngch, m_tof_btof2);
00235 status = m_tuple->addIndexedItem("te_btof2", m_ngch, m_te_btof2);
00236 status = m_tuple->addIndexedItem("tmu_btof2", m_ngch, m_tmu_btof2);
00237 status = m_tuple->addIndexedItem("tpi_btof2", m_ngch, m_tpi_btof2);
00238 status = m_tuple->addIndexedItem("tk_btof2", m_ngch, m_tk_btof2);
00239 status = m_tuple->addIndexedItem("tp_btof2", m_ngch, m_tp_btof2);
00240
00241 } else {
00242 log << MSG::ERROR << "Can not book N-tuple:" << long(m_tuple) << endreq;
00243 }
00244 }
00245
00246
00247
00248 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00249 return StatusCode::SUCCESS;
00250
00251 }
00252
00253
00254 StatusCode DQAKsKpi::execute() {
00255
00256 MsgStream log(msgSvc(), name());
00257 log << MSG::INFO << "in execute()" << endreq;
00258
00259
00260
00261
00262
00263 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00264 m_runNo = eventHeader->runNumber();
00265 m_event=eventHeader->eventNumber();
00266 log << MSG::DEBUG <<"run, evtnum = "
00267 << m_runNo << " , "
00268 << m_event <<endreq;
00269
00270 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00271
00272
00273 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00274 << evtRecEvent->totalCharged() << " , "
00275 << evtRecEvent->totalNeutral() << " , "
00276 << evtRecEvent->totalTracks() <<endreq;
00277
00278 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00279
00280 if(evtRecEvent->totalNeutral()>100) {
00281 return StatusCode::SUCCESS;
00282 }
00283
00284 Vint iGood;
00285 iGood.clear();
00286
00287 Hep3Vector xorigin(0,0,0);
00288
00289 IVertexDbSvc* vtxsvc;
00290 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00291 if(vtxsvc->isVertexValid()){
00292 double* dbv = vtxsvc->PrimaryVertex();
00293 double* vv = vtxsvc->SigmaPrimaryVertex();
00294 xorigin.setX(dbv[0]);
00295 xorigin.setY(dbv[1]);
00296 xorigin.setZ(dbv[2]);
00297 log << MSG::INFO
00298 <<"xorigin.x="<<xorigin.x()<<", "
00299 <<"xorigin.y="<<xorigin.y()<<", "
00300 <<"xorigin.z="<<xorigin.z()<<", "
00301 <<endreq ;
00302 }
00303
00304 int nCharge = 0;
00305 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00306 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00307 if(!(*itTrk)->isMdcTrackValid()) continue;
00308 if (!(*itTrk)->isMdcKalTrackValid()) continue;
00309 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00310
00311 double pch =mdcTrk->p();
00312 double x0 =mdcTrk->x();
00313 double y0 =mdcTrk->y();
00314 double z0 =mdcTrk->z();
00315 double phi0=mdcTrk->helix(1);
00316 double xv=xorigin.x();
00317 double yv=xorigin.y();
00318 double Rxy=fabs((x0-xv)*cos(phi0)+(y0-yv)*sin(phi0));
00319 double vx0 = x0;
00320 double vy0 = y0;
00321 double vz0 = z0-xorigin.z();
00322 double vr0 = Rxy;
00323 double Vct=cos(mdcTrk->theta());
00324
00325 HepVector a = mdcTrk->helix();
00326 HepSymMatrix Ea = mdcTrk->err();
00327 HepPoint3D point0(0.,0.,0.);
00328 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
00329 VFHelix helixip(point0,a,Ea);
00330 helixip.pivot(IP);
00331 HepVector vecipa = helixip.a();
00332 double Rvxy0=fabs(vecipa[0]);
00333 double Rvz0=vecipa[3];
00334 double Rvphi0=vecipa[1];
00335
00336
00337
00338
00339 if(fabs(Rvz0) >= m_vz0cut) continue;
00340 if(Rvxy0 >= m_vr0cut) continue;
00341
00342
00343 iGood.push_back(i);
00344 nCharge += mdcTrk->charge();
00345 }
00346
00347
00348
00349
00350 int m_ngch = iGood.size();
00351 log << MSG::DEBUG << "ngood, totcharge = " << m_ngch << " , " << nCharge << endreq;
00352 if((m_ngch != 4)||(nCharge!=0)){
00353 return StatusCode::SUCCESS;
00354 }
00355
00356
00357
00358 Vint ipip, ipim, ikp, ikm, ipp, ipm;
00359 ipip.clear();
00360 ipim.clear();
00361 ikp.clear();
00362 ikm.clear();
00363 ipp.clear();
00364 ipm.clear();
00365
00366 Vp4 p_pip, p_pim, p_kp, p_km, p_pp, p_pm ;
00367 p_pip.clear();
00368 p_pim.clear();
00369 p_kp.clear();
00370 p_km.clear();
00371 p_pp.clear();
00372 p_pm.clear();
00373
00374 ParticleID *pid = ParticleID::instance();
00375 for(int i = 0; i < m_ngch; i++) {
00376 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00377
00378 pid->init();
00379 pid->setMethod(pid->methodProbability());
00380 pid->setChiMinCut(4);
00381 pid->setRecTrack(*itTrk);
00382 pid->usePidSys(pid->useDedx());
00383
00384 pid->identify(pid->onlyPionKaonProton());
00385
00386
00387
00388 pid->calculate();
00389 if(!(pid->IsPidInfoValid())) continue;
00390
00391 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00392 RecMdcKalTrack* mdcKalTrk = 0 ;
00393 if((*itTrk)->isMdcKalTrackValid()) mdcKalTrk = (*itTrk)->mdcKalTrack();
00394
00395 double prob_pi = pid->probPion();
00396 double prob_K = pid->probKaon();
00397 double prob_p = pid->probProton();
00398
00399 HepLorentzVector ptrk;
00400 ptrk.setPx(mdcTrk->px()) ;
00401 ptrk.setPy(mdcTrk->py()) ;
00402 ptrk.setPz(mdcTrk->pz()) ;
00403 double p3 = ptrk.mag() ;
00404
00405 if (prob_pi > prob_K && prob_pi > prob_p) {
00406 m_pidcode[i]=2;
00407 m_pidprob[i]=pid->prob(2);
00408 m_pidchiDedx[i]=pid->chiDedx(2);
00409 m_pidchiTof1[i]=pid->chiTof1(2);
00410 m_pidchiTof2[i]=pid->chiTof2(2);
00411 ptrk.setE(sqrt(p3*p3 + xmass[2]*xmass[2])) ;
00412 if(mdcTrk->charge() > 0) {
00413 ipip.push_back(iGood[i]);
00414 p_pip.push_back(ptrk);
00415 }
00416 if (mdcTrk->charge() < 0) {
00417 ipim.push_back(iGood[i]);
00418 p_pim.push_back(ptrk);
00419 }
00420 }
00421
00422 if (prob_K > prob_pi && prob_K > prob_p) {
00423 m_pidcode[i]=3;
00424 m_pidprob[i]=pid->prob(3);
00425 m_pidchiDedx[i]=pid->chiDedx(3);
00426 m_pidchiTof1[i]=pid->chiTof1(3);
00427 m_pidchiTof2[i]=pid->chiTof2(3);
00428 ptrk.setE(sqrt(p3*p3 + xmass[3]*xmass[3])) ;
00429 if(mdcTrk->charge() > 0) {
00430 ikp.push_back(iGood[i]);
00431 p_kp.push_back(ptrk);
00432 }
00433 if (mdcTrk->charge() < 0) {
00434 ikm.push_back(iGood[i]);
00435 p_km.push_back(ptrk);
00436 }
00437 }
00438
00439 if (prob_p > prob_pi && prob_p > prob_K) {
00440 m_pidcode[i]=4;
00441 m_pidprob[i]=pid->prob(4);
00442 m_pidchiDedx[i]=pid->chiDedx(4);
00443 m_pidchiTof1[i]=pid->chiTof1(4);
00444 m_pidchiTof2[i]=pid->chiTof2(4);
00445 ptrk.setE(sqrt(p3*p3 + xmass[4]*xmass[4])) ;
00446 if(mdcTrk->charge() > 0) {
00447 ipp.push_back(iGood[i]);
00448 p_pp.push_back(ptrk);
00449 }
00450 if (mdcTrk->charge() < 0) {
00451 ipm.push_back(iGood[i]);
00452 p_pm.push_back(ptrk);
00453 }
00454 }
00455 }
00456
00457 m_npip= ipip.size() ;
00458 m_npim= ipim.size() ;
00459 m_nkp = ikp.size() ;
00460 m_nkm = ikm.size() ;
00461 m_np = ipp.size() ;
00462 m_npb = ipm.size() ;
00463
00464
00465 if( m_npip*m_npim != 2 ) {
00466 return StatusCode::SUCCESS;
00467 }
00468
00469 if( m_nkp+m_nkm != 1 ) {
00470 return StatusCode::SUCCESS;
00471 }
00472
00473
00474
00475 HepPoint3D vx(0., 0., 0.);
00476 HepSymMatrix Evx(3, 0);
00477 double bx = 1E+6;
00478 double by = 1E+6;
00479 double bz = 1E+6;
00480 Evx[0][0] = bx*bx;
00481 Evx[1][1] = by*by;
00482 Evx[2][2] = bz*bz;
00483
00484 VertexParameter vxpar;
00485 vxpar.setVx(vx);
00486 vxpar.setEvx(Evx);
00487
00488 VertexFit *vtxfit_s = VertexFit::instance();
00489 VertexFit *vtxfit_p = VertexFit::instance();
00490 SecondVertexFit *vtxfit2 = SecondVertexFit::instance();
00491
00492 RecMdcKalTrack *pi1KalTrk, *pi2KalTrk, *pi3KalTrk, *k1KalTrk;
00493 RecMdcKalTrack *pipKalTrk, *pimKalTrk, *piKalTrk, *kKalTrk;
00494 WTrackParameter wks,vwks;
00495
00496 double chi_temp = 999.0;
00497 double mks_temp = 10.0 ;
00498 bool okloop=false;
00499 for(unsigned int i1 = 0; i1 < m_npip; i1++) {
00500 RecMdcKalTrack *pi1KalTrk = (*(evtRecTrkCol->begin()+ipip[i1]))-> mdcKalTrack();
00501 pi1KalTrk->setPidType(RecMdcKalTrack::pion);
00502 WTrackParameter wpi1trk(xmass[2], pi1KalTrk->getZHelix(), pi1KalTrk->getZError());
00503
00504 for(unsigned int i2 = 0; i2 < m_npim; i2++) {
00505 RecMdcKalTrack *pi2KalTrk = (*(evtRecTrkCol->begin()+ipim[i2]))-> mdcKalTrack();
00506 pi2KalTrk->setPidType(RecMdcKalTrack::pion);
00507 WTrackParameter wpi2trk(xmass[2], pi2KalTrk->getZHelix(), pi2KalTrk->getZError());
00508
00509 vtxfit_s->init();
00510 vtxfit_s->AddTrack(0, wpi1trk);
00511 vtxfit_s->AddTrack(1, wpi2trk);
00512 vtxfit_s->AddVertex(0, vxpar, 0, 1);
00513
00514 if(!(vtxfit_s->Fit(0))) continue;
00515 vtxfit_s->BuildVirtualParticle(0);
00516 m_vfits_chi = vtxfit_s->chisq(0);
00517 WTrackParameter wkshort = vtxfit_s->wVirtualTrack(0);
00518 VertexParameter vparks = vtxfit_s->vpar(0);
00519
00520 m_vfits_vx = (vparks.Vx())[0];
00521 m_vfits_vy = (vparks.Vx())[1];
00522 m_vfits_vz = (vparks.Vx())[2];
00523 m_vfits_vr = sqrt(m_vfits_vx*m_vfits_vx + m_vfits_vy*m_vfits_vy) ;
00524
00525 if ( m_npip == 2 ) {
00526 int j = i1 ;
00527 int jj = ( i1 == 1 ) ? 0 : 1;
00528 pi3KalTrk = (*(evtRecTrkCol->begin()+ipip[jj]))->mdcKalTrack();
00529 k1KalTrk = (*(evtRecTrkCol->begin()+ikm[0]))->mdcKalTrack();
00530 }
00531 if (m_npim == 2 ) {
00532 int j = i2 ;
00533 int jj = ( i2 == 1 ) ? 0 : 1;
00534 pi3KalTrk = (*(evtRecTrkCol->begin()+ipim[jj]))->mdcKalTrack();
00535 k1KalTrk = (*(evtRecTrkCol->begin()+ikp[0]))->mdcKalTrack();
00536 }
00537
00538 pi3KalTrk->setPidType(RecMdcKalTrack::pion);
00539 WTrackParameter wpi3trk( xmass[2], pi3KalTrk->getZHelix(), pi3KalTrk->getZError());
00540 k1KalTrk->setPidType(RecMdcKalTrack::kaon);
00541 WTrackParameter wk1trk( xmass[3], k1KalTrk->getZHelixK(), k1KalTrk->getZErrorK());
00542
00543 vtxfit_p->init();
00544
00545 vtxfit_p->AddTrack(0, wpi3trk);
00546 vtxfit_p->AddTrack(1, wk1trk);
00547 vtxfit_p->AddVertex(0, vxpar, 0, 1);
00548 if(!(vtxfit_p->Fit(0))) continue;
00549
00550 m_vfitp_chi = vtxfit_p->chisq(0) ;
00551
00552 VertexParameter primaryVpar = vtxfit_p->vpar(0);
00553 m_vfitp_vx = (primaryVpar.Vx())[0];
00554 m_vfitp_vy = (primaryVpar.Vx())[1];
00555 m_vfitp_vz = (primaryVpar.Vx())[2];
00556 m_vfitp_vr = sqrt(m_vfitp_vx*m_vfitp_vx + m_vfitp_vy*m_vfitp_vy);
00557
00558 vtxfit2->init();
00559 vtxfit2->setPrimaryVertex(vtxfit_p->vpar(0));
00560 vtxfit2->AddTrack(0, wkshort);
00561 vtxfit2->setVpar(vtxfit_s->vpar(0));
00562 if(!vtxfit2->Fit()) continue;
00563
00564 if ( fabs(((vtxfit2->wpar()).p()).m()-mks0) < mks_temp ) {
00565 mks_temp = fabs(((vtxfit2->wpar()).p()).m()-mks0) ;
00566
00567 okloop = true;
00568
00569 wks = vtxfit2->wpar();
00570 m_vfit2_mks = (wks.p()).m();
00571 m_vfit2_chi = vtxfit2->chisq();
00572 m_vfit2_ct = vtxfit2->ctau();
00573 m_vfit2_dl = vtxfit2->decayLength();
00574 m_vfit2_dle = vtxfit2->decayLengthError();
00575
00576 pipKalTrk = pi1KalTrk ;
00577 pimKalTrk = pi2KalTrk ;
00578 piKalTrk = pi3KalTrk ;
00579 kKalTrk = k1KalTrk ;
00580
00581 }
00582 }
00583 }
00584
00585 if (! okloop ) {
00586 return StatusCode::SUCCESS;
00587 }
00588
00589 pipKalTrk->setPidType(RecMdcKalTrack::pion);
00590 pimKalTrk->setPidType(RecMdcKalTrack::pion);
00591 piKalTrk->setPidType(RecMdcKalTrack::pion);
00592 kKalTrk->setPidType(RecMdcKalTrack::kaon);
00593
00594 WTrackParameter wpip(xmass[2], pipKalTrk->getZHelix(), pipKalTrk->getZError());
00595 WTrackParameter wpim(xmass[2], pimKalTrk->getZHelix(), pimKalTrk->getZError());
00596
00597 WTrackParameter wpi(xmass[2], piKalTrk->getZHelix(), piKalTrk->getZError());
00598 WTrackParameter wk(xmass[3], kKalTrk->getZHelixK(), kKalTrk->getZErrorK());
00599
00600
00601
00602 int ii ;
00603 for(int j=0; j<m_ngch; j++){
00604 m_charge[j] = 9999.0;
00605 m_vx0[j] = 9999.0;
00606 m_vy0[j] = 9999.0;
00607 m_vz0[j] = 9999.0;
00608 m_vr0[j] = 9999.0;
00609
00610 m_vx[j] = 9999.0;
00611 m_vy[j] = 9999.0;
00612 m_vz[j] = 9999.0;
00613 m_vr[j] = 9999.0;
00614
00615 m_px[j] = 9999.0;
00616 m_py[j] = 9999.0;
00617 m_pz[j] = 9999.0;
00618 m_p[j] = 9999.0;
00619 m_cost[j] = 9999.0;
00620
00621 m_probPH[j] = 9999.0;
00622 m_normPH[j] = 9999.0;
00623 m_chie[j] = 9999.0;
00624 m_chimu[j] = 9999.0;
00625 m_chipi[j] = 9999.0;
00626 m_chik[j] = 9999.0;
00627 m_chip[j] = 9999.0;
00628 m_ghit[j] = 9999.0;
00629 m_thit[j] = 9999.0;
00630
00631 m_e_emc[j] = 9999.0;
00632
00633 m_qual_etof[j] = 9999.0;
00634 m_tof_etof[j] = 9999.0;
00635 m_te_etof[j] = 9999.0;
00636 m_tmu_etof[j] = 9999.0;
00637 m_tpi_etof[j] = 9999.0;
00638 m_tk_etof[j] = 9999.0;
00639 m_tp_etof[j] = 9999.0;
00640
00641 m_qual_btof1[j] = 9999.0;
00642 m_tof_btof1[j] = 9999.0;
00643 m_te_btof1[j] = 9999.0;
00644 m_tmu_btof1[j] = 9999.0;
00645 m_tpi_btof1[j] = 9999.0;
00646 m_tk_btof1[j] = 9999.0;
00647 m_tp_btof1[j] = 9999.0;
00648
00649 m_qual_btof2[j] = 9999.0;
00650 m_tof_btof2[j] = 9999.0;
00651 m_te_btof2[j] = 9999.0;
00652 m_tmu_btof2[j] = 9999.0;
00653 m_tpi_btof2[j] = 9999.0;
00654 m_tk_btof2[j] = 9999.0;
00655 m_tp_btof2[j] = 9999.0;
00656
00657 m_pidcode[j] = 9999.0;
00658 m_pidprob[j] = 9999.0;
00659 m_pidchiDedx[j] = 9999.0;
00660 m_pidchiTof1[j] = 9999.0;
00661 m_pidchiTof2[j] = 99999.0;
00662 }
00663
00664 for(int i = 0; i < m_ngch; i++ ){
00665
00666 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00667 if(!(*itTrk)->isMdcTrackValid()) continue;
00668 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00669 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
00670
00671 if ( mdcKalTrk == pipKalTrk ) {
00672 ii = 0 ;
00673 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
00674 }
00675 if ( mdcKalTrk == pimKalTrk ) {
00676 ii = 1 ;
00677 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
00678 }
00679 if ( mdcKalTrk == piKalTrk ) {
00680 ii = 2 ;
00681 mdcKalTrk->setPidType(RecMdcKalTrack::pion);
00682 }
00683 if ( mdcKalTrk == kKalTrk ) {
00684 ii = 3 ;
00685 mdcKalTrk->setPidType(RecMdcKalTrack::kaon);
00686 }
00687
00688 m_charge[ii] = mdcTrk->charge();
00689 double x0=mdcTrk->x();
00690 double y0=mdcTrk->y();
00691 double z0=mdcTrk->z();
00692 double phi0=mdcTrk->helix(1);
00693 double xv=xorigin.x();
00694 double yv=xorigin.y();
00695 double zv=xorigin.z();
00696 double rv=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
00697
00698 m_vx0[ii] = x0-xv ;
00699 m_vy0[ii] = y0-yv ;
00700 m_vz0[ii] = z0-zv ;
00701 m_vr0[ii] = rv ;
00702
00703 x0=mdcKalTrk->x();
00704 y0=mdcKalTrk->y();
00705 z0=mdcKalTrk->z();
00706 rv=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
00707
00708 m_vx[ii] = x0-xv ;
00709 m_vy[ii] = y0-yv ;
00710 m_vz[ii] = z0-zv ;
00711 m_vr[ii] = rv ;
00712
00713 m_px[ii] = mdcKalTrk->px();
00714 m_py[ii] = mdcKalTrk->py();
00715 m_pz[ii] = mdcKalTrk->pz();
00716 m_p[ii] = mdcKalTrk->p();
00717 m_cost[ii] = mdcKalTrk->pz()/mdcKalTrk->p();
00718
00719 double ptrk = mdcKalTrk->p() ;
00720 if((*itTrk)->isMdcDedxValid()) {
00721 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
00722 m_probPH[ii]= dedxTrk->probPH();
00723 m_normPH[ii]= dedxTrk->normPH();
00724
00725 m_chie[ii] = dedxTrk->chiE();
00726 m_chimu[ii] = dedxTrk->chiMu();
00727 m_chipi[ii] = dedxTrk->chiPi();
00728 m_chik[ii] = dedxTrk->chiK();
00729 m_chip[ii] = dedxTrk->chiP();
00730 m_ghit[ii] = dedxTrk->numGoodHits();
00731 m_thit[ii] = dedxTrk->numTotalHits();
00732 }
00733
00734 if((*itTrk)->isEmcShowerValid()) {
00735 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00736 m_e_emc[ii] = emcTrk->energy();
00737 }
00738
00739 if((*itTrk)->isTofTrackValid()) {
00740 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
00741
00742 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
00743
00744 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
00745 TofHitStatus *status = new TofHitStatus;
00746 status->setStatus((*iter_tof)->status());
00747
00748 if(!(status->is_barrel())){
00749 if( !(status->is_counter()) ) continue;
00750 if( status->layer()!=0 ) continue;
00751 double path=(*iter_tof)->path();
00752 double tof = (*iter_tof)->tof();
00753 double ph = (*iter_tof)->ph();
00754 double rhit = (*iter_tof)->zrhit();
00755 double qual = 0.0 + (*iter_tof)->quality();
00756 double cntr = 0.0 + (*iter_tof)->tofID();
00757 double texp[5];
00758 for(int j = 0; j < 5; j++) {
00759 double gb = ptrk/xmass[j];
00760 double beta = gb/sqrt(1+gb*gb);
00761 texp[j] = path /beta/velc;
00762 }
00763
00764 m_qual_etof[ii] = qual;
00765 m_tof_etof[ii] = tof ;
00766 }
00767 else {
00768 if( !(status->is_counter()) ) continue;
00769 if(status->layer()==1){
00770 double path=(*iter_tof)->path();
00771 double tof = (*iter_tof)->tof();
00772 double ph = (*iter_tof)->ph();
00773 double rhit = (*iter_tof)->zrhit();
00774 double qual = 0.0 + (*iter_tof)->quality();
00775 double cntr = 0.0 + (*iter_tof)->tofID();
00776 double texp[5];
00777 for(int j = 0; j < 5; j++) {
00778 double gb = ptrk/xmass[j];
00779 double beta = gb/sqrt(1+gb*gb);
00780 texp[j] = path /beta/velc;
00781 }
00782
00783 m_qual_btof1[ii] = qual;
00784 m_tof_btof1[ii] = tof ;
00785 }
00786
00787 if(status->layer()==2){
00788 double path=(*iter_tof)->path();
00789 double tof = (*iter_tof)->tof();
00790 double ph = (*iter_tof)->ph();
00791 double rhit = (*iter_tof)->zrhit();
00792 double qual = 0.0 + (*iter_tof)->quality();
00793 double cntr = 0.0 + (*iter_tof)->tofID();
00794 double texp[5];
00795 for(int j = 0; j < 5; j++) {
00796 double gb = ptrk/xmass[j];
00797 double beta = gb/sqrt(1+gb*gb);
00798 texp[j] = path /beta/velc;
00799 }
00800
00801 m_qual_btof2[ii] = qual;
00802 m_tof_btof2[ii] = tof ;
00803 }
00804 }
00805 }
00806 }
00807 }
00808
00809
00810 KinematicFit * kmfit = KinematicFit::instance();
00811
00812 double ecms = 3.097;
00813 double chisq = 9999.;
00814 m_4c_chi2 = 9999.;
00815 m_4c_mks = 10.0;
00816 m_4c_mkspi = 10.0;
00817 m_4c_mksk = 10.0;
00818 m_4c_mkpi = 10.0;
00819 m_4c_ks_px = 10.0;
00820 m_4c_ks_py = 10.0;
00821 m_4c_ks_pz = 10.0;
00822 m_4c_ks_p = 10.0;
00823 m_4c_ks_cos= 10.0;
00824
00825 kmfit->init();
00826 kmfit->AddTrack(0, wpi);
00827 kmfit->AddTrack(1, wk);
00828 kmfit->AddTrack(2, wks);
00829 kmfit->AddFourMomentum(0, p_cms);
00830 bool oksq = kmfit->Fit();
00831 if(oksq) {
00832 chisq = kmfit->chisq();
00833
00834 HepLorentzVector pk = kmfit->pfit(1);
00835 HepLorentzVector pks = kmfit->pfit(2);
00836 HepLorentzVector pkspi = kmfit->pfit(0) + kmfit->pfit(2);
00837 HepLorentzVector pksk = kmfit->pfit(1) + kmfit->pfit(2);
00838 HepLorentzVector pkpi = kmfit->pfit(0) + kmfit->pfit(1);
00839
00840 pks.boost(u_cms);
00841 pkspi.boost(u_cms);
00842 pksk.boost(u_cms);
00843 pkpi.boost(u_cms);
00844
00845 m_4c_chi2 = chisq ;
00846 m_4c_mks = pks.m();
00847 m_4c_mkspi = pkspi.m();
00848 m_4c_mksk = pksk.m();
00849 m_4c_mkpi = pkpi.m();
00850
00851 m_4c_ks_px = pks.px() ;
00852 m_4c_ks_py = pks.py() ;
00853 m_4c_ks_pz = pks.pz() ;
00854 m_4c_ks_p = (pks.vect()).mag() ;
00855 m_4c_ks_cos = m_4c_ks_pz/m_4c_ks_p ;
00856
00857 }
00858
00859 chisq = 9999.;
00860 m_chi2_fs4c = 9999.;
00861 m_mks_fs4c = 10.0;
00862 m_mkspi_fs4c = 10.0;
00863 m_mksk_fs4c = 10.0;
00864 m_mkpi_fs4c = 10.0;
00865
00866 kmfit->init();
00867 kmfit->AddTrack(0, wpip);
00868 kmfit->AddTrack(1, wpim);
00869 kmfit->AddTrack(2, wpi);
00870 kmfit->AddTrack(3, wk);
00871 kmfit->AddFourMomentum(0, p_cms);
00872 oksq = kmfit->Fit();
00873 if(oksq) {
00874 chisq = kmfit->chisq();
00875
00876 HepLorentzVector pks = kmfit->pfit(0) + kmfit->pfit(1);
00877 HepLorentzVector pkspi = pks + kmfit->pfit(2);
00878 HepLorentzVector pksk = pks + kmfit->pfit(3);
00879 HepLorentzVector pkpi = kmfit->pfit(2) + kmfit->pfit(3);
00880
00881 pks.boost(u_cms);
00882 pkspi.boost(u_cms);
00883 pksk.boost(u_cms);
00884 pkpi.boost(u_cms);
00885
00886 m_chi2_fs4c = chisq ;
00887 m_mks_fs4c = pks.m();
00888 m_mkspi_fs4c = pkspi.m();
00889 m_mksk_fs4c = pksk.m();
00890 m_mkpi_fs4c = pkpi.m();
00891 }
00892
00893
00894 if(chisq > 20) { return StatusCode::SUCCESS; }
00895 if(m_vfit2_dl < 0.5) { return StatusCode::SUCCESS; }
00896 if(fabs(m_4c_mks-mks0) > 0.01) { return StatusCode::SUCCESS; }
00897 if(m_4c_mkspi < 1.25) { return StatusCode::SUCCESS; }
00898
00899
00900 TH1 *h(0);
00901 if (m_thsvc->getHist("/DQAHist/DQAKsKpi/hks_dl", h).isSuccess()) {
00902 h->Fill(m_vfit2_dl);
00903 } else {
00904 log << MSG::ERROR << "Couldn't retrieve hks_dl" << endreq;
00905 }
00906
00907 if (m_thsvc->getHist("/DQAHist/DQAKsKpi/hks_m", h).isSuccess()) {
00908 h->Fill(m_4c_mks);
00909 } else {
00910 log << MSG::ERROR << "Couldn't retrieve hks_m" << endreq;
00911 }
00912
00913 if (m_thsvc->getHist("/DQAHist/DQAKsKpi/hkspi_m", h).isSuccess()) {
00914 h->Fill(m_4c_mkspi);
00915 } else {
00916 log << MSG::ERROR << "Couldn't retrieve hkspi_m" << endreq;
00917 }
00918
00919 if (m_thsvc->getHist("/DQAHist/DQAKsKpi/hks_p", h).isSuccess()) {
00920 h->Fill(m_4c_ks_p);
00921 } else {
00922 log << MSG::ERROR << "Couldn't retrieve hks_p" << endreq;
00923 }
00924
00925 if (m_thsvc->getHist("/DQAHist/DQAKsKpi/hkpi_m", h).isSuccess()) {
00926 h->Fill(m_4c_mkpi);
00927 } else {
00928 log << MSG::ERROR << "Couldn't retrieve hkpi_m" << endreq;
00929 }
00930
00932
00933
00934
00935
00936 if(m_npip==2 && m_npim==1){
00937 (*(evtRecTrkCol->begin()+ipip[0]))->setPartId(2);
00938 (*(evtRecTrkCol->begin()+ipip[1]))->setPartId(2);
00939 (*(evtRecTrkCol->begin()+ipim[0]))->setPartId(2);
00940 (*(evtRecTrkCol->begin()+ikm[0]))->setPartId(4);
00941 }
00942 if(m_npip==1 && m_npim==2){
00943 (*(evtRecTrkCol->begin()+ipip[0]))->setPartId(2);
00944 (*(evtRecTrkCol->begin()+ipim[0]))->setPartId(2);
00945 (*(evtRecTrkCol->begin()+ipim[1]))->setPartId(2);
00946 (*(evtRecTrkCol->begin()+ikp[0]))->setPartId(4);
00947 }
00948
00949
00950
00951
00952
00953 if(m_npip==2 && m_npim==1){
00954 (*(evtRecTrkCol->begin()+ipip[0]))->setQuality(1);
00955 (*(evtRecTrkCol->begin()+ipip[1]))->setQuality(1);
00956 (*(evtRecTrkCol->begin()+ipim[0]))->setQuality(1);
00957 (*(evtRecTrkCol->begin()+ikm[0]))->setQuality(1);
00958 }
00959 if(m_npip==1 && m_npim==2){
00960 (*(evtRecTrkCol->begin()+ipip[0]))->setQuality(1);
00961 (*(evtRecTrkCol->begin()+ipim[0]))->setQuality(1);
00962 (*(evtRecTrkCol->begin()+ipim[1]))->setQuality(1);
00963 (*(evtRecTrkCol->begin()+ikp[0]))->setQuality(1);
00964 }
00965
00966
00967
00968 setFilterPassed(true);
00970
00971 m_tuple->write();
00972
00973 return StatusCode::SUCCESS;
00974
00975 }
00976
00977
00978 StatusCode DQAKsKpi::finalize() {
00979
00980 MsgStream log(msgSvc(), name());
00981 log << MSG::INFO << "in finalize()" << endmsg;
00982 return StatusCode::SUCCESS;
00983 }
00984
00985