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 #include "GaudiKernel/ISvcLocator.h"
00010
00011 #include "EventModel/EventModel.h"
00012 #include "EventModel/Event.h"
00013
00014 #include "EvtRecEvent/EvtRecEvent.h"
00015 #include "EvtRecEvent/EvtRecTrack.h"
00016 #include "DstEvent/TofHitStatus.h"
00017 #include "EventModel/EventHeader.h"
00018 #include "McTruth/McParticle.h"
00019
00020 #include "TH1F.h"
00021 #include "TMath.h"
00022 #include "GaudiKernel/INTupleSvc.h"
00023 #include "GaudiKernel/NTuple.h"
00024 #include "GaudiKernel/ITHistSvc.h"
00025
00026 #include "GaudiKernel/Bootstrap.h"
00027 #include "GaudiKernel/IHistogramSvc.h"
00028 #include "CLHEP/Vector/ThreeVector.h"
00029 #include "CLHEP/Vector/LorentzVector.h"
00030 #include "CLHEP/Vector/TwoVector.h"
00031 using CLHEP::Hep3Vector;
00032 using CLHEP::Hep2Vector;
00033 using CLHEP::HepLorentzVector;
00034 #include "CLHEP/Geometry/Point3D.h"
00035 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00036 typedef HepGeom::Point3D<double> HepPoint3D;
00037 #endif
00038 #include "DQARhopiAlg/DQARhopi.h"
00039
00040 #include "VertexFit/KinematicFit.h"
00041 #include "VertexFit/VertexFit.h"
00042 #include "ParticleID/ParticleID.h"
00043
00044 #include <vector>
00045 using namespace Event;
00046
00047
00048 const double mpi = 0.13957;
00049 const double mk = 0.493677;
00050 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00051
00052 const double velc = 299.792458;
00053 typedef std::vector<int> Vint;
00054 typedef std::vector<HepLorentzVector> Vp4;
00055
00056 const HepLorentzVector ecms(0.034,0,0,3.097);
00057 const Hep3Vector u_cms = -ecms.boostVector();
00058
00059 int Ncut0,Ncut1,Ncut2,Ncut3,Ncut4,Ncut5,Ncut6,Ncut7,Ncut8,Ncut9,Ncut10;
00060
00062
00063 DQARhopi::DQARhopi(const std::string& name, ISvcLocator* pSvcLocator) :
00064 Algorithm(name, pSvcLocator) {
00065
00066
00067 declareProperty("Vr0cut", m_vr0cut=1.0);
00068 declareProperty("Vz0cut", m_vz0cut=5.0);
00069 declareProperty("Vctcut", m_cthcut=0.93);
00070 declareProperty("EnergyThreshold", m_energyThreshold=0.05);
00071 declareProperty("GammaAngCut", m_gammaAngCut=25.0);
00072 declareProperty("Test4C", m_test4C = 1);
00073 declareProperty("Test5C", m_test5C = 1);
00074 declareProperty("CheckDedx", m_checkDedx = 1);
00075 declareProperty("CheckTof", m_checkTof = 1);
00076 }
00077
00078
00079 StatusCode DQARhopi::initialize(){
00080 MsgStream log(msgSvc(), name());
00081
00082 log << MSG::INFO << "in initialize()" << endmsg;
00083
00084 StatusCode status;
00085
00086 NTuplePtr nt4(ntupleSvc(), "DQAFILE/Rhopi");
00087 if ( nt4 ) m_tuple4 = nt4;
00088 else {
00089 m_tuple4 = ntupleSvc()->book ("DQAFILE/Rhopi", CLID_ColumnWiseTuple, "ks N-Tuple example");
00090 if ( m_tuple4 ) {
00091 status = m_tuple4->addItem ("run", m_run);
00092 status = m_tuple4->addItem ("rec", m_rec);
00093 status = m_tuple4->addItem ("nch", m_nch);
00094 status = m_tuple4->addItem ("nneu", m_nneu);
00095 status = m_tuple4->addItem ("chi1", m_chi1);
00096 status = m_tuple4->addItem ("mpi0", m_mpi0);
00097 status = m_tuple4->addItem ("mprho0", m_prho0);
00098 status = m_tuple4->addItem ("mprhop", m_prhop);
00099 status = m_tuple4->addItem ("mprhom", m_prhom);
00100 status = m_tuple4->addItem ("mgood", m_good);
00101 status = m_tuple4->addItem ("mgam", m_gam);
00102 status = m_tuple4->addItem ("mpip", m_pip);
00103 status = m_tuple4->addItem ("mpim", m_pim);
00104 status = m_tuple4->addItem ("m2gam", m_2gam);
00105 status = m_tuple4->addItem ("ngch", m_ngch, 0, 10);
00106 status = m_tuple4->addItem ("nggneu", m_nggneu,0, 10);
00107 status = m_tuple4->addItem ("moutpi0",m_outpi0);
00108 status = m_tuple4->addItem ("cosang", m_cosang);
00109 status = m_tuple4->addItem ("moutpip",m_outpip);
00110 status = m_tuple4->addItem ("moutpim",m_outpim);
00111 status = m_tuple4->addItem ("menpip", m_enpip);
00112 status = m_tuple4->addItem ("mdcpip", m_dcpip );
00113 status = m_tuple4->addItem ("menpim", m_enpim );
00114 status = m_tuple4->addItem ("mdcpim", m_dcpim );
00115 status = m_tuple4->addItem ("mpipf", m_pipf);
00116 status = m_tuple4->addItem ("mpimf", m_pimf);
00117 status = m_tuple4->addItem ("mpi0f", m_pi0f);
00118
00119 status = m_tuple4->addItem ("mpmax", m_pmax);
00120 status = m_tuple4->addItem ("mppx", m_ppx);
00121 status = m_tuple4->addItem ("mppy", m_ppy);
00122 status = m_tuple4->addItem ("mppz", m_ppz);
00123 status = m_tuple4->addItem ("mcosthep",m_costhep);
00124 status = m_tuple4->addItem ("mppxkal", m_ppxkal);
00125 status = m_tuple4->addItem ("mppykal", m_ppykal);
00126 status = m_tuple4->addItem ("mppzkal", m_ppzkal);
00127 status = m_tuple4->addItem ("mmpx", m_mpx);
00128 status = m_tuple4->addItem ("mmpy", m_mpy);
00129 status = m_tuple4->addItem ("mmpz", m_mpz);
00130 status = m_tuple4->addItem ("mcosthem",m_costhem);
00131 status = m_tuple4->addItem ("mmpxkal", m_mpxkal);
00132 status = m_tuple4->addItem ("mmpykal", m_mpykal);
00133 status = m_tuple4->addItem ("mmpzkal", m_mpzkal);
00134
00135 status = m_tuple4->addItem ("mvxpin", m_vxpin);
00136 status = m_tuple4->addItem ("mvypin", m_vypin);
00137 status = m_tuple4->addItem ("mvzpin", m_vzpin);
00138 status = m_tuple4->addItem ("mvrpin", m_vrpin);
00139 status = m_tuple4->addItem ("mcosthepin",m_costhepin);
00140 status = m_tuple4->addItem ("mvxmin", m_vxmin);
00141 status = m_tuple4->addItem ("mvymin", m_vymin);
00142 status = m_tuple4->addItem ("mvzmin", m_vzmin);
00143 status = m_tuple4->addItem ("mvrmin", m_vrmin);
00144 status = m_tuple4->addItem ("mcosthemin",m_costhemin);
00145 status = m_tuple4->addItem ("mvxp", m_vxp);
00146 status = m_tuple4->addItem ("mvyp", m_vyp);
00147 status = m_tuple4->addItem ("mvzp", m_vzp);
00148 status = m_tuple4->addItem ("mvrp", m_vrp);
00149 status = m_tuple4->addItem ("mvxm", m_vxm);
00150 status = m_tuple4->addItem ("mvym", m_vym);
00151 status = m_tuple4->addItem ("mvzm", m_vzm);
00152 status = m_tuple4->addItem ("mvrm", m_vrm);
00153 status = m_tuple4->addItem ("nangecc",m_nangecc,0,10);
00154 status = m_tuple4->addIndexedItem ("mdthec", m_nangecc, m_dthec);
00155 status = m_tuple4->addIndexedItem ("mdphic", m_nangecc, m_dphic);
00156 status = m_tuple4->addIndexedItem ("mdangc", m_nangecc, m_dangc);
00157 status = m_tuple4->addIndexedItem ("mspippim", m_nangecc, m_mspippim);
00158
00159 status = m_tuple4->addItem ("angsg", dangsg);
00160 status = m_tuple4->addItem ("thesg", dthesg);
00161 status = m_tuple4->addItem ("phisg", dphisg);
00162 status = m_tuple4->addItem ("cosgth1", cosgth1);
00163 status = m_tuple4->addItem ("cosgth2", cosgth2);
00164
00165 status = m_tuple4->addItem ("mchi5", m_chi5);
00166 status = m_tuple4->addItem ("mkpi0", m_kpi0);
00167 status = m_tuple4->addItem ("mkpkm", m_kpkm);
00168 status = m_tuple4->addItem ("mkpp0", m_kpp0);
00169 status = m_tuple4->addItem ("mkmp0", m_kmp0);
00170 status = m_tuple4->addItem ("mpgam2pi1",m_pgam2pi1);
00171 status = m_tuple4->addItem ("mpgam2pi2",m_pgam2pi2);
00172 status = m_tuple4->addItem ("cosva1", cosva1);
00173 status = m_tuple4->addItem ("cosva2", cosva2);
00174 status = m_tuple4->addItem ("laypi1", m_laypi1);
00175 status = m_tuple4->addItem ("hit1", m_hit1);
00176 status = m_tuple4->addItem ("laypi2", m_laypi2);
00177 status = m_tuple4->addItem ("hit2", m_hit2);
00178 status = m_tuple4->addItem ("anglepm", m_anglepm);
00179
00180 status = m_tuple4->addIndexedItem ("mptrk", m_ngch, m_ptrk);
00181 status = m_tuple4->addIndexedItem ("chie", m_ngch, m_chie);
00182 status = m_tuple4->addIndexedItem ("chimu", m_ngch,m_chimu);
00183 status = m_tuple4->addIndexedItem ("chipi", m_ngch,m_chipi);
00184 status = m_tuple4->addIndexedItem ("chik", m_ngch,m_chik);
00185 status = m_tuple4->addIndexedItem ("chip", m_ngch,m_chip);
00186 status = m_tuple4->addIndexedItem ("probPH", m_ngch,m_probPH);
00187 status = m_tuple4->addIndexedItem ("normPH", m_ngch,m_normPH);
00188 status = m_tuple4->addIndexedItem ("ghit", m_ngch,m_ghit);
00189 status = m_tuple4->addIndexedItem ("thit", m_ngch,m_thit);
00190
00191 status = m_tuple4->addIndexedItem ("ptot_etof", m_ngch,m_ptot_etof);
00192 status = m_tuple4->addIndexedItem ("cntr_etof", m_ngch,m_cntr_etof);
00193 status = m_tuple4->addIndexedItem ("te_etof", m_ngch,m_te_etof);
00194 status = m_tuple4->addIndexedItem ("tmu_etof", m_ngch,m_tmu_etof);
00195 status = m_tuple4->addIndexedItem ("tpi_etof", m_ngch,m_tpi_etof);
00196 status = m_tuple4->addIndexedItem ("tk_etof", m_ngch,m_tk_etof);
00197 status = m_tuple4->addIndexedItem ("tp_etof", m_ngch,m_tp_etof);
00198 status = m_tuple4->addIndexedItem ("ph_etof", m_ngch,m_ph_etof);
00199 status = m_tuple4->addIndexedItem ("rhit_etof", m_ngch,m_rhit_etof);
00200 status = m_tuple4->addIndexedItem ("qual_etof", m_ngch,m_qual_etof);
00201 status = m_tuple4->addIndexedItem ("ec_toff_e", m_ngch,m_ec_toff_e);
00202 status = m_tuple4->addIndexedItem ("ec_toff_mu",m_ngch,m_ec_toff_mu);
00203 status = m_tuple4->addIndexedItem ("ec_toff_pi",m_ngch,m_ec_toff_pi);
00204 status = m_tuple4->addIndexedItem ("ec_toff_k", m_ngch,m_ec_toff_k);
00205 status = m_tuple4->addIndexedItem ("ec_toff_p", m_ngch,m_ec_toff_p);
00206 status = m_tuple4->addIndexedItem ("ec_tsig_e", m_ngch,m_ec_tsig_e);
00207 status = m_tuple4->addIndexedItem ("ec_tsig_mu",m_ngch,m_ec_tsig_mu);
00208 status = m_tuple4->addIndexedItem ("ec_tsig_pi",m_ngch,m_ec_tsig_pi);
00209 status = m_tuple4->addIndexedItem ("ec_tsig_k", m_ngch,m_ec_tsig_k);
00210 status = m_tuple4->addIndexedItem ("ec_tsig_p", m_ngch,m_ec_tsig_p);
00211 status = m_tuple4->addIndexedItem ("ec_tof", m_ngch,m_ec_tof);
00212 status = m_tuple4->addIndexedItem ("ptot_btof1",m_ngch,m_ptot_btof1);
00213 status = m_tuple4->addIndexedItem ("cntr_btof1",m_ngch,m_cntr_btof1);
00214 status = m_tuple4->addIndexedItem ("te_btof1", m_ngch,m_te_btof1);
00215 status = m_tuple4->addIndexedItem ("tmu_btof1", m_ngch,m_tmu_btof1);
00216 status = m_tuple4->addIndexedItem ("tpi_btof1", m_ngch,m_tpi_btof1);
00217 status = m_tuple4->addIndexedItem ("tk_btof1", m_ngch,m_tk_btof1);
00218 status = m_tuple4->addIndexedItem ("tp_btof1", m_ngch,m_tp_btof1);
00219 status = m_tuple4->addIndexedItem ("ph_btof1", m_ngch,m_ph_btof1);
00220 status = m_tuple4->addIndexedItem ("zhit_btof1",m_ngch,m_zhit_btof1);
00221 status = m_tuple4->addIndexedItem ("qual_btof1",m_ngch,m_qual_btof1);
00222 status = m_tuple4->addIndexedItem ("b1_toff_e", m_ngch,m_b1_toff_e);
00223 status = m_tuple4->addIndexedItem ("b1_toff_mu",m_ngch,m_b1_toff_mu);
00224 status = m_tuple4->addIndexedItem ("b1_toff_pi",m_ngch,m_b1_toff_pi);
00225 status = m_tuple4->addIndexedItem ("b1_toff_k", m_ngch,m_b1_toff_k);
00226 status = m_tuple4->addIndexedItem ("b1_toff_p", m_ngch,m_b1_toff_p);
00227 status = m_tuple4->addIndexedItem ("b1_tsig_e", m_ngch,m_b1_tsig_e);
00228 status = m_tuple4->addIndexedItem ("b1_tsig_mu",m_ngch,m_b1_tsig_mu);
00229 status = m_tuple4->addIndexedItem ("b1_tsig_pi",m_ngch,m_b1_tsig_pi);
00230 status = m_tuple4->addIndexedItem ("b1_tsig_k", m_ngch,m_b1_tsig_k);
00231 status = m_tuple4->addIndexedItem ("b1_tsig_p", m_ngch,m_b1_tsig_p);
00232 status = m_tuple4->addIndexedItem ("b1_tof", m_ngch,m_b1_tof);
00233
00234 status = m_tuple4->addIndexedItem ("mdedx_pid", m_ngch,m_dedx_pid);
00235 status = m_tuple4->addIndexedItem ("mtof1_pid", m_ngch,m_tof1_pid);
00236 status = m_tuple4->addIndexedItem ("mtof2_pid", m_ngch,m_tof2_pid);
00237 status = m_tuple4->addIndexedItem ("mprob_pid", m_ngch,m_prob_pid);
00238 status = m_tuple4->addIndexedItem ("mptrk_pid", m_ngch,m_ptrk_pid);
00239 status = m_tuple4->addIndexedItem ("mcost_pid", m_ngch,m_cost_pid);
00240 status = m_tuple4->addItem ("mpnp", m_pnp);
00241 status = m_tuple4->addItem ("mpnm", m_pnm);
00242
00243 status = m_tuple4->addIndexedItem ("numHits", m_nggneu,m_numHits);
00244 status = m_tuple4->addIndexedItem ("secondmoment", m_nggneu,m_secondmoment);
00245 status = m_tuple4->addIndexedItem ("mx", m_nggneu,m_x);
00246 status = m_tuple4->addIndexedItem ("my", m_nggneu,m_y);
00247 status = m_tuple4->addIndexedItem ("mz", m_nggneu,m_z);
00248 status = m_tuple4->addIndexedItem ("cosemc", m_nggneu,m_cosemc);
00249 status = m_tuple4->addIndexedItem ("phiemc", m_nggneu,m_phiemc);
00250 status = m_tuple4->addIndexedItem ("energy", m_nggneu,m_energy);
00251 status = m_tuple4->addIndexedItem ("eseed", m_nggneu,m_eSeed);
00252 status = m_tuple4->addIndexedItem ("me9", m_nggneu,m_e3x3);
00253 status = m_tuple4->addIndexedItem ("me25", m_nggneu,m_e5x5);
00254 status = m_tuple4->addIndexedItem ("mlat", m_nggneu,m_lat);
00255 status = m_tuple4->addIndexedItem ("ma20", m_nggneu,m_a20);
00256 status = m_tuple4->addIndexedItem ("ma42", m_nggneu,m_a42);
00257
00258 }
00259 else {
00260 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple4) << endmsg;
00261 return StatusCode::FAILURE;
00262 }
00263 }
00264
00265 if(service("THistSvc", m_thsvc).isFailure()) {
00266 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
00267 return StatusCode::FAILURE;
00268 }
00269
00270
00271
00272
00273 TH1F* mrho0 = new TH1F("mrho0", "mass for #rho^{0}->#pi^{+}#pi^{-}", 160, 0.0, 3.2);
00274 if(m_thsvc->regHist("/DQAHist/Rhopi/mrho0", mrho0).isFailure()) {
00275 log << MSG::ERROR << "Couldn't register mrho0" << endreq;
00276 }
00277
00278 TH1F* mrhop = new TH1F("mrhop", "mass for #rho^{+}->#pi^{+}#pi^{0}", 160, 0.0,3.2);
00279 if(m_thsvc->regHist("/DQAHist/Rhopi/mrhop", mrhop).isFailure()) {
00280 log << MSG::ERROR << "Couldn't register mrhop" << endreq;
00281 }
00282
00283 TH1F* mrhom = new TH1F("mrhom", "mass for #rho^{-}->#pi^{-}#pi^{0}", 160, 0.0, 3.2);
00284 if(m_thsvc->regHist("/DQAHist/Rhopi/mrhom", mrhom).isFailure()) {
00285 log << MSG::ERROR << "Couldn't register mrhom" << endreq;
00286 }
00287
00288
00289 TH1F* mpi0 = new TH1F("mpi0", "mass for #pi^{0}->#gamma #gamma", 50,0.08, 0.18);
00290 if(m_thsvc->regHist("/DQAHist/Rhopi/mpi0", mpi0).isFailure()) {
00291 log << MSG::ERROR << "Couldn't register mpi0" << endreq;
00292 }
00293
00294
00295
00296
00297
00298 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00299 return StatusCode::SUCCESS;
00300
00301 }
00302
00303
00304 StatusCode DQARhopi::execute() {
00305
00306
00307
00308 MsgStream log(msgSvc(), name());
00309 log << MSG::INFO << "in execute()" << endreq;
00310
00311 setFilterPassed(false);
00312
00313 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00314 int run=eventHeader->runNumber();
00315 int event=eventHeader->eventNumber();
00316 log << MSG::DEBUG <<"run, evtnum = "
00317 << run << " , "
00318 << event <<endreq;
00319
00320 m_run = eventHeader->runNumber();
00321 m_rec = eventHeader->eventNumber();
00322
00323
00324 Ncut0++;
00325 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00326 log << MSG::INFO << "get event tag OK" << endreq;
00327 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00328 << evtRecEvent->totalCharged() << " , "
00329 << evtRecEvent->totalNeutral() << " , "
00330 << evtRecEvent->totalTracks() <<endreq;
00331
00332 m_nch = evtRecEvent->totalCharged();
00333 m_nneu = evtRecEvent->totalNeutral();
00334
00335 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00336
00337
00338
00339
00340
00341 Vint iGood, ipip, ipim, ipnp,ipnm;
00342 iGood.clear();
00343 ipip.clear();
00344 ipim.clear();
00345 ipnp.clear();
00346 ipnm.clear();
00347 Vp4 ppip, ppim , ppnp, ppnm;
00348 ppip.clear();
00349 ppim.clear();
00350 ppnp.clear();
00351 ppnm.clear();
00352
00353 Hep3Vector xorigin(0,0,0);
00354
00355
00356
00357 IVertexDbSvc* vtxsvc;
00358 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00359 if(vtxsvc->isVertexValid()){
00360 double* dbv = vtxsvc->PrimaryVertex();
00361 double* vv = vtxsvc->SigmaPrimaryVertex();
00362 xorigin.setX(dbv[0]);
00363 xorigin.setY(dbv[1]);
00364 xorigin.setZ(dbv[2]);
00365 }
00366
00367 int nCharge = 0;
00368 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00369 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00370 if(!(*itTrk)->isMdcTrackValid()) continue;
00371 if (!(*itTrk)->isMdcKalTrackValid()) continue;
00372 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00373 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
00374
00375 double pch =mdcTrk->p();
00376 double x0 =mdcTrk->x();
00377 double y0 =mdcTrk->y();
00378 double z0 =mdcTrk->z();
00379 double phi0=mdcTrk->helix(1);
00380 double xv=xorigin.x();
00381 double yv=xorigin.y();
00382 double Rxy=fabs((x0-xv)*cos(phi0)+(y0-yv)*sin(phi0));
00383 double m_vx0 = x0;
00384 double m_vy0 = y0;
00385 double m_vz0 = z0-xorigin.z();
00386 double m_vr0 = Rxy;
00387 double m_Vct=cos(mdcTrk->theta());
00388
00389 if(fabs(m_vz0) >= m_vz0cut) continue;
00390 if(m_vr0 >= m_vr0cut) continue;
00391 if(fabs(m_Vct)>=m_cthcut) continue;
00392
00393 iGood.push_back((*itTrk)->trackId());
00394 nCharge += mdcTrk->charge();
00395
00396 }
00397
00398
00399
00400
00401 int nGood = iGood.size();
00402 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
00403 if((nGood != 2)||(nCharge!=0)){
00404 return StatusCode::SUCCESS;
00405 }
00406 Ncut1++;
00407
00408 Vint iGam;
00409 iGam.clear();
00410 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
00411 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00412 if(!(*itTrk)->isEmcShowerValid()) continue;
00413 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00414 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
00415
00416 double dthe = 200.;
00417 double dphi = 200.;
00418 double dang = 200.;
00419 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
00420 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
00421 if(!(*jtTrk)->isExtTrackValid()) continue;
00422 RecExtTrack *extTrk = (*jtTrk)->extTrack();
00423 if(extTrk->emcVolumeNumber() == -1) continue;
00424 Hep3Vector extpos = extTrk->emcPosition();
00425
00426 double angd = extpos.angle(emcpos);
00427 double thed = extpos.theta() - emcpos.theta();
00428 double phid = extpos.deltaPhi(emcpos);
00429 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00430 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00431
00432 if(fabs(thed) < fabs(dthe)) dthe = thed;
00433 if(fabs(phid) < fabs(dphi)) dphi = phid;
00434 if(angd < dang) dang = angd;
00435 }
00436 if(dang>=200) continue;
00437 double eraw = emcTrk->energy();
00438 dthe = dthe * 180 / (CLHEP::pi);
00439 dphi = dphi * 180 / (CLHEP::pi);
00440 dang = dang * 180 / (CLHEP::pi);
00441 double m_dthe = dthe;
00442 double m_dphi = dphi;
00443 double m_dang = dang;
00444 double m_eraw = eraw;
00445
00446 if(eraw < m_energyThreshold) continue;
00447 if(dang < m_gammaAngCut) continue;
00448
00449
00450
00451 iGam.push_back((*itTrk)->trackId());
00452 }
00453
00454
00455
00456
00457 int nGam = iGam.size();
00458
00459 log << MSG::DEBUG << "num Good Photon " << nGam << " , " <<evtRecEvent->totalNeutral()<<endreq;
00460 if(nGam<2){
00461 return StatusCode::SUCCESS;
00462 }
00463 Ncut2++;
00464
00465
00466
00467
00468
00469
00470 Vp4 pGam;
00471 pGam.clear();
00472 for(int i = 0; i < nGam; i++) {
00473 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
00474 RecEmcShower* emcTrk = (*itTrk)->emcShower();
00475 double eraw = emcTrk->energy();
00476 double phi = emcTrk->phi();
00477 double the = emcTrk->theta();
00478 HepLorentzVector ptrk;
00479 ptrk.setPx(eraw*sin(the)*cos(phi));
00480 ptrk.setPy(eraw*sin(the)*sin(phi));
00481 ptrk.setPz(eraw*cos(the));
00482 ptrk.setE(eraw);
00483
00484
00485
00486 pGam.push_back(ptrk);
00487 }
00488
00489 log << MSG::DEBUG << "liuf1 Good Photon " <<endreq;
00490
00491
00492 for(int i = 0; i < nGood; i++) {
00493 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00494 if(!(*itTrk)->isMdcTrackValid()) continue;
00495 if (!(*itTrk)->isMdcKalTrackValid()) continue;
00496 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00497 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
00498 RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);
00499
00500 if(mdcTrk->charge() >0) {
00501 ipip.push_back(iGood[i]);
00502 HepLorentzVector ptrk;
00503 ptrk.setPx(mdcKalTrk->px());
00504 ptrk.setPy(mdcKalTrk->py());
00505 ptrk.setPz(mdcKalTrk->pz());
00506 double p3 = ptrk.mag();
00507 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00508 ppip.push_back(ptrk);
00509 } else {
00510 ipim.push_back(iGood[i]);
00511 HepLorentzVector ptrk;
00512 ptrk.setPx(mdcKalTrk->px());
00513 ptrk.setPy(mdcKalTrk->py());
00514 ptrk.setPz(mdcKalTrk->pz());
00515 double p3 = ptrk.mag();
00516 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00517 ppim.push_back(ptrk);
00518 }
00519 }
00520
00521 int npip = ipip.size();
00522 int npim = ipim.size();
00523 log << MSG::DEBUG << "num of pion "<< ipip.size()<<","<<ipim.size()<<endreq;
00524 Ncut3++;
00525
00526
00527
00528
00529
00530
00531
00532 int langcc=0;
00533 double dthec = 200.;
00534 double dphic = 200.;
00535 double dangc = 200.;
00536 for(int i=0; i < npip; i++) {
00537 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + ipip[i] ;
00538 RecMdcTrack* mdcTrk1 = (*itTrk)->mdcTrack();
00539 RecMdcKalTrack* mdcKalTrk1 = (*itTrk)->mdcKalTrack();
00540 Hep3Vector emcpos(ppip[i].px(), ppip[i].py(), ppip[i].pz());
00541
00542 for(int j = 0; j < npim; j++) {
00543 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + ipim[j];
00544 RecMdcTrack* mdcTrk2 = (*jtTrk)->mdcTrack();
00545 RecMdcKalTrack* mdcKalTrk2 = (*jtTrk)->mdcKalTrack();
00546
00547 HepLorentzVector pippim=ppip[i]+ppim[j];
00548
00549 Hep3Vector extpos(ppim[j].px(), ppim[j].py(), ppim[j].pz());
00550
00551 double angd = extpos.angle(emcpos);
00552 double thed = extpos.theta() - emcpos.theta();
00553 double phid = extpos.deltaPhi(emcpos);
00554 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00555 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00556
00557 m_dthec[langcc] = thed * 180 / (CLHEP::pi);
00558 m_dphic[langcc] = phid * 180 / (CLHEP::pi);
00559 m_dangc[langcc] = angd * 180 / (CLHEP::pi);
00560 m_mspippim[langcc] =pippim.m();
00561 langcc++;
00562 }
00563 }
00564 m_nangecc=langcc;
00565
00566
00567
00568
00569
00570 double m_m2gg,m_momentpi0;
00571 HepLorentzVector pTot,p2g;
00572
00573 log << MSG::DEBUG << "liuf2 Good Photon " <<langcc<<endreq;
00574
00575
00576
00577
00578 double m_momentpip,m_momentpim,extmot1,extmot2;
00579 double emcTg1=0.0;
00580 double emcTg2=0.0;
00581 double nlaypi1=0;
00582 double nhit1=0;
00583 double nlaypi2=0;
00584 double nhit2=0;
00585
00586 EvtRecTrackIterator itTrk11 = evtRecTrkCol->begin() + ipip[0];
00587 RecMdcTrack* mdcTrk11 = (*itTrk11)->mdcTrack();
00588 RecMdcKalTrack* mdcKalTrk11 = (*itTrk11)->mdcKalTrack();
00589 RecEmcShower* emcTrk11 = (*itTrk11)->emcShower();
00590 RecMucTrack *mucTrk11=(*itTrk11)->mucTrack();
00591 double phi01=mdcTrk11->helix(1);
00592
00593 EvtRecTrackIterator itTrk12 = evtRecTrkCol->begin() + ipim[0];
00594 RecMdcTrack* mdcTrk12 = (*itTrk12)->mdcTrack();
00595 RecMdcKalTrack* mdcKalTrk12 = (*itTrk12)->mdcKalTrack();
00596 RecEmcShower* emcTrk12 = (*itTrk12)->emcShower();
00597 RecMucTrack *mucTrk12=(*itTrk12)->mucTrack();
00598 double phi02=mdcTrk12->helix(1);
00599
00600 m_vxpin = mdcTrk11->x();
00601 m_vypin = mdcTrk11->y();
00602 m_vzpin = mdcTrk11->z()-xorigin.z();
00603 m_vrpin = fabs((mdcTrk11->x()-xorigin.x())*cos(phi01)+(mdcTrk11->y()-xorigin.y())*sin(phi01));
00604 m_costhepin =cos(mdcTrk11->theta());
00605
00606 m_momentpip=mdcTrk11->p();
00607 m_ppx =mdcTrk11->px();
00608 m_ppy =mdcTrk11->py();
00609 m_ppz =mdcTrk11->pz();
00610
00611 m_vxp = mdcKalTrk11->x();
00612 m_vyp = mdcKalTrk11->y();
00613 m_vzp = mdcKalTrk11->z()-xorigin.z();
00614 m_vrp = fabs((mdcKalTrk11->x()-xorigin.x())*cos(phi01)+(mdcKalTrk11->y()-xorigin.y())*sin(phi01));
00615 m_costhep =cos(mdcKalTrk11->theta());
00616
00617 extmot1=mdcKalTrk11->p();
00618 m_ppxkal =mdcKalTrk11->px();
00619 m_ppykal =mdcKalTrk11->py();
00620 m_ppzkal =mdcKalTrk11->pz();
00621
00622 m_vxmin = mdcTrk12->x();
00623 m_vymin = mdcTrk12->y();
00624 m_vzmin = mdcTrk12->z();
00625 m_vrmin = fabs((mdcTrk12->x()-xorigin.x())*cos(phi02)+(mdcTrk12->y()-xorigin.y())*sin(phi02));
00626 m_costhemin =cos(mdcTrk12->theta());
00627
00628 m_momentpim=mdcTrk12->p();
00629 m_mpx =mdcTrk12->px();
00630 m_mpy =mdcTrk12->py();
00631 m_mpz =mdcTrk12->pz();
00632
00633 m_vxm = mdcKalTrk12->x();
00634 m_vym = mdcKalTrk12->y();
00635 m_vzm = mdcKalTrk12->z();
00636 m_vrm = fabs((mdcKalTrk12->x()-xorigin.x())*cos(phi02)+(mdcKalTrk12->y()-xorigin.y())*sin(phi02));
00637 m_costhem =cos(mdcKalTrk12->theta());
00638
00639 extmot2 =mdcKalTrk12->p();
00640 m_mpxkal =mdcKalTrk12->px();
00641 m_mpykal =mdcKalTrk12->py();
00642 m_mpzkal =mdcKalTrk12->pz();
00643
00644 if((*itTrk11)->isEmcShowerValid()){
00645 emcTg1 = emcTrk11->energy();
00646 }
00647 if((*itTrk12)->isEmcShowerValid()){
00648 emcTg2 = emcTrk12->energy();
00649 }
00650 if((*itTrk11)->isMucTrackValid()){
00651 nlaypi1=mucTrk11->numLayers();
00652 nhit1 =mucTrk11->numHits();
00653 }
00654 if((*itTrk12)->isMucTrackValid()){
00655 nlaypi2=mucTrk12->numLayers();
00656 nhit2 =mucTrk12->numHits();
00657 }
00658
00659 m_laypi1=nlaypi1;
00660 m_hit1 =nhit1;
00661 m_laypi2=nlaypi2;
00662 m_hit2 =nhit2;
00663
00664 log << MSG::DEBUG << "liuf3 Good Photon " << ipim[0] <<endreq;
00665
00666 RecMdcKalTrack *pipTrk = (*(evtRecTrkCol->begin()+ipip[0]))->mdcKalTrack();
00667 RecMdcKalTrack *pimTrk = (*(evtRecTrkCol->begin()+ipim[0]))->mdcKalTrack();
00668
00669 log << MSG::DEBUG << "liuf4 Good Photon " <<endreq;
00670
00671 WTrackParameter wvpipTrk, wvpimTrk,wvkpTrk, wvkmTrk;
00672 wvpipTrk = WTrackParameter(mpi, pipTrk->getZHelix(), pipTrk->getZError());
00673 wvpimTrk = WTrackParameter(mpi, pimTrk->getZHelix(), pimTrk->getZError());
00674
00675 wvkpTrk = WTrackParameter(mk, pipTrk->getZHelixK(), pipTrk->getZErrorK());
00676 wvkmTrk = WTrackParameter(mk, pimTrk->getZHelixK(), pimTrk->getZErrorK());
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688 HepPoint3D vx(0., 0., 0.);
00689 HepSymMatrix Evx(3, 0);
00690 double bx = 1E+6;
00691 double by = 1E+6;
00692 double bz = 1E+6;
00693 Evx[0][0] = bx*bx;
00694 Evx[1][1] = by*by;
00695 Evx[2][2] = bz*bz;
00696
00697 VertexParameter vxpar;
00698 vxpar.setVx(vx);
00699 vxpar.setEvx(Evx);
00700
00701
00702
00703
00704
00705 VertexFit* vtxfit = VertexFit::instance();
00706
00707
00708
00709
00710 double chi5=9999.0;
00711 double jkpi0=-0.5;
00712 double jkpkm=0.0;
00713 double jkpp0=0.0;
00714 double jkmp0=0.0;
00715 vtxfit->init();
00716 vtxfit->AddTrack(0, wvkpTrk);
00717 vtxfit->AddTrack(1, wvkmTrk);
00718 vtxfit->AddVertex(0, vxpar,0, 1);
00719 if(vtxfit->Fit(0)) {
00720 vtxfit->Swim(0);
00721 WTrackParameter wkp = vtxfit->wtrk(0);
00722 WTrackParameter wkm = vtxfit->wtrk(1);
00723
00724 KinematicFit * kmfit = KinematicFit::instance();
00725
00726
00727
00728
00729
00730 double chisq = 9999.;
00731 for(int i = 0; i < nGam-1; i++) {
00732 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
00733 for(int j = i+1; j < nGam; j++) {
00734 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[j]))->emcShower();
00735 kmfit->init();
00736 kmfit->setDynamicerror(1);
00737 kmfit->AddTrack(0, wkp);
00738 kmfit->AddTrack(1, wkm);
00739 kmfit->AddTrack(2, 0.0, g1Trk);
00740 kmfit->AddTrack(3, 0.0, g2Trk);
00741 kmfit->AddFourMomentum(0, ecms);
00742 bool oksq = kmfit->Fit();
00743 if(oksq) {
00744 double chi2 = kmfit->chisq();
00745 if(chi2 < chi5) {
00746 HepLorentzVector kpi0 = kmfit->pfit(2) + kmfit->pfit(3);
00747 HepLorentzVector kpkm = kmfit->pfit(0) + kmfit->pfit(1);
00748 HepLorentzVector kpp0 = kmfit->pfit(0) + kpi0;
00749 HepLorentzVector kmp0 = kmfit->pfit(1) + kpi0;
00750 chi5 = kmfit->chisq();
00751 jkpi0 = kpi0.m();
00752 jkpkm= kpkm.m();
00753 jkpp0= kpp0.m();
00754 jkmp0= kmp0.m();
00755 }
00756 }
00757 }
00758 }
00759 }
00760
00761
00762
00763
00764
00765 vtxfit->init();
00766 vtxfit->AddTrack(0, wvpipTrk);
00767 vtxfit->AddTrack(1, wvpimTrk);
00768 vtxfit->AddVertex(0, vxpar,0, 1);
00769 if(!vtxfit->Fit(0)) return SUCCESS;
00770 vtxfit->Swim(0);
00771
00772 WTrackParameter wpip = vtxfit->wtrk(0);
00773 WTrackParameter wpim = vtxfit->wtrk(1);
00774
00775 log << MSG::DEBUG << "liuf5 Good Photon " <<endreq;
00776
00777 KinematicFit * kmfit = KinematicFit::instance();
00778
00779
00780
00781
00782
00783 double chisq = 9999.;
00784 int ig1 = -1;
00785 int ig2 = -1;
00786 HepLorentzVector gg1,gg2;
00787 for(int i = 0; i < nGam-1; i++) {
00788 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+iGam[i]))->emcShower();
00789 for(int j = i+1; j < nGam; j++) {
00790 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+iGam[j]))->emcShower();
00791 kmfit->init();
00792 kmfit->setDynamicerror(1);
00793 kmfit->AddTrack(0, wpip);
00794 kmfit->AddTrack(1, wpim);
00795 kmfit->AddTrack(2, 0.0, g1Trk);
00796 kmfit->AddTrack(3, 0.0, g2Trk);
00797 kmfit->AddFourMomentum(0, ecms);
00798 bool oksq = kmfit->Fit();
00799 if(oksq) {
00800 double chi2 = kmfit->chisq();
00801 if(chi2 < chisq) {
00802 chisq = chi2;
00803 ig1 = iGam[i];
00804 ig2 = iGam[j];
00805 gg1 = pGam[i];
00806 gg2 = pGam[j];
00807 }
00808 }
00809 }
00810 }
00811
00812 p2g = gg1 + gg2;
00813 m_pmax=gg1.e()>gg2.e()?gg1.e():gg2.e();
00814 m_m2gg = p2g.m();
00815 m_cosang=(gg1.e()-gg2.e())/p2g.rho();
00816 m_momentpi0=sqrt(p2g.px()*p2g.px()+p2g.py()*p2g.py()+p2g.pz()*p2g.pz());
00817 log << MSG::DEBUG << " chisq for 4c " << chisq <<endreq;
00818 if(chisq > 200) {
00819 return StatusCode::SUCCESS;
00820 }
00821
00822
00823 Vint jGood;
00824 jGood.clear();
00825 jGood.push_back(ipip[0]);
00826 jGood.push_back(ipim[0]);
00827 m_ngch = jGood.size();
00828
00829 Vint jGgam;
00830 jGgam.clear();
00831 jGgam.push_back(ig1);
00832 jGgam.push_back(ig2);
00833 m_nggneu=jGgam.size();
00834
00835 HepLorentzVector ppip1=ppip[0];
00836 HepLorentzVector ppim1=ppim[0];
00837
00838 HepLorentzVector Ppipboost = ppip1.boost(u_cms);
00839 HepLorentzVector Ppimboost = ppim1.boost(u_cms);
00840 Hep3Vector p3pi1 = Ppipboost.vect();
00841 Hep3Vector p3pi2 = Ppimboost.vect();
00842 m_anglepm = (p3pi1.angle(p3pi2))* 180 / (CLHEP::pi);
00843
00844
00845
00846 RecEmcShower *g1Trk = (*(evtRecTrkCol->begin()+ig1))->emcShower();
00847 RecEmcShower *g2Trk = (*(evtRecTrkCol->begin()+ig2))->emcShower();
00848 kmfit->init();
00849 kmfit->setDynamicerror(1);
00850 kmfit->AddTrack(0, wpip);
00851 kmfit->AddTrack(1, wpim);
00852 kmfit->AddTrack(2, 0.0, g1Trk);
00853 kmfit->AddTrack(3, 0.0, g2Trk);
00854 kmfit->AddFourMomentum(0, ecms);
00855 bool oksq = kmfit->Fit();
00856 if(!oksq) return SUCCESS;
00857
00858 HepLorentzVector ppi0 = kmfit->pfit(2) + kmfit->pfit(3);
00859 HepLorentzVector prho0 = kmfit->pfit(0) + kmfit->pfit(1);
00860 HepLorentzVector prhop = kmfit->pfit(0) + ppi0;
00861 HepLorentzVector prhom = kmfit->pfit(1) + ppi0;
00862 HepLorentzVector pgam2pi1 = prho0 + kmfit->pfit(2);
00863 HepLorentzVector pgam2pi2 = prho0 + kmfit->pfit(3);
00864 HepLorentzVector pgam11 =kmfit->pfit(2);
00865 HepLorentzVector pgam12 =kmfit->pfit(3);
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876 m_chi1 = kmfit->chisq();
00877 m_mpi0 = ppi0.m();
00878 m_prho0= prho0.m();
00879 m_prhop= prhop.m();
00880 m_prhom= prhom.m();
00881 m_good = nGood;
00882 m_gam = nGam;
00883
00884 m_pip = npip;
00885 m_pim = npim;
00886 m_2gam = m_m2gg;
00887 m_outpi0=m_momentpi0;
00888 m_outpip=m_momentpip;
00889 m_outpim=m_momentpim;
00890 m_enpip =emcTg1;
00891 m_dcpip =extmot1;
00892 m_enpim =emcTg2;
00893 m_dcpim =extmot2;
00894 m_pipf=kmfit->pfit(0).rho();
00895 m_pimf=kmfit->pfit(1).rho();
00896 m_pi0f=ppi0.rho();
00897
00898 m_chi5=chi5;
00899 m_kpi0=jkpi0;
00900 m_kpkm=jkpkm;
00901 m_kpp0=jkpp0;
00902 m_kmp0=jkmp0;
00903 m_pgam2pi1=pgam2pi1.m();
00904 m_pgam2pi2=pgam2pi2.m();
00905 cosva1=kmfit->pfit(2).rho();
00906 cosva2=kmfit->pfit(3).rho();
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918 for(int ii = 0; ii < m_ngch; ii++) {
00919
00920 m_ptrk[ii] = 9999.0;
00921 m_chie[ii] = 9999.0;
00922 m_chimu[ii] = 9999.0;
00923 m_chipi[ii] = 9999.0;
00924 m_chik[ii] = 9999.0;
00925 m_chip[ii] = 9999.0;
00926 m_ghit[ii] = 9999.0;
00927 m_thit[ii] = 9999.0;
00928 m_probPH[ii] = 9999.0;
00929 m_normPH[ii] = 9999.0;
00930
00931
00932 m_cntr_etof[ii] = 9999.0;
00933 m_ptot_etof[ii] = 9999.0;
00934 m_ph_etof[ii] = 9999.0;
00935 m_rhit_etof[ii] = 9999.0;
00936 m_qual_etof[ii] = 9999.0;
00937 m_te_etof[ii] = 9999.0;
00938 m_tmu_etof[ii] = 9999.0;
00939 m_tpi_etof[ii] = 9999.0;
00940 m_tk_etof[ii] = 9999.0;
00941 m_tp_etof[ii] = 9999.0;
00942 m_ec_tof[ii] = 9999.0;
00943 m_ec_toff_e[ii] = 9999.0;
00944 m_ec_toff_mu[ii] = 9999.0;
00945 m_ec_toff_pi[ii] = 9999.0;
00946 m_ec_toff_k[ii] = 9999.0;
00947 m_ec_toff_p[ii] = 9999.0;
00948 m_ec_tsig_e[ii] = 9999.0;
00949 m_ec_tsig_mu[ii] = 9999.0;
00950 m_ec_tsig_pi[ii] = 9999.0;
00951 m_ec_tsig_k[ii] = 9999.0;
00952 m_ec_tsig_p[ii] = 9999.0;
00953
00954
00955 m_cntr_btof1[ii] = 9999.0;
00956 m_ptot_btof1[ii] = 9999.0;
00957 m_ph_btof1[ii] = 9999.0;
00958 m_zhit_btof1[ii] = 9999.0;
00959 m_qual_btof1[ii] = 9999.0;
00960 m_te_btof1[ii] = 9999.0;
00961 m_tmu_btof1[ii] = 9999.0;
00962 m_tpi_btof1[ii] = 9999.0;
00963 m_tk_btof1[ii] = 9999.0;
00964 m_tp_btof1[ii] = 9999.0;
00965 m_b1_tof[ii] = 9999.0;
00966 m_b1_toff_e[ii] = 9999.0;
00967 m_b1_toff_mu[ii] = 9999.0;
00968 m_b1_toff_pi[ii] = 9999.0;
00969 m_b1_toff_k[ii] = 9999.0;
00970 m_b1_toff_p[ii] = 9999.0;
00971 m_b1_tsig_e[ii] = 9999.0;
00972 m_b1_tsig_mu[ii] = 9999.0;
00973 m_b1_tsig_pi[ii] = 9999.0;
00974 m_b1_tsig_k[ii] = 9999.0;
00975 m_b1_tsig_p[ii] = 9999.0;
00976
00977 m_dedx_pid[ii] = 9999.0;
00978 m_tof1_pid[ii] = 9999.0;
00979 m_tof2_pid[ii] = 9999.0;
00980 m_prob_pid[ii] = 9999.0;
00981 m_ptrk_pid[ii] = 9999.0;
00982 m_cost_pid[ii] = 9999.0;
00983 }
00984
00985
00986 int indx0=0;
00987 for(int i = 0; i < m_ngch; i++) {
00988 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
00989 if(!(*itTrk)->isMdcTrackValid()) continue;
00990 if(!(*itTrk)->isMdcDedxValid())continue;
00991 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00992 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
00993 m_ptrk[indx0] = mdcTrk->p();
00994
00995 m_chie[indx0] = dedxTrk->chiE();
00996 m_chimu[indx0] = dedxTrk->chiMu();
00997 m_chipi[indx0] = dedxTrk->chiPi();
00998 m_chik[indx0] = dedxTrk->chiK();
00999 m_chip[indx0] = dedxTrk->chiP();
01000 m_ghit[indx0] = dedxTrk->numGoodHits();
01001 m_thit[indx0] = dedxTrk->numTotalHits();
01002 m_probPH[indx0] = dedxTrk->probPH();
01003 m_normPH[indx0] = dedxTrk->normPH();
01004 indx0++;
01005 }
01006
01007
01008
01009
01010
01011
01012
01013 int indx1=0;
01014 for(int i = 0; i < m_ngch; i++) {
01015 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
01016 if(!(*itTrk)->isMdcTrackValid()) continue;
01017 if(!(*itTrk)->isTofTrackValid()) continue;
01018
01019 RecMdcTrack * mdcTrk = (*itTrk)->mdcTrack();
01020 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
01021
01022 double ptrk = mdcTrk->p();
01023 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
01024 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
01025 TofHitStatus *status = new TofHitStatus;
01026 status->setStatus((*iter_tof)->status());
01027 if(!(status->is_barrel())){
01028 if( !(status->is_counter()) ) continue;
01029 if( status->layer()!=1 ) continue;
01030 double path=(*iter_tof)->path();
01031 double tof = (*iter_tof)->tof();
01032 double ph = (*iter_tof)->ph();
01033 double rhit = (*iter_tof)->zrhit();
01034 double qual = 0.0 + (*iter_tof)->quality();
01035 double cntr = 0.0 + (*iter_tof)->tofID();
01036 double texp[5];
01037 double tsig[5];
01038 for(int j = 0; j < 5; j++) {
01039 texp[j] = (*iter_tof)->texp(j);
01040
01041
01042 }
01043 m_cntr_etof[indx1] = cntr;
01044 m_ptot_etof[indx1] = ptrk;
01045 m_ph_etof[indx1] = ph;
01046 m_rhit_etof[indx1] = rhit;
01047 m_qual_etof[indx1] = qual;
01048 m_te_etof[indx1] = tof - texp[0];
01049 m_tmu_etof[indx1] = tof - texp[1];
01050 m_tpi_etof[indx1] = tof - texp[2];
01051 m_tk_etof[indx1] = tof - texp[3];
01052 m_tp_etof[indx1] = tof - texp[4];
01053
01054 m_ec_tof[indx1] = tof;
01055
01056 m_ec_toff_e[indx1] = (*iter_tof)->toffset(0);
01057 m_ec_toff_mu[indx1] = (*iter_tof)->toffset(1);
01058 m_ec_toff_pi[indx1] = (*iter_tof)->toffset(2);
01059 m_ec_toff_k[indx1] = (*iter_tof)->toffset(3);
01060 m_ec_toff_p[indx1] = (*iter_tof)->toffset(4);
01061
01062 m_ec_tsig_e[indx1] = (*iter_tof)->sigma(0);
01063 m_ec_tsig_mu[indx1] = (*iter_tof)->sigma(1);
01064 m_ec_tsig_pi[indx1] = (*iter_tof)->sigma(2);
01065 m_ec_tsig_k[indx1] = (*iter_tof)->sigma(3);
01066 m_ec_tsig_p[indx1] = (*iter_tof)->sigma(4);
01067
01068 }
01069 else {
01070 if( !(status->is_cluster()) ) continue;
01071 double path=(*iter_tof)->path();
01072 double tof = (*iter_tof)->tof();
01073 double ph = (*iter_tof)->ph();
01074 double rhit = (*iter_tof)->zrhit();
01075 double qual = 0.0 + (*iter_tof)->quality();
01076 double cntr = 0.0 + (*iter_tof)->tofID();
01077 double texp[5];
01078 for(int j = 0; j < 5; j++) {
01079 texp[j] = (*iter_tof)->texp(j);
01080 }
01081 m_cntr_btof1[indx1] = cntr;
01082 m_ptot_btof1[indx1] = ptrk;
01083 m_ph_btof1[indx1] = ph;
01084 m_zhit_btof1[indx1] = rhit;
01085 m_qual_btof1[indx1] = qual;
01086 m_te_btof1[indx1] = tof - texp[0];
01087 m_tmu_btof1[indx1] = tof - texp[1];
01088 m_tpi_btof1[indx1] = tof - texp[2];
01089 m_tk_btof1[indx1] = tof - texp[3];
01090 m_tp_btof1[indx1] = tof - texp[4];
01091
01092 m_b1_tof[indx1] = tof;
01093
01094 m_b1_toff_e[indx1] = (*iter_tof)->toffset(0);
01095 m_b1_toff_mu[indx1] = (*iter_tof)->toffset(1);
01096 m_b1_toff_pi[indx1] = (*iter_tof)->toffset(2);
01097 m_b1_toff_k[indx1] = (*iter_tof)->toffset(3);
01098 m_b1_toff_p[indx1] = (*iter_tof)->toffset(4);
01099
01100 m_b1_tsig_e[indx1] = (*iter_tof)->sigma(0);
01101 m_b1_tsig_mu[indx1] = (*iter_tof)->sigma(1);
01102 m_b1_tsig_pi[indx1] = (*iter_tof)->sigma(2);
01103 m_b1_tsig_k[indx1] = (*iter_tof)->sigma(3);
01104 m_b1_tsig_p[indx1] = (*iter_tof)->sigma(4);
01105
01106 }
01107 delete status;
01108 }
01109 indx1++;
01110 }
01111
01112
01113
01114
01115 int indx2=0;
01116 ParticleID *pid = ParticleID::instance();
01117 for(int i = 0; i < m_ngch; i++) {
01118 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
01119
01120 pid->init();
01121 pid->setMethod(pid->methodProbability());
01122
01123
01124 pid->setChiMinCut(4);
01125 pid->setRecTrack(*itTrk);
01126 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
01127 pid->identify(pid->onlyPion() | pid->onlyKaon());
01128
01129
01130 pid->calculate();
01131 if(!(pid->IsPidInfoValid())) continue;
01132 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
01133
01134 m_dedx_pid[indx2] = pid->chiDedx(2);
01135 m_tof1_pid[indx2] = pid->chiTof1(2);
01136 m_tof2_pid[indx2] = pid->chiTof2(2);
01137 m_prob_pid[indx2] = pid->probPion();
01138
01139
01140
01141
01142
01143 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
01144 RecMdcKalTrack::setPidType (RecMdcKalTrack::pion);
01145
01146 m_ptrk_pid[indx2] = mdcKalTrk->p();
01147 m_cost_pid[indx2] = cos(mdcTrk->theta());
01148
01149 if(mdcTrk->charge() >0 && pid->probPion() > pid->probKaon()) {
01150 ipnp.push_back(jGood[i]);
01151 HepLorentzVector ptrk;
01152 ptrk.setPx(mdcKalTrk->px());
01153 ptrk.setPy(mdcKalTrk->py());
01154 ptrk.setPz(mdcKalTrk->pz());
01155 double p3 = ptrk.mag();
01156 ptrk.setE(sqrt(p3*p3+mpi*mpi));
01157
01158
01159 ppnp.push_back(ptrk);
01160 }
01161 if(mdcTrk->charge() <0 && pid->probPion() > pid->probKaon()) {
01162 ipnm.push_back(jGood[i]);
01163 HepLorentzVector ptrk;
01164 ptrk.setPx(mdcKalTrk->px());
01165 ptrk.setPy(mdcKalTrk->py());
01166 ptrk.setPz(mdcKalTrk->pz());
01167 double p3 = ptrk.mag();
01168 ptrk.setE(sqrt(p3*p3+mpi*mpi));
01169
01170
01171
01172 ppnm.push_back(ptrk);
01173 }
01174 }
01175 int npnp = ipnp.size();
01176 int npnm = ipnm.size();
01177
01178 m_pnp = npnp;
01179 m_pnm = npnm;
01180
01181 int iphoton = 0;
01182 for (int i=0; i<m_nggneu; i++)
01183 {
01184 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + jGgam[i];
01185 if (!(*itTrk)->isEmcShowerValid()) continue;
01186 RecEmcShower *emcTrk = (*itTrk)->emcShower();
01187 m_numHits[iphoton] = emcTrk->numHits();
01188 m_secondmoment[iphoton] = emcTrk->secondMoment();
01189 m_x[iphoton] = emcTrk->x();
01190 m_y[iphoton]= emcTrk->y();
01191 m_z[iphoton]= emcTrk->z();
01192 m_cosemc[iphoton] = cos(emcTrk->theta());
01193 m_phiemc[iphoton] = emcTrk->phi();
01194 m_energy[iphoton] = emcTrk->energy();
01195 m_eSeed[iphoton] = emcTrk->eSeed();
01196 m_e3x3[iphoton] = emcTrk->e3x3();
01197 m_e5x5[iphoton] = emcTrk->e5x5();
01198 m_lat[iphoton] = emcTrk->latMoment();
01199 m_a20[iphoton] = emcTrk->a20Moment();
01200 m_a42[iphoton] = emcTrk->a42Moment();
01201 iphoton++;
01202 }
01203 m_tuple4->write();
01204 Ncut4++;
01205
01206 if(kmfit->chisq()>=chi5) return SUCCESS;
01207 if(pgam2pi2.m()<=1.0) return SUCCESS;
01208 if(pgam2pi1.m()<=1.0) return SUCCESS;
01209 if(nGam!=2) return SUCCESS;
01210 if(emcTg1/extmot1>=0.8) return SUCCESS;
01211 if(emcTg2/extmot2>=0.8) return SUCCESS;
01212 Ncut5++;
01213
01214
01215 TH1 *h(0);
01216 if (m_thsvc->getHist("/DQAHist/Rhopi/mpi0", h).isSuccess()) {
01217 h->Fill(ppi0.m());
01218 } else {
01219 log << MSG::ERROR << "Couldn't retrieve mpi0" << endreq;
01220 }
01221
01222 if(fabs(ppi0.m()-0.135)>=0.015) return SUCCESS;
01223 Ncut6++;
01224
01225 if (m_thsvc->getHist("/DQAHist/Rhopi/mrho0", h).isSuccess()) {
01226 h->Fill(prho0.m());
01227 } else {
01228 log << MSG::ERROR << "Couldn't retrieve mrho0" << endreq;
01229 }
01230
01231
01232 if (m_thsvc->getHist("/DQAHist/Rhopi/mrhop", h).isSuccess()) {
01233 h->Fill(prhop.m());
01234 } else {
01235 log << MSG::ERROR << "Couldn't retrieve mrhop" << endreq;
01236 }
01237
01238 if (m_thsvc->getHist("/DQAHist/Rhopi/mrhom", h).isSuccess()) {
01239 h->Fill(prhom.m());
01240 } else {
01241 log << MSG::ERROR << "Couldn't retrieve mrhom" << endreq;
01242 }
01243
01244
01245
01246
01247 (*(evtRecTrkCol->begin()+ipip[0]))->setPartId(3);
01248 (*(evtRecTrkCol->begin()+ipim[0]))->setPartId(3);
01249
01250
01251
01252
01253
01254
01255
01256 (*(evtRecTrkCol->begin()+ipip[0]))->setQuality(0);
01257 (*(evtRecTrkCol->begin()+ipim[0]))->setQuality(0);
01258
01259 setFilterPassed(true);
01260
01261 return StatusCode::SUCCESS;
01262 }
01263
01264
01265
01266 StatusCode DQARhopi::finalize() {
01267 cout<<"total number: "<<Ncut0<<endl;
01268 cout<<"nGood==2, nCharge==0: "<<Ncut1<<endl;
01269 cout<<"nGam>=2: "<<Ncut2<<endl;
01270 cout<<"Pass Pid: "<<Ncut3<<endl;
01271 cout<<"Pass 4C: "<<Ncut4<<endl;
01272 cout<<"final cut without pi0:"<<Ncut5<<endl;
01273 cout<<"final cut with pi0: "<<Ncut6<<endl;
01274 MsgStream log(msgSvc(), name());
01275 log << MSG::INFO << "in finalize()" << endmsg;
01276 return StatusCode::SUCCESS;
01277 }
01278