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 #include "EventModel/EventHeader.h"
00011
00012 #include "EvtRecEvent/EvtRecEvent.h"
00013 #include "EvtRecEvent/EvtRecTrack.h"
00014 #include "DstEvent/TofHitStatus.h"
00015
00016
00017 #include "TH1F.h"
00018 #include "VertexFit/KinematicFit.h"
00019 #include "VertexFit/VertexFit.h"
00020 #include "VertexFit/SecondVertexFit.h"
00021 #include "VertexFit/IVertexDbSvc.h"
00022
00023 #include "ParticleID/ParticleID.h"
00024
00025
00026 #include "TMath.h"
00027 #include "GaudiKernel/INTupleSvc.h"
00028 #include "GaudiKernel/NTuple.h"
00029 #include "GaudiKernel/Bootstrap.h"
00030 #include "GaudiKernel/ISvcLocator.h"
00031 #include "GaudiKernel/ITHistSvc.h"
00032
00033 #include "CLHEP/Vector/ThreeVector.h"
00034 #include "CLHEP/Vector/LorentzVector.h"
00035 #include "CLHEP/Vector/TwoVector.h"
00036 using CLHEP::Hep3Vector;
00037 using CLHEP::Hep2Vector;
00038 using CLHEP::HepLorentzVector;
00039 #include "CLHEP/Geometry/Point3D.h"
00040 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00041 typedef HepGeom::Point3D<double> HepPoint3D;
00042 #endif
00043 #include "DQAJpsi2PPbarAlg/DQAJpsi2PPbarAlg.h"
00044
00045 #include "VertexFit/KinematicFit.h"
00046 #include "VertexFit/VertexFit.h"
00047 #include "ParticleID/ParticleID.h"
00048
00049 #include <vector>
00050
00051
00052 const double mpi0 = 0.134977;
00053 const double mks0 = 0.497648;
00054 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00055
00056 const double velc = 299.792458;
00057 typedef std::vector<int> Vint;
00058 typedef std::vector<HepLorentzVector> Vp4;
00059
00060
00061
00062 const HepLorentzVector p_cms(0.034067,0.0,0.0,3.097);
00063 const Hep3Vector u_cms = -p_cms.boostVector();
00064
00065
00066
00067 static int counter[10]={0,0,0,0,0,0,0,0,0,0};
00069
00070 DQAJpsi2PPbarAlg::DQAJpsi2PPbarAlg(const std::string& name, ISvcLocator* pSvcLocator) :
00071 Algorithm(name, pSvcLocator) {
00072
00073
00074 declareProperty("Vr0cut", m_vr0cut=5.0);
00075 declareProperty("Vz0cut", m_vz0cut=20.0);
00076 declareProperty("Vr1cut", m_vr1cut=1.0);
00077 declareProperty("Vz1cut", m_vz1cut=10.0);
00078 declareProperty("Vctcut", m_cthcut=0.93);
00079 declareProperty("EnergyThreshold", m_energyThreshold=0.04);
00080 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
00081 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
00082 declareProperty("Test4C", m_test4C = 1);
00083 declareProperty("Test5C", m_test5C = 1);
00084 declareProperty("CheckDedx", m_checkDedx = 1);
00085 declareProperty("CheckTof", m_checkTof = 1);
00086 declareProperty("GammaAngCut", m_gammaAngCut=20.0);
00087 }
00088
00089
00090 StatusCode DQAJpsi2PPbarAlg::initialize(){
00091 MsgStream log(msgSvc(), name());
00092
00093 log << MSG::INFO << "in initialize()" << endmsg;
00094
00095 StatusCode status;
00096
00097
00098
00099 if(service("THistSvc", m_thsvc).isFailure()) {
00100 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
00101 return StatusCode::FAILURE;
00102 }
00103
00104
00105
00106 TH1F* hbst_p = new TH1F("bst_p", "bst_p", 80, 1.15, 1.31);
00107 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/hbst_p", hbst_p).isFailure()) {
00108 log << MSG::ERROR << "Couldn't register bst_p" << endreq;
00109 }
00110
00111 TH1F* hbst_cos = new TH1F("bst_cos", "bst_cos", 20, -1.0, 1.0);
00112 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/hbst_cos", hbst_cos).isFailure()) {
00113 log << MSG::ERROR << "Couldn't register bst_cos" << endreq;
00114 }
00115
00116 TH1F* hmpp = new TH1F("mpp", "mpp", 100, 3.05, 3.15);
00117 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/hmpp", hmpp).isFailure()) {
00118 log << MSG::ERROR << "Couldn't register mpp" << endreq;
00119 }
00120
00121 TH1F* hangle = new TH1F("angle", "angle", 50, 175.0, 180.);
00122 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/hangle", hangle).isFailure()) {
00123 log << MSG::ERROR << "Couldn't register angle" << endreq;
00124 }
00125
00126 TH1F* hdeltatof = new TH1F("deltatof", "deltatof", 50, -10., 10.);
00127 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/hdeltatof", hdeltatof).isFailure()) {
00128 log << MSG::ERROR << "Couldn't register deltatof" << endreq;
00129 }
00130
00131 TH1F* he_emc1 = new TH1F("e_emc1", "e_emc1", 100, 0.0, 2.0);
00132 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/he_emc1", he_emc1).isFailure()) {
00133 log << MSG::ERROR << "Couldn't register e_emc1" << endreq;
00134 }
00135
00136 TH1F* he_emc2 = new TH1F("e_emc2", "e_emc2", 100, 0.0, 2.0);
00137 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/he_emc2", he_emc2).isFailure()) {
00138 log << MSG::ERROR << "Couldn't register e_emc2" << endreq;
00139 }
00140
00141
00142
00143
00144
00145 NTuplePtr nt1(ntupleSvc(), "DQAFILE/DQAJpsi2PPbar");
00146 if ( nt1 ) m_tuple = nt1;
00147 else {
00148 m_tuple = ntupleSvc()->book ("DQAFILE/DQAJpsi2PPbar", CLID_ColumnWiseTuple, "N-Tuple");
00149 if ( m_tuple ) {
00150 status = m_tuple->addItem ("runNo", m_runNo);
00151 status = m_tuple->addItem ("event", m_event);
00152 status = m_tuple->addItem ("Nchrg", m_nchrg);
00153 status = m_tuple->addItem ("Nneu", m_nneu);
00154 status = m_tuple->addItem ("ngch", m_ngch, 0, 10);
00155
00156 status = m_tuple->addIndexedItem ("charge",m_ngch, m_charge);
00157 status = m_tuple->addIndexedItem ("vx0", m_ngch, m_vx0);
00158 status = m_tuple->addIndexedItem ("vy0", m_ngch, m_vy0);
00159 status = m_tuple->addIndexedItem ("vz0", m_ngch, m_vz0);
00160 status = m_tuple->addIndexedItem ("vr0", m_ngch, m_vr0);
00161
00162 status = m_tuple->addIndexedItem ("vx", m_ngch, m_vx);
00163 status = m_tuple->addIndexedItem ("vy", m_ngch, m_vy);
00164 status = m_tuple->addIndexedItem ("vz", m_ngch, m_vz);
00165 status = m_tuple->addIndexedItem ("vr", m_ngch, m_vr);
00166
00167 status = m_tuple->addIndexedItem ("px", m_ngch, m_px);
00168 status = m_tuple->addIndexedItem ("py", m_ngch, m_py);
00169 status = m_tuple->addIndexedItem ("pz", m_ngch, m_pz);
00170 status = m_tuple->addIndexedItem ("p", m_ngch, m_p );
00171 status = m_tuple->addIndexedItem ("cos", m_ngch, m_cos);
00172
00173 status = m_tuple->addIndexedItem ("bst_px", m_ngch, m_bst_px) ;
00174 status = m_tuple->addIndexedItem ("bst_py", m_ngch, m_bst_py) ;
00175 status = m_tuple->addIndexedItem ("bst_pz", m_ngch, m_bst_pz) ;
00176 status = m_tuple->addIndexedItem ("bst_p", m_ngch, m_bst_p) ;
00177 status = m_tuple->addIndexedItem ("bst_cos", m_ngch, m_bst_cos);
00178
00179
00180 status = m_tuple->addIndexedItem ("chie" , m_ngch, m_chie) ;
00181 status = m_tuple->addIndexedItem ("chimu" , m_ngch, m_chimu) ;
00182 status = m_tuple->addIndexedItem ("chipi" , m_ngch, m_chipi) ;
00183 status = m_tuple->addIndexedItem ("chik" , m_ngch, m_chik) ;
00184 status = m_tuple->addIndexedItem ("chip" , m_ngch, m_chip) ;
00185 status = m_tuple->addIndexedItem ("ghit" , m_ngch, m_ghit) ;
00186 status = m_tuple->addIndexedItem ("thit" , m_ngch, m_thit) ;
00187 status = m_tuple->addIndexedItem ("probPH" , m_ngch, m_probPH) ;
00188 status = m_tuple->addIndexedItem ("normPH" , m_ngch, m_normPH) ;
00189
00190 status = m_tuple->addIndexedItem ("e_emc" , m_ngch, m_e_emc) ;
00191
00192
00193 status = m_tuple->addIndexedItem ("qual_etof" , m_ngch, m_qual_etof );
00194 status = m_tuple->addIndexedItem ("tof_etof" , m_ngch, m_tof_etof );
00195 status = m_tuple->addIndexedItem ("te_etof" , m_ngch, m_te_etof );
00196 status = m_tuple->addIndexedItem ("tmu_etof" , m_ngch, m_tmu_etof );
00197 status = m_tuple->addIndexedItem ("tpi_etof" , m_ngch, m_tpi_etof );
00198 status = m_tuple->addIndexedItem ("tk_etof" , m_ngch, m_tk_etof );
00199 status = m_tuple->addIndexedItem ("tp_etof" , m_ngch, m_tp_etof );
00200
00201 status = m_tuple->addIndexedItem ("qual_btof1", m_ngch, m_qual_btof1 );
00202 status = m_tuple->addIndexedItem ("tof_btof1" , m_ngch, m_tof_btof1 );
00203 status = m_tuple->addIndexedItem ("te_btof1" , m_ngch, m_te_btof1 );
00204 status = m_tuple->addIndexedItem ("tmu_btof1" , m_ngch, m_tmu_btof1 );
00205 status = m_tuple->addIndexedItem ("tpi_btof1" , m_ngch, m_tpi_btof1 );
00206 status = m_tuple->addIndexedItem ("tk_btof1" , m_ngch, m_tk_btof1 );
00207 status = m_tuple->addIndexedItem ("tp_btof1" , m_ngch, m_tp_btof1 );
00208
00209 status = m_tuple->addIndexedItem ("qual_btof2", m_ngch, m_qual_btof2 );
00210 status = m_tuple->addIndexedItem ("tof_btof2" , m_ngch, m_tof_btof2 );
00211 status = m_tuple->addIndexedItem ("te_btof2" , m_ngch, m_te_btof2 );
00212 status = m_tuple->addIndexedItem ("tmu_btof2" , m_ngch, m_tmu_btof2 );
00213 status = m_tuple->addIndexedItem ("tpi_btof2" , m_ngch, m_tpi_btof2 );
00214 status = m_tuple->addIndexedItem ("tk_btof2" , m_ngch, m_tk_btof2 );
00215 status = m_tuple->addIndexedItem ("tp_btof2" , m_ngch, m_tp_btof2 );
00216
00217 status = m_tuple->addIndexedItem ("dedx_pid", m_ngch, m_dedx_pid);
00218 status = m_tuple->addIndexedItem ("tof1_pid", m_ngch, m_tof1_pid);
00219 status = m_tuple->addIndexedItem ("tof2_pid", m_ngch, m_tof2_pid);
00220 status = m_tuple->addIndexedItem ("prob_pi", m_ngch, m_prob_pi );
00221 status = m_tuple->addIndexedItem ("prob_k", m_ngch, m_prob_k );
00222 status = m_tuple->addIndexedItem ("prob_p", m_ngch, m_prob_p );
00223
00224 status = m_tuple->addItem ( "np", m_np );
00225 status = m_tuple->addItem ( "npb", m_npb );
00226
00227 status = m_tuple->addItem ( "m2p", m_m2p );
00228 status = m_tuple->addItem ( "angle", m_angle );
00229 status = m_tuple->addItem ( "deltatof", m_deltatof );
00230
00231 status = m_tuple->addItem ( "vtx_m2p", m_vtx_m2p );
00232 status = m_tuple->addItem ( "vtx_angle", m_vtx_angle );
00233
00234 status = m_tuple->addItem ( "m_chi2_4c", m_chi2_4c );
00235 status = m_tuple->addItem ( "m_m2p_4c", m_m2p_4c );
00236 status = m_tuple->addItem ( "m_angle_4c", m_angle_4c );
00237
00238 }
00239 else {
00240 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple) << endmsg;
00241 return StatusCode::FAILURE;
00242 }
00243 }
00244
00245
00246
00247
00248
00249 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00250 return StatusCode::SUCCESS;
00251
00252 }
00253
00254
00255 StatusCode DQAJpsi2PPbarAlg::execute() {
00256
00257
00258 MsgStream log(msgSvc(), name());
00259 log << MSG::INFO << "in execute()" << endreq;
00260
00261 setFilterPassed(false);
00262
00263 counter[0]++;
00264 log << MSG::INFO << "counter[0]=" << counter[0]<< endreq;
00265
00266 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00267 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00268 m_runNo = eventHeader->runNumber();
00269 m_event = eventHeader->eventNumber();
00270 m_nchrg = evtRecEvent->totalCharged();
00271 m_nneu = evtRecEvent->totalNeutral();
00272
00273 log << MSG::INFO <<"ncharg, nneu, tottks = "
00274 << evtRecEvent->totalCharged() << " , "
00275 << evtRecEvent->totalNeutral() << " , "
00276 << evtRecEvent->totalTracks() <<endreq;
00277
00278 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00279
00280
00281
00282
00283
00284 Vint iGood, iptrk, imtrk;
00285 iGood.clear();
00286 iptrk.clear();
00287 imtrk.clear();
00288 int nCharge = 0;
00289 Hep3Vector xorigin(0,0,0);
00290
00291
00292 IVertexDbSvc* vtxsvc;
00293 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00294 if (vtxsvc->isVertexValid()){
00295 double* dbv = vtxsvc->PrimaryVertex();
00296 double* vv = vtxsvc->SigmaPrimaryVertex();
00297
00298
00299 xorigin.setX(dbv[0]);
00300 xorigin.setY(dbv[1]);
00301 xorigin.setZ(dbv[2]);
00302 log << MSG::INFO
00303 <<"xorigin.x="<<xorigin.x()<<", "
00304 <<"xorigin.y="<<xorigin.y()<<", "
00305 <<"xorigin.z="<<xorigin.z()<<", "
00306 <<endreq ;
00307 }
00308
00309 for (int i = 0; i < evtRecEvent->totalCharged(); i++){
00310 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00311 if(!(*itTrk)->isMdcTrackValid()) continue;
00312 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00313 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00314 double x0=mdcTrk->x();
00315 double y0=mdcTrk->y();
00316 double z0=mdcTrk->z();
00317 double phi0=mdcTrk->helix(1);
00318 double xv=xorigin.x();
00319 double yv=xorigin.y();
00320 double zv=xorigin.z();
00321 double rv=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
00322 double cost = cos(mdcTrk->theta());
00323 if(fabs(z0-zv) >= m_vz1cut) continue;
00324 if(fabs(rv) >= m_vr1cut) continue;
00325
00326
00327 iGood.push_back(i);
00328 nCharge += mdcTrk->charge();
00329
00330 if (mdcTrk->charge() > 0) {
00331 iptrk.push_back(i);
00332 } else {
00333 imtrk.push_back(i);
00334 }
00335 }
00336
00337
00338
00339 int nGood = iGood.size();
00340 m_ngch = iGood.size();
00341 log << MSG::INFO << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
00342 if((nGood != 2) || (nCharge!=0) ){
00343 return StatusCode::SUCCESS;
00344 }
00345 counter[1]++;
00346 log << MSG::INFO << "counter[1]=" << counter[1]<< endreq;
00347
00349
00351
00352 Vint ipp, ipm;
00353 ipp.clear();
00354 ipm.clear();
00355
00356 Vp4 p_pp, p_pm;
00357 p_pp.clear();
00358 p_pm.clear();
00359
00360 int ii=-1 ;
00361
00362 ParticleID *pid = ParticleID::instance();
00363 for(int i = 0; i < nGood; i++) {
00364
00365 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00366 if(!(*itTrk)->isMdcTrackValid()) continue;
00367 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00368 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00369 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
00370
00371
00372 pid->init();
00373 pid->setMethod(pid->methodProbability());
00374 pid->setChiMinCut(4);
00375 pid->setRecTrack(*itTrk);
00376
00377
00378 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
00379
00380 pid->identify(pid->onlyPionKaonProton());
00381
00382
00383
00384 pid->calculate();
00385 if(!(pid->IsPidInfoValid())) continue;
00386
00387 double prob_pi = pid->probPion();
00388 double prob_k = pid->probKaon();
00389 double prob_p = pid->probProton();
00390
00391
00392
00393 if (prob_p > prob_pi && prob_p > prob_k) {
00394
00395 HepLorentzVector pkaltrk;
00396
00397 mdcKalTrk->setPidType(RecMdcKalTrack::proton);
00398 pkaltrk.setPx(mdcKalTrk->px());
00399 pkaltrk.setPy(mdcKalTrk->py());
00400 pkaltrk.setPz(mdcKalTrk->pz());
00401 double p3 = pkaltrk.mag();
00402 pkaltrk.setE(sqrt(p3*p3+xmass[4]*xmass[4]));
00403
00404 if(mdcTrk->charge() >0) {
00405 ii = 0;
00406 ipp.push_back(iGood[i]);
00407 p_pp.push_back(pkaltrk);
00408 } else {
00409 ii = 1 ;
00410 ipm.push_back(iGood[i]);
00411 p_pm.push_back(pkaltrk);
00412 }
00413
00414 m_dedx_pid[ii] = pid->chiDedx(2);
00415 m_tof1_pid[ii] = pid->chiTof1(2);
00416 m_tof2_pid[ii] = pid->chiTof2(2);
00417 m_prob_pi[ii] = pid->probPion();
00418 m_prob_k[ii] = pid->probKaon();
00419 m_prob_p[ii] = pid->probProton();
00420
00421 }
00422 }
00423
00424 m_np = ipp.size();
00425 m_npb = ipm.size();
00426
00427
00428
00429 counter[2]++;
00430 log << MSG::INFO << "counter[2]=" << counter[2]<< endreq;
00431
00433
00435 Vp4 p_ptrk, p_mtrk;
00436 p_ptrk.clear(), p_mtrk.clear();
00437 RecMdcKalTrack *ppKalTrk = 0 , *pmKalTrk = 0 ;
00438
00439 ii = -1 ;
00440 for(int i = 0; i < m_ngch; i++ ){
00441 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + iGood[i];
00442 if(!(*itTrk)->isMdcTrackValid()) continue;
00443 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00444 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00445 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
00446
00447 mdcKalTrk->setPidType(RecMdcKalTrack::proton);
00448
00449 if (mdcTrk->charge() > 0 ) {
00450 ii = 0 ;
00451 ppKalTrk = mdcKalTrk ;
00452 }else{
00453 ii = 1 ;
00454 pmKalTrk = mdcKalTrk ;
00455 }
00456
00457 m_charge[ii] = mdcTrk->charge();
00458
00459 double x0=mdcKalTrk->x();
00460 double y0=mdcKalTrk->y();
00461 double z0=mdcKalTrk->z();
00462 double phi0=mdcTrk->helix(1);
00463 double xv=xorigin.x();
00464 double yv=xorigin.y();
00465 double zv=xorigin.z();
00466 double rv=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
00467
00468 m_vx0[ii] = x0 ;
00469 m_vy0[ii] = y0 ;
00470 m_vz0[ii] = z0 ;
00471 m_vr0[ii] = sqrt((x0*x0)+(y0*y0)) ;
00472
00473 m_vx[ii] = x0-xv ;
00474 m_vy[ii] = y0-yv ;
00475 m_vz[ii] = fabs(z0-zv) ;
00476 m_vr[ii] = fabs(rv) ;
00477
00478 mdcKalTrk->setPidType(RecMdcKalTrack::proton);
00479 m_px[ii] = mdcKalTrk->px();
00480 m_py[ii] = mdcKalTrk->py();
00481 m_pz[ii] = mdcKalTrk->pz();
00482 m_p[ii] = mdcKalTrk->p();
00483 m_cos[ii] = mdcKalTrk->pz()/mdcKalTrk->p();
00484 double temp_e = sqrt(m_p[ii]*m_p[ii] + xmass[4]*xmass[4]);
00485 HepLorentzVector temp_p(m_px[ii],m_py[ii],m_pz[ii],temp_e);
00486 if (i==0){
00487 p_ptrk.push_back(temp_p);
00488 }else{
00489 p_mtrk.push_back(temp_p);
00490 }
00491
00492
00493 double ptrk = mdcKalTrk->p() ;
00494 if ((*itTrk)->isMdcDedxValid()) {
00495 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
00496 m_chie[ii] = dedxTrk->chiE();
00497 m_chimu[ii] = dedxTrk->chiMu();
00498 m_chipi[ii] = dedxTrk->chiPi();
00499 m_chik[ii] = dedxTrk->chiK();
00500 m_chip[ii] = dedxTrk->chiP();
00501 m_ghit[ii] = dedxTrk->numGoodHits();
00502 m_thit[ii] = dedxTrk->numTotalHits();
00503 m_probPH[ii]= dedxTrk->probPH();
00504 m_normPH[ii]= dedxTrk->normPH();
00505
00506 }
00507
00508 if((*itTrk)->isEmcShowerValid()) {
00509 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00510 m_e_emc[ii] = emcTrk->energy();
00511 }
00512
00513
00514 if((*itTrk)->isTofTrackValid()) {
00515
00516 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
00517
00518 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
00519 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
00520 TofHitStatus *status = new TofHitStatus;
00521 status->setStatus((*iter_tof)->status());
00522
00523 if(!(status->is_barrel())){
00524 if( !(status->is_counter()) ) continue;
00525 if( status->layer()!=0 ) continue;
00526 double path=(*iter_tof)->path();
00527 double tof = (*iter_tof)->tof();
00528 double ph = (*iter_tof)->ph();
00529 double rhit = (*iter_tof)->zrhit();
00530 double qual = 0.0 + (*iter_tof)->quality();
00531 double cntr = 0.0 + (*iter_tof)->tofID();
00532 double texp[5];
00533 for(int j = 0; j < 5; j++) {
00534 double gb = ptrk/xmass[j];
00535 double beta = gb/sqrt(1+gb*gb);
00536 texp[j] = path /beta/velc;
00537 }
00538
00539 m_qual_etof[ii] = qual;
00540 m_tof_etof[ii] = tof ;
00541 m_te_etof[ii] = tof - texp[0];
00542 m_tmu_etof[ii] = tof - texp[1];
00543 m_tpi_etof[ii] = tof - texp[2];
00544 m_tk_etof[ii] = tof - texp[3];
00545 m_tp_etof[ii] = tof - texp[4];
00546 }
00547 else {
00548 if( !(status->is_counter()) ) continue;
00549 if(status->layer()==1){
00550 double path=(*iter_tof)->path();
00551 double tof = (*iter_tof)->tof();
00552 double ph = (*iter_tof)->ph();
00553 double rhit = (*iter_tof)->zrhit();
00554 double qual = 0.0 + (*iter_tof)->quality();
00555 double cntr = 0.0 + (*iter_tof)->tofID();
00556 double texp[5];
00557 for(int j = 0; j < 5; j++) {
00558 double gb = ptrk/xmass[j];
00559 double beta = gb/sqrt(1+gb*gb);
00560 texp[j] = path /beta/velc;
00561 }
00562
00563 m_qual_btof1[ii] = qual;
00564 m_tof_btof1[ii] = tof ;
00565 m_te_btof1[ii] = tof - texp[0];
00566 m_tmu_btof1[ii] = tof - texp[1];
00567 m_tpi_btof1[ii] = tof - texp[2];
00568 m_tk_btof1[ii] = tof - texp[3];
00569 m_tp_btof1[ii] = tof - texp[4];
00570
00571 }
00572
00573 if(status->layer()==2){
00574 double path=(*iter_tof)->path();
00575 double tof = (*iter_tof)->tof();
00576 double ph = (*iter_tof)->ph();
00577 double rhit = (*iter_tof)->zrhit();
00578 double qual = 0.0 + (*iter_tof)->quality();
00579 double cntr = 0.0 + (*iter_tof)->tofID();
00580 double texp[5];
00581 for(int j = 0; j < 5; j++) {
00582 double gb = ptrk/xmass[j];
00583 double beta = gb/sqrt(1+gb*gb);
00584 texp[j] = path /beta/velc;
00585 }
00586
00587 m_qual_btof2[ii] = qual;
00588 m_tof_btof2[ii] = tof ;
00589 m_te_btof2[ii] = tof - texp[0];
00590 m_tmu_btof2[ii] = tof - texp[1];
00591 m_tpi_btof2[ii] = tof - texp[2];
00592 m_tk_btof2[ii] = tof - texp[3];
00593 m_tp_btof2[ii] = tof - texp[4];
00594 }
00595 }
00596 delete status;
00597 }
00598 }
00599 }
00600 counter[3]++;
00601 log << MSG::INFO << "counter[3]=" << counter[3]<< endreq;
00602
00603
00604
00605
00606
00607
00608 p_ptrk[0].boost(u_cms);
00609 p_mtrk[0].boost(u_cms);
00610 for (int i=0; i<=1; i++ ) {
00611 HepLorentzVector p;
00612 if (i==0) p = p_ptrk[0];
00613 if (i==1) p = p_mtrk[0];
00614
00615 m_bst_px[i] = p.px();
00616 m_bst_py[i] = p.py();
00617 m_bst_pz[i] = p.pz();
00618 m_bst_p[i] = p.rho();
00619 m_bst_cos[i]= p.cosTheta();
00620 }
00621
00622 m_m2p = (p_ptrk[0] + p_mtrk[0]).m();
00623
00624
00625
00626 m_angle = (p_ptrk[0].vect()).angle((p_mtrk[0]).vect()) * 180.0/(CLHEP::pi) ;
00627
00628
00629
00630
00631 double deltatof = 0.0 ;
00632 if(m_tof_btof1[0]*m_tof_btof1[1] != 0.0) deltatof+=m_tof_btof1[0]-m_tof_btof1[1];
00633 if(m_tof_btof2[0]*m_tof_btof2[1] != 0.0) deltatof+=m_tof_btof2[0]-m_tof_btof2[1];
00634 m_deltatof = deltatof ;
00635
00636
00637
00638
00639
00640
00641
00642
00643 if (fabs(m_bst_p[0]-1.232)>0.02 ) {return StatusCode::SUCCESS ; }
00644 if (fabs(m_bst_p[1]-1.232)>0.02 ) {return StatusCode::SUCCESS ; }
00645 if (fabs(m_deltatof)>4.0) {return StatusCode::SUCCESS ; }
00646 if (m_angle<178.0) {return StatusCode::SUCCESS ; }
00647 if (m_e_emc[0]>0.7) {return StatusCode::SUCCESS ; }
00648
00649 counter[4]++;
00650 log << MSG::INFO << "counter[4]=" << counter[4]<< endreq;
00651
00652
00653
00654 (*(evtRecTrkCol->begin()+iptrk[0]))->tagProton();
00655 (*(evtRecTrkCol->begin()+imtrk[0]))->tagProton();
00656
00657
00658
00659
00660
00661
00662 (*(evtRecTrkCol->begin()+iptrk[0]))->setQuality(0);
00663 (*(evtRecTrkCol->begin()+imtrk[0]))->setQuality(0);
00664
00665 setFilterPassed(true);
00666
00667 TH1 *h(0);
00668 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hbst_p", h).isSuccess()) {
00669 h->Fill(m_bst_p[0]);
00670 } else {
00671 log << MSG::ERROR << "Couldn't retrieve hbst_p" << endreq;
00672 }
00673
00674 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hbst_cos", h).isSuccess()) {
00675 h->Fill(m_bst_cos[0]);
00676 } else {
00677 log << MSG::ERROR << "Couldn't retrieve hbst_cos" << endreq;
00678 }
00679
00680 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hmpp", h).isSuccess()) {
00681 h->Fill(m_m2p);
00682 } else {
00683 log << MSG::ERROR << "Couldn't retrieve hmpp" << endreq;
00684 }
00685
00686 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hangle", h).isSuccess()) {
00687 h->Fill(m_angle);
00688 } else {
00689 log << MSG::ERROR << "Couldn't retrieve hangle" << endreq;
00690 }
00691
00692 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hdeltatof", h).isSuccess()) {
00693 h->Fill(m_deltatof);
00694 } else {
00695 log << MSG::ERROR << "Couldn't retrieve hdeltatof" << endreq;
00696 }
00697
00698 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/he_emc1", h).isSuccess()) {
00699 h->Fill(m_e_emc[0]);
00700 } else {
00701 log << MSG::ERROR << "Couldn't retrieve he_emc1" << endreq;
00702 }
00703
00704 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/he_emc2", h).isSuccess()) {
00705 h->Fill(m_e_emc[1]);
00706 } else {
00707 log << MSG::ERROR << "Couldn't retrieve he_emc2" << endreq;
00708 }
00709
00710
00711
00712
00713
00714 m_tuple -> write();
00715
00716 return StatusCode::SUCCESS;
00717 }
00718
00719
00720
00721 StatusCode DQAJpsi2PPbarAlg::finalize() {
00722
00723 MsgStream log(msgSvc(), name());
00724 log << MSG::INFO << "in finalize()" << endmsg;
00725
00726 std::cout << "*********** Singal.cxx *****************" << std::endl;
00727 std::cout << "the totale events is " << counter[0] << std::endl;
00728 std::cout << "select good charged track " << counter[1] << std::endl;
00729 std::cout << "PID " << counter[2] << std::endl;
00730 std::cout << "inform. for charged track " << counter[3] << std::endl;
00731 std::cout << "after all selections " << counter[4] << std::endl;
00732 std::cout << "****************************************" << std::endl;
00733
00734 return StatusCode::SUCCESS;
00735 }
00736