00001 #include "GaudiKernel/MsgStream.h"
00002 #include "GaudiKernel/AlgFactory.h"
00003 #include "GaudiKernel/ISvcLocator.h"
00004 #include "GaudiKernel/SmartDataPtr.h"
00005 #include "GaudiKernel/IDataProviderSvc.h"
00006 #include "GaudiKernel/PropertyMgr.h"
00007
00008 #include "EventModel/EventModel.h"
00009 #include "EventModel/Event.h"
00010 #include "EventModel/EventHeader.h"
00011 #include "EvtRecEvent/EvtRecEvent.h"
00012 #include "EvtRecEvent/EvtRecTrack.h"
00013
00014 #include "TMath.h"
00015 #include "GaudiKernel/INTupleSvc.h"
00016 #include "GaudiKernel/NTuple.h"
00017 #include "GaudiKernel/Bootstrap.h"
00018 #include "GaudiKernel/IHistogramSvc.h"
00019 #include "CLHEP/Vector/ThreeVector.h"
00020 #include "CLHEP/Vector/LorentzVector.h"
00021 #include "CLHEP/Vector/TwoVector.h"
00022 #include "MdcRawEvent/MdcDigi.h"
00023 using CLHEP::Hep3Vector;
00024 using CLHEP::Hep2Vector;
00025 using CLHEP::HepLorentzVector;
00026 #include "CLHEP/Geometry/Point3D.h"
00027 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00028 typedef HepGeom::Point3D<double> HepPoint3D;
00029 #endif
00030 #include "DigammaPreSelect/DigammaPreSelect.h"
00031 int selected =0;
00032
00033 #include <vector>
00034
00035 typedef std::vector<int> Vint;
00036 typedef std::vector<HepLorentzVector> Vp4;
00037 const double PI = 3.1415927;
00038 const double EBeam=1.548;
00040
00041 DigammaPreSelect::DigammaPreSelect(const std::string& name, ISvcLocator* pSvcLocator) :
00042 Algorithm(name, pSvcLocator) {
00043
00044
00045
00046 declareProperty ("Vr0cut", m_vr0cut=1.0);
00047 declareProperty ("Vz0cut", m_vz0cut=5.0);
00048
00049 declareProperty ("lowEnergyShowerCut", m_lowEnergyShowerCut=0.1);
00050 declareProperty ("highEnergyShowerCut", m_highEnergyShowerCut=0.5);
00051 declareProperty ("matchThetaCut", m_matchThetaCut = 0.2);
00052 declareProperty ("matchPhiCut", m_matchPhiCut = 0.2);
00053
00054 declareProperty ("highMomentumCut", m_highMomentumCut = 0.5);
00055 declareProperty ("EoPMaxCut", m_EoPMaxCut =1.3);
00056 declareProperty ("EoPMinCut", m_EoPMinCut = 0.6);
00057 declareProperty ("minAngShEnergyCut", m_minAngShEnergyCut = 0.2);
00058 declareProperty ("minAngCut", m_minAngCut = 0.3);
00059 declareProperty ("acolliCut", m_acolliCut = 0.03);
00060 declareProperty ("eNormCut", m_eNormCut = 0.5);
00061 declareProperty ("pNormCut", m_pNormCut = 0.5);
00062 declareProperty ("BarrelOrEndcap", m_BarrelOrEndcap = 1);
00063 declareProperty ("Output", m_output = false);
00064 declareProperty ("oneProngMomentumCut", m_oneProngMomentumCut =1.2);
00065
00066 }
00067
00068
00069 StatusCode DigammaPreSelect::initialize(){
00070 MsgStream log(msgSvc(), name());
00071
00072 log << MSG::INFO << "in initialize()" << endmsg;
00073
00074 if(m_output) {
00075 StatusCode status;
00076 NTuplePtr nt1(ntupleSvc(), "FILE1/bhabha");
00077 if ( nt1 ) m_tuple1 = nt1;
00078 else {
00079 m_tuple1 = ntupleSvc()->book ("FILE1/bhabha", CLID_ColumnWiseTuple, "N-Tuple example");
00080 if ( m_tuple1 ) {
00081 status = m_tuple1->addItem ("sh1_ene",m_sh1_ene);
00082 status = m_tuple1->addItem ("sh1_theta", m_sh1_theta);
00083 status = m_tuple1->addItem ("sh1_phi", m_sh1_phi);
00084 status = m_tuple1->addItem ("sh2_ene",m_sh2_ene);
00085 status = m_tuple1->addItem ("sh2_theta", m_sh2_theta);
00086 status = m_tuple1->addItem ("sh2_phi", m_sh2_phi);
00087 status = m_tuple1->addItem ("di_phi", m_di_phi);
00088 status = m_tuple1->addItem ("di_the", m_di_the);
00089 status = m_tuple1->addItem ("acolli", m_acolli);
00090 status = m_tuple1->addItem ("mdc_hit", m_mdc_hit);
00091 status = m_tuple1->addItem ("etot", m_etot);
00092
00093 }
00094 else {
00095 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00096 return StatusCode::FAILURE;
00097 }
00098 }
00099
00100
00101 NTuplePtr nt2(ntupleSvc(), "FILE1/bha1");
00102 if ( nt2 ) m_tuple2 = nt2;
00103 else {
00104 m_tuple2 = ntupleSvc()->book ("FILE1/bha1", CLID_ColumnWiseTuple, "N-Tuple example");
00105 if ( m_tuple2 ) {
00106 status = m_tuple2->addItem ("sh_ene",m_sh_ene);
00107 status = m_tuple2->addItem ("sh_theta", m_sh_theta);
00108 status = m_tuple2->addItem ("sh_phi", m_sh_phi);
00109 }
00110 else {
00111 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
00112 return StatusCode::FAILURE;
00113 }
00114 }
00115 }
00116
00117
00118
00119
00120
00121
00122 m_rejected=0;
00123 m_events=0;
00124 m_oneProngsSelected=0;
00125 m_twoProngsMatchedSelected=0;
00126 m_twoProngsOneMatchedSelected=0;
00127 m_selectedTrkID1=999;
00128 m_selectedTrkID2=999;
00129
00130 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00131 return StatusCode::SUCCESS;
00132
00133 }
00134
00135
00136 StatusCode DigammaPreSelect::execute() {
00137
00138 setFilterPassed(false);
00139
00140 MsgStream log(msgSvc(), name());
00141 log << MSG::INFO << "in execute()" << endreq;
00142
00143 m_events++;
00144
00145 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00146 if(!eventHeader)
00147 {
00148 cout<<" eventHeader "<<endl;
00149 return StatusCode::FAILURE;
00150 }
00151
00152 int run=eventHeader->runNumber();
00153 int event=eventHeader->eventNumber();
00154 if(event%1000==0) cout<<" run,event: "<<run<<","<<event<<endl;
00155
00156 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00157 if(!evtRecEvent ) {
00158 cout<<" evtRecEvent "<<endl;
00159 return StatusCode::FAILURE;
00160 }
00161
00162 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00163 << evtRecEvent->totalCharged() << " , "
00164 << evtRecEvent->totalNeutral() << " , "
00165 << evtRecEvent->totalTracks() <<endreq;
00166 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00167 if(!evtRecTrkCol){
00168 cout<<" evtRecTrkCol "<<endl;
00169 return StatusCode::FAILURE;
00170 }
00171
00172 Vint iGood;
00173 iGood.clear();
00174
00175 double ene5x5,theta,phi,eseed;
00176 double showerX,showerY,showerZ;
00177 long int thetaIndex,phiIndex;
00178
00179 RecEmcID showerId;
00180 unsigned int npart;
00181
00182 Vint ipip, ipim;
00183 ipip.clear();
00184 ipim.clear();
00185 Vp4 ppip, ppim;
00186 ppip.clear();
00187 ppim.clear();
00188
00189 vector<RecEmcShower* > GoodShowers;
00190 GoodShowers.clear();
00191 double etot=0;
00192 for(int i = 0; i< evtRecEvent->totalTracks(); i++) {
00193 if(i>=evtRecTrkCol->size()) break;
00194 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00195 if((*itTrk)->isEmcShowerValid()) {
00196 RecEmcShower *theShower = (*itTrk)->emcShower();
00197 etot+=theShower->e5x5();
00198 showerId = theShower->getShowerId();
00199 npart = EmcID::barrel_ec(showerId);
00200 ene5x5=theShower->e5x5();
00201 thetaIndex=EmcID::theta_module(showerId);
00202 phiIndex=EmcID::phi_module(showerId);
00203 if (ene5x5>0.4&&ene5x5<4&&npart==1&&(m_BarrelOrEndcap==1)){
00204 GoodShowers.push_back(theShower);
00205 iGood.push_back((*itTrk)->trackId());
00206 }
00207 else if (ene5x5>0.4&&ene5x5<4&&(npart==2||npart==0)&&(m_BarrelOrEndcap==2)){
00208 GoodShowers.push_back(theShower);
00209 iGood.push_back((*itTrk)->trackId());
00210 }
00211 else if (ene5x5>0.4&&ene5x5<4&&(m_BarrelOrEndcap==0)){
00212 GoodShowers.push_back(theShower);
00213 iGood.push_back((*itTrk)->trackId());
00214 }
00215 }
00216
00217 }
00218
00219
00220
00221 double MaxE(0), MaxPhi, MaxThe;
00222 int MaxId;
00223 double SecE(0), SecPhi, SecThe;
00224 for(int i=0; i<GoodShowers.size(); i++) {
00225 RecEmcShower *theShower = GoodShowers[i];
00226 double eraw = theShower->energy();
00227 if(eraw> MaxE) {
00228 MaxId =i;
00229 MaxE =eraw;
00230 MaxPhi = theShower->phi();
00231 MaxThe = theShower->theta();
00232 }
00233 }
00234 for(int i=0; i<GoodShowers.size(); i++) {
00235 RecEmcShower *theShower = GoodShowers[i];
00236 double eraw = theShower->energy();
00237 if(i!=MaxId&&eraw>SecE) {
00238 SecE =eraw;
00239 SecPhi = theShower->phi();
00240 SecThe = theShower->theta();
00241 }
00242 }
00243
00244 double dphi = (fabs(MaxPhi-SecPhi)-PI)*180./PI;
00245 double dthe = (fabs(MaxThe+SecThe)-PI)*180./PI;
00246 if(GoodShowers.size()>=2&&SecE>1.0&&dphi>-4&&dphi<2&&abs(dthe)<3&&etot>2.7) {
00247 setFilterPassed(true);
00248 selected++;
00249 }
00250
00251 if(m_output) {
00252 m_etot=etot;
00253 m_sh1_ene = MaxE;
00254 m_sh1_theta = MaxThe;
00255 m_sh1_phi = MaxPhi;
00256 m_sh2_ene = SecE;
00257 m_sh2_theta = SecThe;
00258 m_sh2_phi = SecPhi;
00259 m_di_phi = dphi;
00260 m_di_the = dthe;
00261 m_tuple1->write();
00262 }
00263
00264 return StatusCode::SUCCESS;
00265
00266 }
00267
00268
00269 StatusCode DigammaPreSelect::finalize() {
00270
00271 MsgStream log(msgSvc(), name());
00272 log << MSG::INFO << "in finalize()" << endmsg;
00273 cout<<"total events: "<<m_events<<endl;
00274 cout<<"selected digamma: "<<selected<<endl;
00275
00276 return StatusCode::SUCCESS;
00277 }
00278