00001 #include "GaudiKernel/MsgStream.h"
00002 #include "GaudiKernel/AlgFactory.h"
00003 #include "GaudiKernel/SmartDataPtr.h"
00004 #include "GaudiKernel/IDataProviderSvc.h"
00005 #include "GaudiKernel/PropertyMgr.h"
00006
00007 #include "EventModel/EventModel.h"
00008 #include "EventModel/Event.h"
00009 #include "EventModel/EventHeader.h"
00010 #include "EvtRecEvent/EvtRecEvent.h"
00011 #include "EvtRecEvent/EvtRecTrack.h"
00012
00013 #include "TMath.h"
00014 #include "GaudiKernel/INTupleSvc.h"
00015 #include "GaudiKernel/NTuple.h"
00016 #include "GaudiKernel/Bootstrap.h"
00017 #include "GaudiKernel/IHistogramSvc.h"
00018 #include "CLHEP/Vector/ThreeVector.h"
00019 #include "CLHEP/Vector/LorentzVector.h"
00020 #include "CLHEP/Vector/TwoVector.h"
00021 #include "EmcRawEvent/EmcDigi.h"
00022 #include "RawEvent/RawDataUtil.h"
00023 #include "MdcRawEvent/MdcDigi.h"
00024
00025 #include "GaudiKernel/Bootstrap.h"
00026 #include "GaudiKernel/ISvcLocator.h"
00027 using CLHEP::Hep3Vector;
00028 using CLHEP::Hep2Vector;
00029 using CLHEP::HepLorentzVector;
00030 #include "CLHEP/Geometry/Point3D.h"
00031 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00032 typedef HepGeom::Point3D<double> HepPoint3D;
00033 #endif
00034 #include "EventPreSelect/EventPreSelect.h"
00035 #include "EventPreSelect/DimuPreSelect.h"
00036
00037 #include <vector>
00038
00039 typedef std::vector<int> Vint;
00040 typedef std::vector<HepLorentzVector> Vp4;
00041
00042 const double mpi = 0.13957;
00043
00044 int EventPreSelect::idmax[43]={40,44,48,56,64,72,80,80,76,76,
00045 88,88,100,100,112,112,128,128,140,140,
00046 160,160,160,160,176,176,176,176,208,208,
00047 208,208,240,240,240,240,256,256,256,256,
00048 288,288,288};
00050
00051 EventPreSelect::EventPreSelect(const std::string& name, ISvcLocator* pSvcLocator) :
00052 Algorithm(name, pSvcLocator) {
00053
00054
00055 declareProperty("Output", m_output = false);
00056 declareProperty("Ecm", m_ecm=3.686);
00057 declareProperty("SelectBhabha", m_selectBhabha=false);
00058 declareProperty("SelectDimu", m_selectDimu=false);
00059 declareProperty("SelectHadron", m_selectHadron=false);
00060 declareProperty("SelectDiphoton", m_selectDiphoton=false);
00061 declareProperty("WriteDst", m_writeDst=false);
00062 declareProperty("WriteRec", m_writeRec=false);
00063 declareProperty("Vr0cut", m_vr0cut=1.0);
00064 declareProperty("Vz0cut", m_vz0cut=5.0);
00065 declareProperty("Pt0HighCut", m_pt0HighCut=5.0);
00066 declareProperty("Pt0LowCut", m_pt0LowCut=0.05);
00067 declareProperty("EnergyThreshold", m_energyThreshold=0.05);
00068 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
00069 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
00070
00071 declareProperty("BhabhaEmcECut", m_bhabhaEmcECut=0.7*m_ecm);
00072 declareProperty("BhabhaMaxECut", m_bhabhaMaxECut=0.3*m_ecm);
00073 declareProperty("BhabhaSecECut", m_bhabhaSecECut=0.1*m_ecm);
00074 declareProperty("BhabhaDTheCut", m_bhabhaDTheCut=3.);
00075 declareProperty("BhabhaDPhiCut1", m_bhabhaDPhiCut1=-25.);
00076 declareProperty("BhabhaDPhiCut2", m_bhabhaDPhiCut2=-4.);
00077 declareProperty("BhabhaDPhiCut3", m_bhabhaDPhiCut3=2.);
00078 declareProperty("BhabhaDPhiCut4", m_bhabhaDPhiCut4=20.);
00079 declareProperty("BhabhaMdcHitCutB", m_bhabhaMdcHitCutB=15);
00080 declareProperty("BhabhaMdcHitCutE", m_bhabhaMdcHitCutE=5);
00081
00082 declareProperty("DimuEHighCut", m_dimuEHighCut=0.27*m_ecm);
00083 declareProperty("DimuELowCut", m_dimuELowCut=0.027*m_ecm);
00084 declareProperty("DimuDTheCut", m_dimuDTheCut=3.);
00085 declareProperty("DimuDPhiCut", m_dimuDPhiCut=23.);
00086
00087 declareProperty("HadronChaECut", m_hadronChaECut=0.3*m_ecm);
00088 declareProperty("HadronNeuECut", m_hadronNeuECut=0.3*m_ecm);
00089
00090 declareProperty("DiphotonEmcECut", m_diphotonEmcECut=0.7*m_ecm);
00091 declareProperty("DiphotonSecECut", m_diphotonSecECut=0.3*m_ecm);
00092 declareProperty("DiphotonDTheCut", m_diphotonDTheCut=3.);
00093 declareProperty("DiphotonDPhiCut1", m_diphotonDPhiCut1=-4.);
00094 declareProperty("DiphotonDPhiCut2", m_diphotonDPhiCut2=2.);
00095 }
00096
00097
00098 StatusCode EventPreSelect::initialize() {
00099 MsgStream log(msgSvc(), name());
00100
00101 log << MSG::INFO << "in initialize()" << endmsg;
00102
00103 m_bhabhaEmcECut=0.7*m_ecm;
00104 m_bhabhaMaxECut=0.3*m_ecm;
00105 m_bhabhaSecECut=0.1*m_ecm;
00106 m_dimuEHighCut=0.27*m_ecm;
00107 m_dimuELowCut=0.027*m_ecm;
00108 m_diphotonEmcECut=0.7*m_ecm;
00109 m_diphotonSecECut=0.3*m_ecm;
00110 m_hadronChaECut=0.3*m_ecm;
00111 m_hadronNeuECut=0.3*m_ecm;
00112
00113
00114 if (m_ecm==2.4) {
00115 m_bhabhaEmcECut=0.7*m_ecm;
00116 m_bhabhaMaxECut=0.15*m_ecm;
00117 m_bhabhaSecECut=0.05*m_ecm;
00118 m_diphotonEmcECut=0.7*m_ecm;
00119 m_diphotonSecECut=0.2*m_ecm;
00120
00121 }
00122
00123
00124
00125 StatusCode sc;
00126
00127 log << MSG::INFO << "creating sub-algorithms...." << endreq;
00128
00129 if(m_selectBhabha) {
00130 sc = createSubAlgorithm( "EventWriter", "SelectBarrelBhabha", m_subalg1);
00131 if( sc.isFailure() ) {
00132 log << MSG::ERROR << "Error creating Sub-Algorithm SelectBhabhaBarrel" <<endreq;
00133 return sc;
00134 } else {
00135 log << MSG::INFO << "Success creating Sub-Algorithm SelectBhabhaBarrel" <<endreq;
00136 }
00137 sc = createSubAlgorithm( "EventWriter", "SelectEndcapBhabha", m_subalg2);
00138 if( sc.isFailure() ) {
00139 log << MSG::ERROR << "Error creating Sub-Algorithm SelectBhabhaEndcap" <<endreq;
00140 return sc;
00141 } else {
00142 log << MSG::INFO << "Success creating Sub-Algorithm SelectBhabhaEndcap" <<endreq;
00143 }
00144 }
00145
00146 if(m_selectDimu) {
00147 m_dimuAlg = new DimuPreSelect;
00148 sc = createSubAlgorithm( "EventWriter", "SelectBarrelDimu", m_subalg3);
00149 if( sc.isFailure() ) {
00150 log << MSG::ERROR << "Error creating Sub-Algorithm SelectDimuBarrel" <<endreq;
00151 return sc;
00152 } else {
00153 log << MSG::INFO << "Success creating Sub-Algorithm SelectDimuBarrel" <<endreq;
00154 }
00155 sc = createSubAlgorithm( "EventWriter", "SelectEndcapDimu", m_subalg4);
00156 if( sc.isFailure() ) {
00157 log << MSG::ERROR << "Error creating Sub-Algorithm SelectDimuEndcap" <<endreq;
00158 return sc;
00159 } else {
00160 log << MSG::INFO << "Success creating Sub-Algorithm SelectDimuEndap" <<endreq;
00161 }
00162 }
00163
00164 if(m_selectHadron) {
00165 sc = createSubAlgorithm( "EventWriter", "SelectHadron", m_subalg5);
00166 if( sc.isFailure() ) {
00167 log << MSG::ERROR << "Error creating Sub-Algorithm SelectHadron" <<endreq;
00168 return sc;
00169 } else {
00170 log << MSG::INFO << "Success creating Sub-Algorithm SelectHadron" <<endreq;
00171 }
00172 }
00173
00174 if(m_selectDiphoton) {
00175 sc = createSubAlgorithm( "EventWriter", "SelectBarrelDiphoton", m_subalg6);
00176 if( sc.isFailure() ) {
00177 log << MSG::ERROR << "Error creating Sub-Algorithm SelectDiphotonBarrel" <<endreq;
00178 return sc;
00179 } else {
00180 log << MSG::INFO << "Success creating Sub-Algorithm SelectDiphotonBarrel" <<endreq;
00181 }
00182 sc = createSubAlgorithm( "EventWriter", "SelectEndcapDiphoton", m_subalg7);
00183 if( sc.isFailure() ) {
00184 log << MSG::ERROR << "Error creating Sub-Algorithm SelectDiphotonEndcap" <<endreq;
00185 return sc;
00186 } else {
00187 log << MSG::INFO << "Success creating Sub-Algorithm SelectDiphotonEndcap" <<endreq;
00188 }
00189 }
00190
00191 if(m_writeDst) {
00192 sc = createSubAlgorithm( "EventWriter", "WriteDst", m_subalg8);
00193 if( sc.isFailure() ) {
00194 log << MSG::ERROR << "Error creating Sub-Algorithm WriteDst" <<endreq;
00195 return sc;
00196 } else {
00197 log << MSG::INFO << "Success creating Sub-Algorithm WriteDst" <<endreq;
00198 }
00199 }
00200
00201 if(m_writeRec) {
00202 sc = createSubAlgorithm( "EventWriter", "WriteRec", m_subalg9);
00203 if( sc.isFailure() ) {
00204 log << MSG::ERROR << "Error creating Sub-Algorithm WriteRec" <<endreq;
00205 return sc;
00206 } else {
00207 log << MSG::INFO << "Success creating Sub-Algorithm WriteRec" <<endreq;
00208 }
00209 }
00210
00211 if(m_output) {
00212 StatusCode status;
00213 NTuplePtr nt0(ntupleSvc(), "FILE1/hadron");
00214 if ( nt0 ) m_tuple0 = nt0;
00215 else {
00216 m_tuple0 = ntupleSvc()->book ("FILE1/hadron", CLID_ColumnWiseTuple, "N-Tuple example");
00217 if ( m_tuple0 ) {
00218 status = m_tuple0->addItem ("esum", m_esum);
00219 status = m_tuple0->addItem ("eemc", m_eemc);
00220 status = m_tuple0->addItem ("etot", m_etot);
00221 status = m_tuple0->addItem ("nGood", m_nGood);
00222 status = m_tuple0->addItem ("nCharge", m_nCharge);
00223 status = m_tuple0->addItem ("nGam", m_nGam);
00224 status = m_tuple0->addItem ("ptot", m_ptot);
00225 status = m_tuple0->addItem ("pp", m_pp);
00226 status = m_tuple0->addItem ("pm", m_pm);
00227 status = m_tuple0->addItem ("run", m_runnb);
00228 status = m_tuple0->addItem ("event", m_evtnb);
00229 status = m_tuple0->addItem ("maxE", m_maxE);
00230 status = m_tuple0->addItem ("secE", m_secE);
00231 status = m_tuple0->addItem ("dthe", m_dThe);
00232 status = m_tuple0->addItem ("dphi", m_dPhi);
00233 status = m_tuple0->addItem ("mdcHit1", m_mdcHit1);
00234 status = m_tuple0->addItem ("mdcHit2", m_mdcHit2);
00235 }
00236 else {
00237 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple0) << endmsg;
00238 return StatusCode::FAILURE;
00239 }
00240 }
00241
00242 NTuplePtr nt1(ntupleSvc(), "FILE1/vxyz");
00243 if ( nt1 ) m_tuple1 = nt1;
00244 else {
00245 m_tuple1 = ntupleSvc()->book ("FILE1/vxyz", CLID_ColumnWiseTuple, "ks N-Tuple example");
00246 if ( m_tuple1 ) {
00247 status = m_tuple1->addItem ("vx0", m_vx0);
00248 status = m_tuple1->addItem ("vy0", m_vy0);
00249 status = m_tuple1->addItem ("vz0", m_vz0);
00250 status = m_tuple1->addItem ("vr0", m_vr0);
00251 status = m_tuple1->addItem ("theta0", m_theta0);
00252 status = m_tuple1->addItem ("p0", m_p0);
00253 status = m_tuple1->addItem ("pt0", m_pt0);
00254 }
00255 else {
00256 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00257 return StatusCode::FAILURE;
00258 }
00259 }
00260
00261 NTuplePtr nt2(ntupleSvc(), "FILE1/photon");
00262 if ( nt2 ) m_tuple2 = nt2;
00263 else {
00264 m_tuple2 = ntupleSvc()->book ("FILE1/photon", CLID_ColumnWiseTuple, "ks N-Tuple example");
00265 if ( m_tuple2 ) {
00266 status = m_tuple2->addItem ("dthe", m_dthe);
00267 status = m_tuple2->addItem ("dphi", m_dphi);
00268 status = m_tuple2->addItem ("dang", m_dang);
00269 status = m_tuple2->addItem ("eraw", m_eraw);
00270 }
00271 else {
00272 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
00273 return StatusCode::FAILURE;
00274 }
00275 }
00276
00277 if(m_selectDimu) {
00278 NTuplePtr nt3(ntupleSvc(), "FILE1/dimu");
00279 if ( nt3 ) m_tuple3 = nt3;
00280 else {
00281 m_tuple3 = ntupleSvc()->book ("FILE1/dimu", CLID_ColumnWiseTuple, "ks N-Tuple example");
00282 if ( m_tuple3 ) {
00283 m_dimuAlg->BookNtuple(m_tuple3);
00284 }
00285 else {
00286 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple3) << endmsg;
00287 return StatusCode::FAILURE;
00288 }
00289 }
00290 }
00291 }
00292
00293
00294
00295
00296 m_events=0;
00297 m_barrelBhabhaNumber=0;
00298 m_endcapBhabhaNumber=0;
00299 m_barrelDimuNumber=0;
00300 m_endcapDimuNumber=0;
00301 m_hadronNumber=0;
00302 m_barrelDiphotonNumber=0;
00303 m_endcapDiphotonNumber=0;
00304
00305 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00306 return StatusCode::SUCCESS;
00307
00308 }
00309
00310
00311 StatusCode EventPreSelect::execute() {
00312
00313
00314
00315 MsgStream log(msgSvc(), name());
00316 log << MSG::INFO << "in execute()" << endreq;
00317
00318 if( m_writeDst) m_subalg8->execute();
00319 if( m_writeRec) m_subalg9->execute();
00320
00321 m_isBarrelBhabha = false;
00322 m_isEndcapBhabha = false;
00323 m_isBarrelDimu = false;
00324 m_isEndcapDimu = false;
00325 m_isHadron = false;
00326 m_isBarrelDiphoton = false;
00327 m_isEndcapDiphoton = false;
00328
00329 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00330 if(!eventHeader)
00331 {
00332 cout<<" eventHeader "<<endl;
00333 return StatusCode::FAILURE;
00334 }
00335
00336 int run=eventHeader->runNumber();
00337 int event=eventHeader->eventNumber();
00338
00339 if(m_events%1000==0) cout<< m_events << " -------- run,event: "<<run<<","<<event<<endl;
00340 m_events++;
00341
00342 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00343 if(!evtRecEvent ) {
00344 cout<<" evtRecEvent "<<endl;
00345 return StatusCode::FAILURE;
00346 }
00347
00348 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00349 << evtRecEvent->totalCharged() << " , "
00350 << evtRecEvent->totalNeutral() << " , "
00351 << evtRecEvent->totalTracks() <<endreq;
00352 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00353 if(!evtRecTrkCol){
00354 cout<<" evtRecTrkCol "<<endl;
00355 return StatusCode::FAILURE;
00356 }
00357
00358 if(evtRecEvent->totalTracks()!=evtRecTrkCol->size()) return StatusCode::SUCCESS;
00359
00360
00361 Vint iGood;
00362 iGood.clear();
00363 int nCharge = 0;
00364 for(int i = 0;i < evtRecEvent->totalCharged(); i++)
00365 {
00366
00367 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00368 if(!(*itTrk)->isMdcTrackValid()) continue;
00369 RecMdcTrack *mdcTrk =(*itTrk)->mdcTrack();
00370 double vx0 = mdcTrk->x();
00371 double vy0 = mdcTrk->y();
00372 double vz0 = mdcTrk->z();
00373 double vr0 = mdcTrk->r();
00374 double theta0 = mdcTrk->theta();
00375 double p0 = mdcTrk->p();
00376 double pt0 = mdcTrk->pxy();
00377
00378 if(m_output) {
00379 m_vx0 = vx0;
00380 m_vy0 = vy0;
00381 m_vz0 = vz0;
00382 m_vr0 = vr0;
00383 m_theta0 = theta0;
00384 m_p0 = p0;
00385 m_pt0 = pt0;
00386 m_tuple1->write();
00387 }
00388
00389 if(fabs(vz0) >= m_vz0cut) continue;
00390 if(vr0 >= m_vr0cut) continue;
00391 if(pt0 >= m_pt0HighCut) continue;
00392 if(pt0 <= m_pt0LowCut) continue;
00393
00394 iGood.push_back((*itTrk)->trackId());
00395 nCharge += mdcTrk->charge();
00396 }
00397 int nGood = iGood.size();
00398
00399
00400 Vint iGam;
00401 iGam.clear();
00402 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
00403
00404 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00405 if(!(*itTrk)->isEmcShowerValid()) continue;
00406 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00407 double eraw = emcTrk->energy();
00408 if(m_output) {
00409 m_eraw = eraw;
00410 m_tuple2->write();
00411 }
00412 if(eraw < m_energyThreshold) continue;
00413 iGam.push_back((*itTrk)->trackId());
00414 }
00415 int nGam = iGam.size();
00416
00417
00418 Vint ipip, ipim;
00419 ipip.clear();
00420 ipim.clear();
00421 Vp4 ppip, ppim;
00422 ppip.clear();
00423 ppim.clear();
00424
00425
00426 double echarge = 0.;
00427 double ptot = 0.;
00428 double etot = 0.;
00429 double eemc = 0.;
00430 double pp = 0.;
00431 double pm = 0.;
00432
00433 for(int i = 0; i < nGood; i++) {
00434 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00435
00436 if((*itTrk)->isMdcTrackValid()) {
00437
00438
00439
00440
00441
00442 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00443
00444 ptot += mdcTrk->p();
00445
00446 HepLorentzVector ptrk;
00447 ptrk.setPx(mdcTrk->px());
00448 ptrk.setPy(mdcTrk->py());
00449 ptrk.setPz(mdcTrk->pz());
00450 double p3 = ptrk.mag();
00451 ptrk.setE(sqrt(p3*p3+mpi*mpi));
00452 ptrk = ptrk.boost(-0.011,0,0);
00453
00454 echarge += ptrk.e();
00455 etot += ptrk.e();
00456
00457 if(mdcTrk->charge() >0) {
00458 ppip.push_back(ptrk);
00459 pp = ptrk.rho();
00460 } else {
00461 ppim.push_back(ptrk);
00462 pm = ptrk.rho();
00463 }
00464
00465 }
00466
00467 if((*itTrk)->isEmcShowerValid()) {
00468
00469 RecEmcShower* emcTrk = (*itTrk)->emcShower();
00470 double eraw = emcTrk->energy();
00471 double phi = emcTrk->phi();
00472 double the = emcTrk->theta();
00473
00474 HepLorentzVector ptrk;
00475 ptrk.setPx(eraw*sin(the)*cos(phi));
00476 ptrk.setPy(eraw*sin(the)*sin(phi));
00477 ptrk.setPz(eraw*cos(the));
00478 ptrk.setE(eraw);
00479 ptrk = ptrk.boost(-0.011,0,0);
00480
00481 eemc += ptrk.e();
00482 etot += ptrk.e();
00483
00484 }
00485
00486 }
00487
00488
00489
00490
00491 Vp4 pGam;
00492 pGam.clear();
00493 double eneu=0;
00494 for(int i = 0; i < nGam; i++) {
00495 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
00496 RecEmcShower* emcTrk = (*itTrk)->emcShower();
00497 double eraw = emcTrk->energy();
00498 double phi = emcTrk->phi();
00499 double the = emcTrk->theta();
00500 HepLorentzVector ptrk;
00501 ptrk.setPx(eraw*sin(the)*cos(phi));
00502 ptrk.setPy(eraw*sin(the)*sin(phi));
00503 ptrk.setPz(eraw*cos(the));
00504 ptrk.setE(eraw);
00505 ptrk = ptrk.boost(-0.011,0,0);
00506 pGam.push_back(ptrk);
00507 eneu += ptrk.e();
00508 etot += ptrk.e();
00509 eemc += ptrk.e();
00510 }
00511
00512 double esum = echarge + eneu;
00513
00514
00515 double maxE = 0.;
00516 double secE = 0.;
00517 double maxThe = 999.;
00518 double maxPhi = 999.;
00519 double secThe = 999.;
00520 double secPhi = 999.;
00521 int npart = 999.;
00522 double dphi = 999.;
00523 double dthe = 999.;
00524 int mdcHit1 = 0.;
00525 int mdcHit2 = 0.;
00526
00527 SmartDataPtr<RecEmcShowerCol> aShowerCol(eventSvc(),"/Event/Recon/RecEmcShowerCol");
00528 if (!aShowerCol) {
00529 log << MSG::WARNING << "Could not find RecEmcShowerCol" << endreq;
00530 return( StatusCode::SUCCESS);
00531 }
00532
00533 int ishower = 0;
00534 RecEmcShowerCol::iterator iShowerCol;
00535 for(iShowerCol=aShowerCol->begin();
00536 iShowerCol!=aShowerCol->end();
00537 iShowerCol++) {
00538
00539 if(ishower == 0) {
00540 maxE = (*iShowerCol)->e5x5();
00541 maxThe = (*iShowerCol)->theta();
00542 maxPhi = (*iShowerCol)->phi();
00543 npart = (*iShowerCol)->module();
00544 } else if(ishower == 1) {
00545 secE = (*iShowerCol)->e5x5();
00546 secThe = (*iShowerCol)->theta();
00547 secPhi = (*iShowerCol)->phi();
00548 } else if(ishower == 2) {
00549 break;
00550 }
00551
00552 ishower++;
00553
00554 }
00555
00556 if(aShowerCol->size() >= 2) {
00557
00558 dphi = (fabs(maxPhi-secPhi)-CLHEP::pi)*180./CLHEP::pi;
00559 dthe = (fabs(maxThe+secThe)-CLHEP::pi)*180./CLHEP::pi;
00560
00561 double phi1 = maxPhi<0 ? maxPhi+CLHEP::twopi : maxPhi;
00562 double phi2 = secPhi<0 ? secPhi+CLHEP::twopi : secPhi;
00563
00564
00565 double phi11=min(phi1,phi2);
00566 double phi22=max(phi1,phi2);
00567 double phi12=(phi11+phi22-CLHEP::pi)*0.5;
00568 double phi21=(phi11+phi22+CLHEP::pi)*0.5;
00569 if(phi12<0.) phi12 += CLHEP::twopi;
00570 if(phi21>CLHEP::twopi) phi21 -= CLHEP::twopi;
00571
00572 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc(), "/Event/Digi/MdcDigiCol");
00573 if (!mdcDigiCol) {
00574 log << MSG::FATAL << "Could not find MdcDigiCol!" << endreq;
00575 return StatusCode::FAILURE;
00576 }
00577 int hitnum = mdcDigiCol->size();
00578 for (int i = 0;i< hitnum;i++ ) {
00579 MdcDigi* digi= dynamic_cast<MdcDigi*>(mdcDigiCol->containedObject(i));
00580 double time = digi->getTimeChannel();
00581 double charge = digi->getChargeChannel();
00582 if (time == 0x7FFFFFFF || charge == 0x7FFFFFFF) continue;
00583 Identifier id(digi->identify());
00584 unsigned int iphi=MdcID::wire(id);
00585 unsigned int ilayer=MdcID::layer(id);
00586 if(ilayer>=43)
00587 log << MSG::ERROR << "MDC(" << ilayer <<","<<iphi <<")"<<endreq;
00588 double phi=CLHEP::twopi*iphi/idmax[ilayer];
00589 if(WhetherSector(phi,phi11,phi12)) mdcHit1++;
00590 else if(WhetherSector(phi,phi21,phi22)) mdcHit2++;
00591
00592
00593 }
00594 }
00595
00596
00597
00598
00599
00600
00601 if( eemc>m_bhabhaEmcECut && maxE>m_bhabhaMaxECut && secE>m_bhabhaSecECut
00602 && abs(dthe)<m_bhabhaDTheCut &&
00603 ( (dphi>m_bhabhaDPhiCut1&&dphi<m_bhabhaDPhiCut2) ||
00604 (dphi>m_bhabhaDPhiCut3&&dphi<m_bhabhaDPhiCut4) ) ) {
00605 if( npart==1 && mdcHit1>m_bhabhaMdcHitCutB && mdcHit2>m_bhabhaMdcHitCutB ) {
00606 m_isBarrelBhabha = true;
00607 m_barrelBhabhaNumber++;
00608 } else if( ( npart==0 || npart==2 )
00609 && mdcHit1>m_bhabhaMdcHitCutE && mdcHit2>m_bhabhaMdcHitCutE ) {
00610 m_isEndcapBhabha = true;
00611 m_endcapBhabhaNumber++;
00612 }
00613 }
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627 if(m_selectDimu) {
00628 if( m_dimuAlg->IsDimu() == 1 ) {
00629 m_isBarrelDimu = true;
00630 m_barrelDimuNumber++;
00631 } else if(m_dimuAlg->IsDimu() == 2 ) {
00632 m_isEndcapDimu = true;
00633 m_endcapDimuNumber++;
00634 }
00635 }
00636
00637
00638 if( (nGood>=1 && esum>m_hadronChaECut)
00639 || (nGood==0 && esum>m_hadronNeuECut) ) {
00640 m_isHadron = true;
00641 m_hadronNumber++;
00642 }
00643
00644
00645 if( nGood==0 && eemc>m_diphotonEmcECut && secE>m_diphotonSecECut
00646 && fabs(dthe)<m_diphotonDTheCut
00647 && dphi>m_diphotonDPhiCut1 && dphi<m_diphotonDPhiCut2 ) {
00648 if( npart==1 ) {
00649 m_isBarrelDiphoton = true;
00650 m_barrelDiphotonNumber++;
00651 } else if( npart==0 || npart==2 ) {
00652 m_isEndcapDiphoton = true;
00653 m_endcapDiphotonNumber++;
00654 }
00655 }
00656
00657
00658 if( m_selectBhabha && m_isBarrelBhabha ) m_subalg1->execute();
00659 if( m_selectBhabha && m_isEndcapBhabha ) m_subalg2->execute();
00660 if( m_selectDimu && m_isBarrelDimu ) m_subalg3->execute();
00661 if( m_selectDimu && m_isEndcapDimu ) m_subalg4->execute();
00662 if( m_selectHadron && m_isHadron ) m_subalg5->execute();
00663 if( m_selectDiphoton && m_isBarrelDiphoton ) m_subalg6->execute();
00664 if( m_selectDiphoton && m_isEndcapDiphoton ) m_subalg7->execute();
00665
00666
00667 if(m_output) {
00668 m_runnb = run;
00669 m_evtnb = event;
00670 m_esum = esum;
00671 m_eemc = eemc;
00672 m_etot = etot;
00673 m_nCharge = nCharge;
00674 m_nGood = nGood;
00675 m_nGam = nGam;
00676 m_ptot = ptot;
00677 m_pp = pp;
00678 m_pm = pm;
00679 m_maxE = maxE;
00680 m_secE = secE;
00681 m_dThe = dthe;
00682 m_dPhi = dphi;
00683 m_mdcHit1 = mdcHit1;
00684 m_mdcHit2 = mdcHit2;
00685 m_tuple0->write();
00686 }
00687
00688 return StatusCode::SUCCESS;
00689
00690 }
00691
00692
00693 StatusCode EventPreSelect::finalize() {
00694
00695 MsgStream log(msgSvc(), name());
00696 log << MSG::INFO << "in finalize()" << endmsg;
00697
00698 if(m_selectDimu&&m_output) {
00699 m_dimuAlg->Print();
00700 }
00701
00702 cout<<"Total events: "<<m_events<<endl;
00703
00704 if(m_selectBhabha) {
00705 cout << "Selected barrel bhabha: " << m_barrelBhabhaNumber << endl;
00706 cout << "Selected endcap bhabha: " << m_endcapBhabhaNumber << endl;
00707 }
00708
00709 if(m_selectDimu) {
00710 delete m_dimuAlg;
00711 cout << "Selected barrel dimu: " << m_barrelDimuNumber << endl;
00712 cout << "Selected endcap dimu: " << m_endcapDimuNumber << endl;
00713 }
00714
00715 if(m_selectHadron) {
00716 cout << "Selected hadron: " << m_hadronNumber << endl;
00717 }
00718
00719 if(m_selectDiphoton) {
00720 cout << "Selected barrel diphoton: " << m_barrelDiphotonNumber << endl;
00721 cout << "Selected endcap diphoton: " << m_endcapDiphotonNumber << endl;
00722 }
00723
00724 return StatusCode::SUCCESS;
00725 }
00726
00727 bool EventPreSelect::WhetherSector(double ph,double ph1,double ph2) {
00728 double phi1=min(ph1,ph2);
00729 double phi2=max(ph1,ph2);
00730 double delta=0.0610865;
00731 if((phi2-phi1)<CLHEP::pi){
00732 phi1 -=delta;
00733 phi2 +=delta;
00734 if(phi1<0.) phi1 += CLHEP::twopi;
00735 if(phi2>CLHEP::twopi) phi2 -= CLHEP::twopi;
00736 double tmp1=min(phi1,phi2);
00737 double tmp2=max(phi1,phi2);
00738 phi1=tmp1;
00739 phi2=tmp2;
00740 }
00741 else{
00742 phi1 +=delta;
00743 phi2 -=delta;
00744 }
00745 if((phi2-phi1)<CLHEP::pi){
00746 if(ph<=phi2&&ph>=phi1) return true;
00747 else return false;
00748 }
00749 else{
00750 if(ph>=phi2||ph<=phi1) return true;
00751 else return false;
00752 }
00753 }
00754
00755