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 "BhabhaPreSelect/BhabhaPreSelect.h"
00031 #include "BhabhaPreSelect/BhabhaType.h"
00032 int selected =0;
00033
00034 #include <vector>
00035
00036 typedef std::vector<int> Vint;
00037 typedef std::vector<HepLorentzVector> Vp4;
00038 const double PI = 3.1415927;
00039 const double EBeam=1.548;
00040 int BhabhaPreSelect::idmax[43]={40,44,48,56,64,72,80,80,76,76,
00041 88,88,100,100,112,112,128,128,140,140,
00042 160,160,160,160,176,176,176,176,208,208,
00043 208,208,240,240,240,240,256,256,256,256,
00044 288,288,288};
00045
00047
00048 BhabhaPreSelect::BhabhaPreSelect(const std::string& name, ISvcLocator* pSvcLocator) :
00049 Algorithm(name, pSvcLocator) {
00050
00051
00052
00053 declareProperty ("BarrelOrEndcap", m_BarrelOrEndcap = 1);
00054 declareProperty ("Output", m_output = false);
00055
00056
00057 }
00058
00059
00060 StatusCode BhabhaPreSelect::initialize(){
00061 MsgStream log(msgSvc(), name());
00062
00063 log << MSG::INFO << "in initialize()" << endmsg;
00064 cout<<" BarrelOrEndcap "<<m_BarrelOrEndcap<<endl;
00065
00066 if(m_output) {
00067 StatusCode status;
00068 NTuplePtr nt1(ntupleSvc(), "FILE1/bhabha");
00069 if ( nt1 ) m_tuple1 = nt1;
00070 else {
00071 m_tuple1 = ntupleSvc()->book ("FILE1/bhabha", CLID_ColumnWiseTuple, "N-Tuple example");
00072 if ( m_tuple1 ) {
00073 status = m_tuple1->addItem ("sh1_ene",m_sh1_ene);
00074 status = m_tuple1->addItem ("sh1_theta", m_sh1_theta);
00075 status = m_tuple1->addItem ("sh1_phi", m_sh1_phi);
00076 status = m_tuple1->addItem ("sh2_ene",m_sh2_ene);
00077 status = m_tuple1->addItem ("sh2_theta", m_sh2_theta);
00078 status = m_tuple1->addItem ("sh2_phi", m_sh2_phi);
00079 status = m_tuple1->addItem ("di_phi", m_di_phi);
00080 status = m_tuple1->addItem ("di_the", m_di_the);
00081 status = m_tuple1->addItem ("acolli", m_acolli);
00082 status = m_tuple1->addItem ("etot", m_etot);
00083 status = m_tuple1->addItem ("mdc_hit1", m_mdc_hit1);
00084 status = m_tuple1->addItem ("mdc_hit2", m_mdc_hit2);
00085
00086 }
00087 else {
00088 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00089 return StatusCode::FAILURE;
00090 }
00091 }
00092
00093
00094 NTuplePtr nt2(ntupleSvc(), "FILE1/bha1");
00095 if ( nt2 ) m_tuple2 = nt2;
00096 else {
00097 m_tuple2 = ntupleSvc()->book ("FILE1/bha1", CLID_ColumnWiseTuple, "N-Tuple example");
00098 if ( m_tuple2 ) {
00099 status = m_tuple2->addItem ("sh_ene",m_sh_ene);
00100 status = m_tuple2->addItem ("sh_theta", m_sh_theta);
00101 status = m_tuple2->addItem ("sh_phi", m_sh_phi);
00102 }
00103 else {
00104 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
00105 return StatusCode::FAILURE;
00106 }
00107 }
00108 }
00109
00110
00111
00112
00113
00114
00115 m_rejected=0;
00116 m_events=0;
00117 m_oneProngsSelected=0;
00118 m_twoProngsMatchedSelected=0;
00119 m_twoProngsOneMatchedSelected=0;
00120 m_selectedTrkID1=999;
00121 m_selectedTrkID2=999;
00122
00123 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00124 return StatusCode::SUCCESS;
00125
00126 }
00127
00128
00129 StatusCode BhabhaPreSelect::execute() {
00130
00131 MsgStream log(msgSvc(), name());
00132 log << MSG::INFO << "in execute()" << endreq;
00133
00134 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00135
00136
00137
00138
00139
00140
00141 int run=eventHeader->runNumber();
00142 int event=eventHeader->eventNumber();
00143
00144 if(event%1000==0) cout<<" run,event: "<<run<<","<<event<<endl;
00145
00146 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00147 if(!evtRecEvent ) {
00148 cout<<" evtRecEvent "<<endl;
00149 return StatusCode::FAILURE;
00150 }
00151
00152 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00153 << evtRecEvent->totalCharged() << " , "
00154 << evtRecEvent->totalNeutral() << " , "
00155 << evtRecEvent->totalTracks() <<endreq;
00156 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00157 if(!evtRecTrkCol){
00158 cout<<" evtRecTrkCol "<<endl;
00159 return StatusCode::FAILURE;
00160 }
00161
00162 m_events++;
00163 setFilterPassed(false);
00164
00165 Vint iGood;
00166 iGood.clear();
00167
00168 double ene5x5,theta,phi,eseed;
00169 double showerX,showerY,showerZ;
00170 long int thetaIndex,phiIndex;
00171
00172 RecEmcID showerId;
00173 unsigned int npart;
00174
00175 Vint ipip, ipim;
00176 ipip.clear();
00177 ipim.clear();
00178 Vp4 ppip, ppim;
00179 ppip.clear();
00180 ppim.clear();
00181
00182 vector<RecEmcShower* > GoodShowers;
00183 GoodShowers.clear();
00184 double etot=0;
00185 for(int i = 0; i< evtRecEvent->totalTracks(); i++) {
00186 if(i>=evtRecTrkCol->size()) break;
00187 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00188 if((*itTrk)->isEmcShowerValid()) {
00189 RecEmcShower *theShower = (*itTrk)->emcShower();
00190 etot+=theShower->e5x5();
00191 showerId = theShower->getShowerId();
00192 npart = EmcID::barrel_ec(showerId);
00193 ene5x5=theShower->e5x5();
00194 thetaIndex=EmcID::theta_module(showerId);
00195 phiIndex=EmcID::phi_module(showerId);
00196 if (ene5x5>0.4&&ene5x5<4&&npart==1&&(m_BarrelOrEndcap==1)){
00197 GoodShowers.push_back(theShower);
00198 iGood.push_back((*itTrk)->trackId());
00199 }
00200 else if (ene5x5>0.4&&ene5x5<4&&(npart==2||npart==0)&&(m_BarrelOrEndcap==2)){
00201 GoodShowers.push_back(theShower);
00202 iGood.push_back((*itTrk)->trackId());
00203 }
00204 else if (ene5x5>0.4&&ene5x5<4&&(m_BarrelOrEndcap==0)){
00205 GoodShowers.push_back(theShower);
00206 iGood.push_back((*itTrk)->trackId());
00207 }
00208 }
00209
00210 }
00211
00212
00213
00214 double MaxE(0), MaxPhi, MaxThe;
00215 int MaxId;
00216 double SecE(0), SecPhi, SecThe;
00217 for(int i=0; i<GoodShowers.size(); i++) {
00218 RecEmcShower *theShower = GoodShowers[i];
00219 double eraw = theShower->energy();
00220 if(eraw> MaxE) {
00221 MaxId =i;
00222 MaxE =eraw;
00223 MaxPhi = theShower->phi();
00224 MaxThe = theShower->theta();
00225 }
00226 }
00227 for(int i=0; i<GoodShowers.size(); i++) {
00228 RecEmcShower *theShower = GoodShowers[i];
00229 double eraw = theShower->energy();
00230 if(i!=MaxId&&eraw>SecE) {
00231 SecE =eraw;
00232 SecPhi = theShower->phi();
00233 SecThe = theShower->theta();
00234 }
00235 }
00236
00237 double dphi = (fabs(MaxPhi-SecPhi)-PI)*180./PI;
00238 double dthe = (fabs(MaxThe+SecThe)-PI)*180./PI;
00239
00240 int total1 = 0;
00241 int total2 = 0;
00242 if(GoodShowers.size()>=2 && MaxE>1.0 && abs(dthe)<3
00243 && ((dphi>-25&&dphi<-4)||(dphi>2&&dphi<20)) && etot>2.7) {
00244
00245 double phi1 = MaxPhi<0 ? MaxPhi+CLHEP::twopi : MaxPhi;
00246 double phi2 = SecPhi<0 ? SecPhi+CLHEP::twopi : SecPhi;
00247
00248
00249 double phi11=min(phi1,phi2);
00250 double phi22=max(phi1,phi2);
00251 double phi12=(phi11+phi22-CLHEP::pi)*0.5;
00252 double phi21=(phi11+phi22+CLHEP::pi)*0.5;
00253 if(phi12<0.) phi12 += CLHEP::twopi;
00254 if(phi21>CLHEP::twopi) phi21 -= CLHEP::twopi;
00255
00256 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc(), "/Event/Digi/MdcDigiCol");
00257 if (!mdcDigiCol) {
00258 log << MSG::FATAL << "Could not find MdcDigiCol!" << endreq;
00259 return StatusCode::FAILURE;
00260 }
00261 int hitnum = mdcDigiCol->size();
00262 for (int i = 0;i< hitnum;i++ ) {
00263 MdcDigi* digi= dynamic_cast<MdcDigi*>(mdcDigiCol->containedObject(i));
00264 double time = digi->getTimeChannel();
00265 double charge = digi->getChargeChannel();
00266 if (time == 0x7FFFFFFF || charge == 0x7FFFFFFF) continue;
00267 Identifier id(digi->identify());
00268 unsigned int iphi=MdcID::wire(id);
00269 unsigned int ilayer=MdcID::layer(id);
00270 if(ilayer>=43)
00271 log << MSG::ERROR << "MDC(" << ilayer <<","<<iphi <<")"<<endreq;
00272 double phi=CLHEP::twopi*iphi/idmax[ilayer];
00273 if(WhetherSector(phi,phi11,phi12)) total1++;
00274 else if(WhetherSector(phi,phi21,phi22)) total2++;
00275
00276
00277 }
00278
00279 if ( (m_BarrelOrEndcap==1&&total1>15 &&total2>15)
00280 || (m_BarrelOrEndcap!=1&&total1>5 &&total2>5) ) {
00281 setFilterPassed(true);
00282 selected++;
00283 }
00284 }
00285
00286 if(m_output) {
00287 m_etot=etot;
00288 m_mdc_hit1=total1;
00289 m_mdc_hit2=total2;
00290 m_sh1_ene = MaxE;
00291 m_sh1_theta = MaxThe;
00292 m_sh1_phi = MaxPhi;
00293 m_sh2_ene = SecE;
00294 m_sh2_theta = SecThe;
00295 m_sh2_phi = SecPhi;
00296 m_di_phi = dphi;
00297 m_di_the = dthe;
00298 m_tuple1->write();
00299 }
00300
00301 return StatusCode::SUCCESS;
00302
00303 }
00304
00305
00306 StatusCode BhabhaPreSelect::finalize() {
00307
00308 MsgStream log(msgSvc(), name());
00309 log << MSG::INFO << "in finalize()" << endmsg;
00310 cout<<"m_events "<<m_events<<endl;
00311 cout<<" selected "<<selected<<endl;
00312 cout<<"m_rejected "<<m_rejected<<endl;
00313 cout<<"m_oneProngsSelected "<<m_oneProngsSelected<<endl;
00314 cout<<"m_twoProngsMatchedSelected "<<m_twoProngsMatchedSelected<<endl;
00315 cout<<"m_twoProngsOneMatchedSelected "<<m_twoProngsOneMatchedSelected<<endl;
00316
00317 return StatusCode::SUCCESS;
00318 }
00319
00320 bool BhabhaPreSelect::WhetherSector(double ph,double ph1,double ph2){
00321 double phi1=min(ph1,ph2);
00322 double phi2=max(ph1,ph2);
00323 double delta=0.0610865;
00324 if((phi2-phi1)<CLHEP::pi){
00325 phi1 -=delta;
00326 phi2 +=delta;
00327 if(phi1<0.) phi1 += CLHEP::twopi;
00328 if(phi2>CLHEP::twopi) phi2 -= CLHEP::twopi;
00329 double tmp1=min(phi1,phi2);
00330 double tmp2=max(phi1,phi2);
00331 phi1=tmp1;
00332 phi2=tmp2;
00333 }
00334 else{
00335 phi1 +=delta;
00336 phi2 -=delta;
00337 }
00338 if((phi2-phi1)<CLHEP::pi){
00339 if(ph<=phi2&&ph>=phi1) return true;
00340 else return false;
00341 }
00342 else{
00343 if(ph>=phi2||ph<=phi1) return true;
00344 else return false;
00345 }
00346 }
00347