00001
00002
00003 #include "GaudiKernel/MsgStream.h"
00004 #include "GaudiKernel/AlgFactory.h"
00005 #include "GaudiKernel/ISvcLocator.h"
00006 #include "GaudiKernel/SmartDataPtr.h"
00007 #include "GaudiKernel/IDataProviderSvc.h"
00008 #include "GaudiKernel/PropertyMgr.h"
00009
00010 #include "EventModel/EventModel.h"
00011 #include "EventModel/EventHeader.h"
00012 #include "EventModel/Event.h"
00013 #include "TrigEvent/TrigEvent.h"
00014 #include "TrigEvent/TrigData.h"
00015 #include "McTruth/McParticle.h"
00016
00017 #include "EvtRecEvent/EvtRecEvent.h"
00018 #include "EvtRecEvent/EvtRecTrack.h"
00019 #include "DstEvent/TofHitStatus.h"
00020
00021 #include "TMath.h"
00022 #include "GaudiKernel/INTupleSvc.h"
00023 #include "GaudiKernel/NTuple.h"
00024 #include "GaudiKernel/Bootstrap.h"
00025 #include "GaudiKernel/IHistogramSvc.h"
00026 #include "CLHEP/Vector/ThreeVector.h"
00027 #include "CLHEP/Vector/LorentzVector.h"
00028 #include "CLHEP/Vector/TwoVector.h"
00029 using CLHEP::Hep3Vector;
00030 using CLHEP::Hep2Vector;
00031 using CLHEP::HepLorentzVector;
00032 #include "CLHEP/Geometry/Point3D.h"
00033 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00034 typedef HepGeom::Point3D<double> HepPoint3D;
00035 #endif
00036 #include "VertexFit/KinematicFit.h"
00037 #include "VertexFit/VertexFit.h"
00038 #include "VertexFit/SecondVertexFit.h"
00039 #include "VertexFit/WTrackParameter.h"
00040 #include "ParticleID/ParticleID.h"
00041 #include "MdcRecEvent/RecMdcKalTrack.h"
00042
00043 #include "PipiJpsiAlg/PipiJpsi.h"
00044
00045 #include <vector>
00046
00047
00048 const double me = 0.000511;
00049 const double mpi = 0.13957;
00050 const double mproton = 0.938272;
00051 const double mmu = 0.105658;
00052 const double mpsip = 3.686;
00053 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00054 const double velc = 29.9792458;
00055 const double PI = 3.1415926;
00056
00057 typedef std::vector<int> Vint;
00058 typedef std::vector<HepLorentzVector> Vp4;
00059
00060
00061 static long m_cout_all(0), m_cout_col(0), m_cout_charge(0), m_cout_nGood(0), m_cout_mom(0), m_cout_recoil(0), m_cout_everat(0);
00063
00064 PipiJpsi::PipiJpsi(const std::string& name, ISvcLocator* pSvcLocator) :
00065 Algorithm(name, pSvcLocator) {
00066
00067 declareProperty("Vr0cut", m_vr0cut=1.0);
00068 declareProperty("Vz0cut", m_vz0cut=5.0);
00069 declareProperty("TrackCosThetaCut",m_cosThetaCut=0.93);
00070 declareProperty("PipiDangCut",m_pipi_dang_cut=0.98);
00071
00072 declareProperty("CheckDedx", m_checkDedx = true);
00073 declareProperty("CheckTof", m_checkTof = true);
00074
00075 declareProperty("Subsample", m_subsample_flag=false);
00076 declareProperty("Trigger", m_trigger_flag=false);
00077 declareProperty("DistinEMuon", m_distin_emuon=2.0);
00078
00079 declareProperty("EventRate", m_eventrate=false);
00080 declareProperty("ChanDet", m_chan_det=1);
00081 }
00082
00083
00084 StatusCode PipiJpsi::initialize(){
00085 MsgStream log(msgSvc(), name());
00086
00087 log << MSG::INFO << "in initialize()" << endmsg;
00088
00089 StatusCode status;
00090
00091 NTuplePtr nt1(ntupleSvc(), "FILE1/vxyz");
00092 if ( nt1 ) m_tuple1 = nt1;
00093 else {
00094 m_tuple1 = ntupleSvc()->book ("FILE1/vxyz", CLID_ColumnWiseTuple, "ks N-Tuple example");
00095 if ( m_tuple1 ) {
00096 status = m_tuple1->addItem ("vx0", m_vx0);
00097 status = m_tuple1->addItem ("vy0", m_vy0);
00098 status = m_tuple1->addItem ("vz0", m_vz0);
00099 status = m_tuple1->addItem ("vr0", m_vr0);
00100 }
00101 else {
00102 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00103 return StatusCode::FAILURE;
00104 }
00105 }
00106
00107 NTuplePtr nt2(ntupleSvc(), "FILE1/photon");
00108 if ( nt2 ) m_tuple2 = nt2;
00109 else {
00110 m_tuple2 = ntupleSvc()->book ("FILE1/photon", CLID_ColumnWiseTuple, "ks N-Tuple example");
00111 if ( m_tuple2 ) {
00112 status = m_tuple2->addItem ("dthe", m_dthe);
00113 status = m_tuple2->addItem ("dphi", m_dphi);
00114 status = m_tuple2->addItem ("dang", m_dang);
00115 status = m_tuple2->addItem ("eraw", m_eraw);
00116 status = m_tuple2->addItem ("nGam", m_nGam);
00117 }
00118 else {
00119 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
00120 return StatusCode::FAILURE;
00121 }
00122 }
00123
00124
00125 if(m_checkDedx) {
00126 NTuplePtr nt3(ntupleSvc(), "FILE1/dedx");
00127 if ( nt3 ) m_tuple3 = nt3;
00128 else {
00129 m_tuple3 = ntupleSvc()->book ("FILE1/dedx", CLID_ColumnWiseTuple, "ks N-Tuple example");
00130 if ( m_tuple3 ) {
00131 status = m_tuple3->addItem ("ptrk", m_ptrk);
00132 status = m_tuple3->addItem ("chie", m_chie);
00133 status = m_tuple3->addItem ("chimu", m_chimu);
00134 status = m_tuple3->addItem ("chipi", m_chipi);
00135 status = m_tuple3->addItem ("chik", m_chik);
00136 status = m_tuple3->addItem ("chip", m_chip);
00137 status = m_tuple3->addItem ("probPH", m_probPH);
00138 status = m_tuple3->addItem ("normPH", m_normPH);
00139 status = m_tuple3->addItem ("ghit", m_ghit);
00140 status = m_tuple3->addItem ("thit", m_thit);
00141 }
00142 else {
00143 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple3) << endmsg;
00144 return StatusCode::FAILURE;
00145 }
00146 }
00147 }
00148
00149 if(m_checkTof) {
00150 NTuplePtr nt4(ntupleSvc(), "FILE1/tofe");
00151 if ( nt4 ) m_tuple4 = nt4;
00152 else {
00153 m_tuple4 = ntupleSvc()->book ("FILE1/tofe",CLID_ColumnWiseTuple, "ks N-Tuple example");
00154 if ( m_tuple4 ) {
00155 status = m_tuple4->addItem ("ptrk", m_ptot_etof);
00156 status = m_tuple4->addItem ("cntr", m_cntr_etof);
00157 status = m_tuple4->addItem ("path", m_path_etof);
00158 status = m_tuple4->addItem ("ph", m_ph_etof);
00159 status = m_tuple4->addItem ("rhit", m_rhit_etof);
00160 status = m_tuple4->addItem ("qual", m_qual_etof);
00161 status = m_tuple4->addItem ("tof", m_tof_etof);
00162 status = m_tuple4->addItem ("te", m_te_etof);
00163 status = m_tuple4->addItem ("tmu", m_tmu_etof);
00164 status = m_tuple4->addItem ("tpi", m_tpi_etof);
00165 status = m_tuple4->addItem ("tk", m_tk_etof);
00166 status = m_tuple4->addItem ("tp", m_tp_etof);
00167 }
00168 else {
00169 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple4) << endmsg;
00170 return StatusCode::FAILURE;
00171 }
00172 }
00173 }
00174
00175
00176
00177 if(m_checkTof) {
00178 NTuplePtr nt5(ntupleSvc(), "FILE1/tof1");
00179 if ( nt5 ) m_tuple5 = nt5;
00180 else {
00181 m_tuple5 = ntupleSvc()->book ("FILE1/tof1", CLID_ColumnWiseTuple, "ks N-Tuple example");
00182 if ( m_tuple5 ) {
00183 status = m_tuple5->addItem ("ptrk", m_ptot_btof1);
00184 status = m_tuple5->addItem ("cntr", m_cntr_btof1);
00185 status = m_tuple5->addItem ("path", m_path_btof1);
00186 status = m_tuple5->addItem ("ph", m_ph_btof1);
00187 status = m_tuple5->addItem ("zhit", m_zhit_btof1);
00188 status = m_tuple5->addItem ("qual", m_qual_btof1);
00189 status = m_tuple5->addItem ("tof", m_tof_btof1);
00190 status = m_tuple5->addItem ("te", m_te_btof1);
00191 status = m_tuple5->addItem ("tmu", m_tmu_btof1);
00192 status = m_tuple5->addItem ("tpi", m_tpi_btof1);
00193 status = m_tuple5->addItem ("tk", m_tk_btof1);
00194 status = m_tuple5->addItem ("tp", m_tp_btof1);
00195 }
00196 else {
00197 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple5) << endmsg;
00198 return StatusCode::FAILURE;
00199 }
00200 }
00201 }
00202
00203
00204 if(m_checkTof) {
00205 NTuplePtr nt6(ntupleSvc(), "FILE1/tof2");
00206 if ( nt6 ) m_tuple6 = nt6;
00207 else {
00208 m_tuple6 = ntupleSvc()->book ("FILE1/tof2", CLID_ColumnWiseTuple, "ks N-Tuple example");
00209 if ( m_tuple6 ) {
00210 status = m_tuple6->addItem ("ptrk", m_ptot_btof2);
00211 status = m_tuple6->addItem ("cntr", m_cntr_btof2);
00212 status = m_tuple6->addItem ("path", m_path_btof2);
00213 status = m_tuple6->addItem ("ph", m_ph_btof2);
00214 status = m_tuple6->addItem ("zhit", m_zhit_btof2);
00215 status = m_tuple6->addItem ("qual", m_qual_btof2);
00216 status = m_tuple6->addItem ("tof", m_tof_btof2);
00217 status = m_tuple6->addItem ("te", m_te_btof2);
00218 status = m_tuple6->addItem ("tmu", m_tmu_btof2);
00219 status = m_tuple6->addItem ("tpi", m_tpi_btof2);
00220 status = m_tuple6->addItem ("tk", m_tk_btof2);
00221 status = m_tuple6->addItem ("tp", m_tp_btof2);
00222 }
00223 else {
00224 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple6) << endmsg;
00225 return StatusCode::FAILURE;
00226 }
00227 }
00228 }
00229
00230 NTuplePtr nt8(ntupleSvc(), "FILE1/infmom");
00231 if ( nt8 ) m_tuple8 = nt8;
00232 else {
00233 m_tuple8 = ntupleSvc()->book ("FILE1/infmom", CLID_ColumnWiseTuple, "information with momentum method");
00234 if ( m_tuple8 ) {
00235 status = m_tuple8->addItem ("momlepp", m_mom_lepp );
00236 status = m_tuple8->addItem ("momlepmm", m_mom_lepm );
00237 status = m_tuple8->addItem ("mompionm", m_mom_pionm );
00238 status = m_tuple8->addItem ("mompionp", m_mom_pionp );
00239 status = m_tuple8->addItem ("pipidang", m_pipi_dang );
00240 status = m_tuple8->addItem ("cmslepp", m_cms_lepp );
00241 status = m_tuple8->addItem ("cmslepm", m_cms_lepm );
00242 status = m_tuple8->addItem ("invtwopi", m_mass_twopi );
00243 status = m_tuple8->addItem ("invjpsi", m_mass_jpsi );
00244 status = m_tuple8->addItem ("recoil", m_mass_recoil);
00245 status = m_tuple8->addItem ("invmass", m_inv_mass );
00246 status = m_tuple8->addItem ("totene", m_tot_e );
00247 status = m_tuple8->addItem ("totpx", m_tot_px );
00248 status = m_tuple8->addItem ("totpy", m_tot_py );
00249 status = m_tuple8->addItem ("totpz", m_tot_pz );
00250 status = m_tuple8->addItem ("epratio", m_ep_ratio );
00251 status = m_tuple8->addItem ("eveflag", m_event_flag );
00252 status = m_tuple8->addItem ("tplepratiom", m_trans_ratio_lepm );
00253 status = m_tuple8->addItem ("tplepratiop", m_trans_ratio_lepp );
00254 status = m_tuple8->addItem ("tppionratiom", m_trans_ratio_pionm );
00255 status = m_tuple8->addItem ("tppionratiop", m_trans_ratio_pionp );
00256 status = m_tuple8->addItem ("run", m_run );
00257 status = m_tuple8->addItem ("event", m_event );
00258 status = m_tuple8->addItem ("ntrack", m_index, 0, 4 );
00259 status = m_tuple8->addIndexedItem ("costhe", m_index, m_cos_theta );
00260 status = m_tuple8->addIndexedItem ("phi", m_index, m_phi );
00261 status = m_tuple8->addIndexedItem ("fourmom", m_index, 4, m_four_mom);
00262 status = m_tuple8->addItem ("pionmat", m_pion_matched);
00263 status = m_tuple8->addItem ("lepmat", m_lep_matched);
00264
00265 status = m_tuple8->addItem("indexmc", m_idxmc, 0, 100);
00266 status = m_tuple8->addIndexedItem("pdgid", m_idxmc, m_pdgid);
00267 status = m_tuple8->addIndexedItem("motheridx", m_idxmc, m_motheridx);
00268 status = m_tuple8->addItem("truepp", m_true_pionp);
00269 status = m_tuple8->addItem("truepm", m_true_pionm);
00270 }
00271 else {
00272 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple8) << endmsg;
00273 return StatusCode::FAILURE;
00274 }
00275 }
00276
00277
00278
00279
00280
00281 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00282 return StatusCode::SUCCESS;
00283
00284 }
00285
00286
00287 StatusCode PipiJpsi::execute() {
00288
00289
00290
00291 MsgStream log(msgSvc(), name());
00292 log << MSG::INFO << "in execute()" << endreq;
00293 m_cout_all ++;
00294 StatusCode sc=StatusCode::SUCCESS;
00295
00296 setFilterPassed(false);
00297
00298 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00299 if(!eventHeader){
00300 log << MSG::ERROR << "EventHeader not found" << endreq;
00301 return StatusCode::SUCCESS;
00302 }
00303 int run(eventHeader->runNumber());
00304 int event(eventHeader->eventNumber());
00305 if(event%1000==0) cout << "run: " << run << " event: " << event << endl;
00306
00307 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), "/Event/EvtRec/EvtRecEvent");
00308 if(!evtRecEvent){
00309 log << MSG::ERROR << "EvtRecEvent not found" << endreq;
00310 return StatusCode::SUCCESS;
00311 }
00312 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00313 << evtRecEvent->totalCharged() << " , "
00314 << evtRecEvent->totalNeutral() << " , "
00315 << evtRecEvent->totalTracks() <<endreq;
00316
00317 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), "/Event/EvtRec/EvtRecTrackCol");
00318 if(!evtRecTrkCol){
00319 log << MSG::ERROR << "EvtRecTrackCol" << endreq;
00320 return StatusCode::SUCCESS;
00321 }
00322
00323 if(m_trigger_flag){
00324 SmartDataPtr<TrigData> trigData(eventSvc(),EventModel::Trig::TrigData);
00325 if (!trigData) {
00326 log << MSG::FATAL << "Could not find Trigger Data for physics analysis" << endreq;
00327 return StatusCode::FAILURE;
00328 }
00330 log << MSG::DEBUG << "Trigger conditions: " << endreq;
00331 for(int i=0; i < 48; i++){
00332 log << MSG::DEBUG << "i:" << i << " name:" << trigData->getTrigCondName(i) << " cond:" << trigData->getTrigCondition(i) << endreq;
00333 }
00334
00335 int m_trig_tot(0), m_trig_which(-1);
00336 if(m_eventrate){
00337 for(int j=0; j<16; j++){
00338 if(trigData->getTrigChannel(j)){
00339 m_trig_tot ++;
00340 m_trig_which = j+1;
00341 }
00342 }
00343 if(m_trig_tot==1 && m_trig_which==m_chan_det) m_cout_everat++;
00344 return sc;
00345 }
00346 }
00347
00348 m_cout_col ++;
00349 if(evtRecEvent->totalCharged()<3 || evtRecTrkCol->size()<3 || evtRecEvent->totalTracks()>99 || evtRecTrkCol->size()>99) return StatusCode::SUCCESS;
00350 m_cout_charge ++;
00351
00352
00353 Vint iGood; iGood.clear();
00354 int m_num[4]={0,0,0,0};
00355 int nCharge = 0;
00356 m_pion_matched = 0; m_lep_matched = 0;
00357 HepLorentzVector m_lv_pionp, m_lv_pionm, m_lv_lepp, m_lv_lepm, m_lv_ele, m_lv_pos, m_lv_mum, m_lv_mup;
00358
00359 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00360 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00361 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00362 RecMdcKalTrack* mdcTrk = (*itTrk)->mdcKalTrack();
00363
00364 m_vx0 = mdcTrk->x();
00365 m_vy0 = mdcTrk->y();
00366 m_vz0 = mdcTrk->z();
00367 m_vr0 = mdcTrk->r();
00368 if(fabs(m_vz0) >= m_vz0cut) continue;
00369 if(m_vr0 >= m_vr0cut) continue;
00370 iGood.push_back(i);
00371 nCharge += mdcTrk->charge();
00372 if(mdcTrk->p()<1.0){ if((*itTrk)->isEmcShowerValid()) m_pion_matched ++; }
00373 else{ if((*itTrk)->isEmcShowerValid()) m_lep_matched ++; }
00374
00375 if(mdcTrk->charge()>0){
00376 if(mdcTrk->p()<1.0){
00377 mdcTrk->setPidType(RecMdcKalTrack::pion);
00378
00379 m_lv_pionp = mdcTrk->p4(xmass[2]);
00380 m_num[0] ++;
00381 }
00382 else{
00383 mdcTrk->setPidType(RecMdcKalTrack::electron);
00384 m_lv_pos = mdcTrk->p4(xmass[0]);
00385 mdcTrk->setPidType(RecMdcKalTrack::muon);
00386 m_lv_mup = mdcTrk->p4(xmass[1]);
00387 m_num[2] ++;
00388 }
00389 }
00390 else{
00391 if(mdcTrk->p()<1.0){
00392 mdcTrk->setPidType(RecMdcKalTrack::pion);
00393 m_lv_pionm = mdcTrk->p4(xmass[2]);
00394 m_num[1] ++;
00395 }
00396 else{
00397 mdcTrk->setPidType(RecMdcKalTrack::electron);
00398 m_lv_ele = mdcTrk->p4(xmass[0]);
00399 mdcTrk->setPidType(RecMdcKalTrack::muon);
00400 m_lv_mum = mdcTrk->p4(xmass[1]);
00401 m_num[3] ++;
00402 }
00403 }
00404 }
00405
00406 int nGood = iGood.size();
00407 log << MSG::DEBUG << "With KalmanTrack, ngood, totcharge = " << nGood << " , " << nCharge << endreq;
00408 if(nGood<3 || nGood>4) return sc;
00409 m_cout_nGood ++;
00410
00411 m_ep_ratio = 0;
00412 for(int i=0; i< evtRecEvent->totalTracks(); i++){
00413 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00414 if(!(*itTrk)->isEmcShowerValid()) continue;
00415 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00416 m_ep_ratio += emcTrk->energy();
00417 }
00418
00419 if(m_ep_ratio > m_distin_emuon){
00420 m_lv_lepp = m_lv_pos;
00421 m_lv_lepm = m_lv_ele;
00422 }
00423 else{
00424 m_lv_lepp = m_lv_mup;
00425 m_lv_lepm = m_lv_mum;
00426 }
00427
00428 HepLorentzVector m_lv_lab(0.04,0,0,3.686);
00429 if(nGood==4){
00430 if(nCharge) return sc;
00431 m_event_flag = 4;
00432 }
00433 else{
00434 if(m_num[0]>1 || m_num[1]>1 || m_num[2]>1 || m_num[3]>1) return sc;
00435 if(m_num[0]==0){
00436 if(nCharge != -1) return StatusCode::SUCCESS;
00437 m_lv_pionp = m_lv_lab - m_lv_pionm - m_lv_lepp - m_lv_lepm;
00438 if(m_lv_pionp.vect().cosTheta()>m_cosThetaCut)return StatusCode::SUCCESS;
00439 m_event_flag = 0;
00440 }
00441 if(m_num[1]==0){
00442 if(nCharge != 1) return StatusCode::SUCCESS;
00443 m_lv_pionm = m_lv_lab - m_lv_pionp - m_lv_lepp - m_lv_lepm;
00444 if(m_lv_pionm.vect().cosTheta()>m_cosThetaCut)return StatusCode::SUCCESS;
00445 m_event_flag = 1;
00446 }
00447 if(m_num[2]==0){
00448 if(nCharge != -1) return StatusCode::SUCCESS;
00449 m_lv_lepp = m_lv_lab - m_lv_pionp - m_lv_pionm - m_lv_lepm;
00450 if(m_lv_lepp.vect().cosTheta()>m_cosThetaCut)return StatusCode::SUCCESS;
00451 m_event_flag = 2;
00452 }
00453 if(m_num[3]==0){
00454 if(nCharge != 1) return StatusCode::SUCCESS;
00455 m_lv_lepm = m_lv_lab - m_lv_pionp - m_lv_pionm - m_lv_lepp;
00456 if(m_lv_lepm.vect().cosTheta()>m_cosThetaCut)return StatusCode::SUCCESS;
00457 m_event_flag = 3;
00458 }
00459 }
00460 m_cout_mom ++;
00461
00462
00463
00464 HepLorentzVector m_lv_recoil, m_lv_jpsi;
00465 m_lv_recoil = m_lv_lab - m_lv_pionp - m_lv_pionm;
00466 m_lv_jpsi = m_lv_lepp + m_lv_lepm;
00467
00468 m_mass_twopi = (m_lv_pionp + m_lv_pionm).m();
00469 m_mass_recoil = m_lv_recoil.m();
00470 m_mass_jpsi = m_lv_jpsi.m();
00471
00472
00473 if( m_mass_recoil < 3.05 || m_mass_recoil > 3.15 ) return sc;
00474 if( m_mass_jpsi < 3.0 || m_mass_jpsi > 3.2 ) return sc;
00475 m_cout_recoil ++;
00476
00477 HepLorentzVector m_ttm(m_lv_jpsi + m_lv_pionp + m_lv_pionm);
00478 if(m_ttm.m()>4 || m_ttm.m()<3) return sc;
00479
00480
00481 m_pipi_dang = m_lv_pionp.vect().cosTheta(m_lv_pionm.vect());
00482
00483 m_mom_pionp = m_lv_pionp.vect().mag();
00484 m_mom_pionm = m_lv_pionm.vect().mag();
00485 m_mom_lepp = m_lv_lepp.vect().mag();
00486 m_mom_lepm = m_lv_lepm.vect().mag();
00487 m_trans_ratio_lepp = m_lv_lepp.vect().perp()/m_lv_lepp.vect().mag();
00488 m_trans_ratio_lepm = m_lv_lepm.vect().perp()/m_lv_lepm.vect().mag();
00489 m_trans_ratio_pionp = m_lv_pionp.vect().perp()/m_lv_pionp.vect().mag();
00490 m_trans_ratio_pionm = m_lv_pionm.vect().perp()/m_lv_pionm.vect().mag();
00491
00492 Hep3Vector m_boost_jpsi(m_lv_recoil.boostVector());
00493 HepLorentzVector m_lv_cms_lepp(boostOf(m_lv_lepp,-m_boost_jpsi));
00494 HepLorentzVector m_lv_cms_lepm(boostOf(m_lv_lepm,-m_boost_jpsi));
00495 m_cms_lepm = m_lv_cms_lepm.vect().mag();
00496 m_cms_lepp = m_lv_cms_lepp.vect().mag();
00497 log << MSG::DEBUG << "jpsi four momentum in cms " << m_lv_cms_lepp + m_lv_cms_lepm <<endreq;
00498
00499 m_inv_mass = m_ttm.m();
00500 m_tot_e = m_ttm.e();
00501 m_tot_px = m_ttm.px();
00502 m_tot_py = m_ttm.py();
00503 m_tot_pz = m_ttm.pz();
00504 m_run = run;
00505 m_event = event;
00506 HepLorentzVector m_lv_book(0,0,0,0);
00507 for(m_index=0; m_index<4; m_index++){
00508 switch(m_index){
00509 case 0: m_lv_book = m_lv_pionp;
00510 break;
00511 case 1: m_lv_book = m_lv_pionm;
00512 break;
00513 case 2: m_lv_book = m_lv_lepp;
00514 break;
00515 case 3: m_lv_book = m_lv_lepm;
00516 break;
00517 default: m_lv_book.setE(2008);
00518 }
00519 m_cos_theta[m_index] = m_lv_book.vect().cosTheta();
00520 m_phi[m_index] = m_lv_book.vect().phi();
00521 m_four_mom[m_index][0] = m_lv_book.e();
00522 m_four_mom[m_index][1] = m_lv_book.px();
00523 m_four_mom[m_index][2] = m_lv_book.py();
00524 m_four_mom[m_index][3] = m_lv_book.pz();
00525 }
00526
00527 if(m_subsample_flag) setFilterPassed(true);
00528 else if(m_mass_recoil>3.08 && m_mass_recoil<3.12 && m_mass_jpsi>3.0 && m_mass_jpsi<3.2 && m_cms_lepp<1.7 && m_cms_lepp>1.4 && m_cms_lepm<1.7 && m_cms_lepm>1.4 && m_event_flag==4 && m_pipi_dang<m_pipi_dang_cut) setFilterPassed(true);
00529
00530
00531
00532 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
00533 if(m_run<0){
00534 int m_numParticle(0), m_true_pid(0);
00535 if(!mcParticleCol){
00536 log << MSG::ERROR << "Could not retrieve McParticelCol" << endreq;
00537 return StatusCode::FAILURE;
00538 }
00539 else{
00540 bool psipDecay(false);
00541 int rootIndex(-1);
00542 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
00543 for (; iter_mc != mcParticleCol->end(); iter_mc++){
00544 if ((*iter_mc)->primaryParticle()) continue;
00545 if (!(*iter_mc)->decayFromGenerator()) continue;
00546
00547 if ((*iter_mc)->particleProperty()==100443){
00548 psipDecay = true;
00549 rootIndex = (*iter_mc)->trackIndex();
00550 }
00551 if (!psipDecay) continue;
00552 int mcidx = ((*iter_mc)->mother()).trackIndex() - rootIndex;
00553 int pdgid = (*iter_mc)->particleProperty();
00554 m_pdgid[m_numParticle] = pdgid;
00555 m_motheridx[m_numParticle] = mcidx;
00556 m_numParticle ++;
00557
00558
00559 if((*iter_mc)->particleProperty() == 211) m_true_pionp = (*iter_mc)->initialFourMomentum().vect().mag();
00560 if((*iter_mc)->particleProperty() == -211) m_true_pionm = (*iter_mc)->initialFourMomentum().vect().mag();
00561 }
00562 m_idxmc = m_numParticle;
00563 }
00564 }
00565
00566 m_tuple1->write();
00567 m_tuple8->write();
00568
00569
00570
00571 Vint iGam; iGam.clear();
00572 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
00573 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00574 if(!(*itTrk)->isEmcShowerValid()) continue;
00575 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00576 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
00577
00578 double dthe = 200.;
00579 double dphi = 200.;
00580 double dang = 200.;
00581 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
00582 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
00583 if(!(*jtTrk)->isExtTrackValid()) continue;
00584 RecExtTrack *extTrk = (*jtTrk)->extTrack();
00585 if(extTrk->emcVolumeNumber() == -1) continue;
00586 Hep3Vector extpos = extTrk->emcPosition();
00587
00588 double angd = extpos.angle(emcpos);
00589 double thed = extpos.theta() - emcpos.theta();
00590 double phid = extpos.deltaPhi(emcpos);
00591 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00592 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00593
00594 if(fabs(thed) < fabs(dthe)) dthe = thed;
00595 if(fabs(phid) < fabs(dphi)) dphi = phid;
00596 if(angd < dang) dang = angd;
00597 }
00598 if(dang>=200) continue;
00599 double eraw = emcTrk->energy();
00600 dthe = dthe * 180 / (CLHEP::pi);
00601 dphi = dphi * 180 / (CLHEP::pi);
00602 dang = dang * 180 / (CLHEP::pi);
00603 m_dthe = dthe;
00604 m_dphi = dphi;
00605 m_dang = dang;
00606 m_eraw = eraw;
00607
00608
00609
00610 iGam.push_back((*itTrk)->trackId());
00611 }
00612
00613 m_nGam = iGam.size();
00614 log << MSG::DEBUG << "num Good Photon " << m_nGam << " , " <<evtRecEvent->totalNeutral()<<endreq;
00615 m_tuple2->write();
00616
00617
00618
00619
00620 if(m_checkDedx) {
00621 int m_dedx_cout(0);
00622 for(int i = 0; i < nGood; i++) {
00623 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00624 if(!(*itTrk)->isMdcDedxValid())continue;
00625 RecMdcKalTrack *mdcTrk = (*itTrk)->mdcKalTrack();
00626 RecMdcDedx *dedxTrk = (*itTrk)->mdcDedx();
00627
00628 m_ptrk = mdcTrk->p();
00629 m_chie = dedxTrk->chiE();
00630 m_chimu = dedxTrk->chiMu();
00631 m_chipi = dedxTrk->chiPi();
00632 m_chik = dedxTrk->chiK();
00633 m_chip = dedxTrk->chiP();
00634 m_ghit = dedxTrk->numGoodHits();
00635 m_thit = dedxTrk->numTotalHits();
00636 m_probPH = dedxTrk->probPH();
00637 m_normPH = dedxTrk->normPH();
00638
00639 m_tuple3->write();
00640 }
00641 }
00642
00643
00644
00645
00646 if(m_checkTof) {
00647 int m_endcap_cout(0), m_layer1_cout(0), m_layer2_cout(0);
00648 for(int i = 0; i < nGood; i++) {
00649 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00650 if(!(*itTrk)->isTofTrackValid()) continue;
00651
00652 RecMdcKalTrack *mdcTrk = (*itTrk)->mdcKalTrack();
00653 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
00654
00655 double ptrk = mdcTrk->p();
00656
00657 for( SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin() ;iter_tof != tofTrkCol.end(); iter_tof++ ) {
00658 TofHitStatus *status = new TofHitStatus;
00659 status->setStatus((*iter_tof)->status());
00660 if(!(status->is_barrel())){
00661 if( !(status->is_counter()) ) continue;
00662 if( status->layer()!=0 ) continue;
00663 double path = (*iter_tof)->path();
00664 double tof = (*iter_tof)->tof();
00665 double ph = (*iter_tof)->ph();
00666 double rhit = (*iter_tof)->zrhit();
00667 double qual = 0.0 + (*iter_tof)->quality();
00668 double cntr = 0.0 + (*iter_tof)->tofID();
00669 double texp[5];
00670 for(int j = 0; j < 5; j++) {
00671 double gb = xmass[j]/ptrk;
00672 double beta = sqrt(1+gb*gb);
00673 texp[j] = path*beta/velc;
00674 }
00675 m_cntr_etof = cntr;
00676 m_ptot_etof = ptrk;
00677 m_path_etof = path;
00678 m_ph_etof = ph;
00679 m_rhit_etof = rhit;
00680 m_qual_etof = qual;
00681 m_tof_etof = tof;
00682 m_te_etof = tof - texp[0];
00683 m_tmu_etof = tof - texp[1];
00684 m_tpi_etof = tof - texp[2];
00685 m_tk_etof = tof - texp[3];
00686 m_tp_etof = tof - texp[4];
00687
00688 m_tuple4->write();
00689 }
00690 else {
00691 if( !(status->is_counter()) ) continue;
00692 if(status->layer()==1){
00693 double path=(*iter_tof)->path();
00694 double tof = (*iter_tof)->tof();
00695 double ph = (*iter_tof)->ph();
00696 double rhit = (*iter_tof)->zrhit();
00697 double qual = 0.0 + (*iter_tof)->quality();
00698 double cntr = 0.0 + (*iter_tof)->tofID();
00699 double texp[5];
00700 for(int j = 0; j < 5; j++) {
00701 double gb = xmass[j]/ptrk;
00702 double beta = sqrt(1+gb*gb);
00703 texp[j] = path*beta/velc;
00704 }
00705
00706 m_cntr_btof1 = cntr;
00707 m_ptot_btof1 = ptrk;
00708 m_path_btof1 = path;
00709 m_ph_btof1 = ph;
00710 m_zhit_btof1 = rhit;
00711 m_qual_btof1 = qual;
00712 m_tof_btof1 = tof;
00713 m_te_btof1 = tof - texp[0];
00714 m_tmu_btof1 = tof - texp[1];
00715 m_tpi_btof1 = tof - texp[2];
00716 m_tk_btof1 = tof - texp[3];
00717 m_tp_btof1 = tof - texp[4];
00718
00719 m_tuple5->write();
00720 }
00721
00722 if(status->layer()==2){
00723 double path=(*iter_tof)->path();
00724 double tof = (*iter_tof)->tof();
00725 double ph = (*iter_tof)->ph();
00726 double rhit = (*iter_tof)->zrhit();
00727 double qual = 0.0 + (*iter_tof)->quality();
00728 double cntr = 0.0 + (*iter_tof)->tofID();
00729 double texp[5];
00730 for(int j = 0; j < 5; j++) {
00731 double gb = xmass[j]/ptrk;
00732 double beta = sqrt(1+gb*gb);
00733 texp[j] = path*beta/velc;
00734 }
00735
00736 m_cntr_btof2 = cntr;
00737 m_ptot_btof2 = ptrk;
00738 m_path_btof2 = path;
00739 m_ph_btof2 = ph;
00740 m_zhit_btof2 = rhit;
00741 m_qual_btof2 = qual;
00742 m_tof_btof2 = tof;
00743 m_te_btof2 = tof - texp[0];
00744 m_tmu_btof2 = tof - texp[1];
00745 m_tpi_btof2 = tof - texp[2];
00746 m_tk_btof2 = tof - texp[3];
00747 m_tp_btof2 = tof - texp[4];
00748
00749 m_tuple6->write();
00750 }
00751 }
00752
00753 delete status;
00754 }
00755 }
00756 }
00757
00758
00759 return StatusCode::SUCCESS;
00760 }
00761
00762
00763 StatusCode PipiJpsi::finalize() {
00764
00765 MsgStream log(msgSvc(), name());
00766 log << MSG::INFO << "in finalize()" << endmsg;
00767 if(m_eventrate) cout << "all event: " << m_cout_all << endl << "only channel " << m_chan_det << ": " << m_cout_everat << endl;
00768 cout << "all event: " << m_cout_all << endl << "all collection point is OK: " << m_cout_col << endl << "charged tracks >=3: " << m_cout_charge << endl << "good charged tracks [3,4]: " << m_cout_nGood << endl << "after momentum assign: " << m_cout_mom << endl << "after recoild mass cut: " << m_cout_recoil << endl;
00769 return StatusCode::SUCCESS;
00770 }
00771
00772