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/EventHeader.h"
00010 #include "EventModel/Event.h"
00011 #include "TrigEvent/TrigEvent.h"
00012 #include "TrigEvent/TrigData.h"
00013 #include "McTruth/McParticle.h"
00014 #include "EvTimeEvent/RecEsTime.h"
00015
00016 #include "EvtRecEvent/EvtRecEvent.h"
00017 #include "EvtRecEvent/EvtRecTrack.h"
00018 #include "DstEvent/TofHitStatus.h"
00019
00020 #include "TMath.h"
00021 #include "GaudiKernel/INTupleSvc.h"
00022 #include "GaudiKernel/NTuple.h"
00023 #include "GaudiKernel/Bootstrap.h"
00024 #include "GaudiKernel/IHistogramSvc.h"
00025 #include "CLHEP/Vector/ThreeVector.h"
00026 using CLHEP::Hep3Vector;
00027 #include "CLHEP/Vector/LorentzVector.h"
00028 using CLHEP::HepLorentzVector;
00029 #include "CLHEP/Vector/TwoVector.h"
00030 using CLHEP::Hep2Vector;
00031 #include "CLHEP/Geometry/Point3D.h"
00032 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00033 typedef HepGeom::Point3D<double> HepPoint3D;
00034 #endif
00035 #include "VertexFit/KinematicFit.h"
00036 #include "VertexFit/VertexFit.h"
00037 #include "VertexFit/SecondVertexFit.h"
00038 #include "ParticleID/ParticleID.h"
00039 #include "HltEvent/DstHltInf.h"
00040 #include "HltEvent/HltInf.h"
00041
00042 #include "BbEmcAlg/BbEmc.h"
00043 #include <vector>
00044 typedef std::vector<int> Vint;
00045 typedef std::vector<HepLorentzVector> Vp4;
00046 const double velc = 299.792458;
00047 const double xmass[5] = {0.000511, 0.105658, 0.139570, 0.493677, 0.938272};
00048 const double pai=3.1415926;
00049
00050 static long m_cout_all(0), m_cout_tracks(0), m_cout_good(0), m_cout_dang(0);
00052 BbEmc::BbEmc(const std::string& name, ISvcLocator* pSvcLocator) :
00053 Algorithm(name, pSvcLocator) {
00054
00055
00056 declareProperty("hist", m_hist=false);
00057 declareProperty("Trigger", m_trigger_flag=true);
00058 declareProperty("Hlt", m_hlt_flag=true);
00059 declareProperty("Estime", m_est_flag=true);
00060
00061 declareProperty("KalTrk", m_kalTrk_flag=true);
00062
00063 declareProperty("EneCut", m_energy_cut=1.2);
00064 declareProperty("MomCut", m_mom_cut=0.04);
00065 declareProperty("DangCut", m_dangCut= 2.5);
00066
00067 declareProperty("Vr0cut", m_vr0cut=1.0);
00068 declareProperty("Vz0cut", m_vz0cut=5.0);
00069 }
00070
00071
00072 StatusCode BbEmc::initialize() {
00073 MsgStream log(msgSvc(), name());
00074 log << MSG::INFO << "in initialize()" << endmsg;
00075
00076 StatusCode status;
00077
00078 if(m_hist){
00079 NTuplePtr nt1(ntupleSvc(), "FILEBbEmc/bbEmc");
00080 if ( nt1 ) m_tuple1 = nt1;
00081 else {
00082 m_tuple1 = ntupleSvc()->book ("FILEBbEmc/bbEmc", CLID_ColumnWiseTuple, "BbEmc N-Tuple example");
00083 if ( m_tuple1 ){
00084 status = m_tuple1->addItem ("run", m_run);
00085 status = m_tuple1->addItem ("event", m_event);
00086 status = m_tuple1->addItem ("numCtrk", m_num_Ctrk);
00087 status = m_tuple1->addItem ("numNtrk", m_num_Ntrk);
00088 status = m_tuple1->addItem ("numGtrk", m_num_Gtrk);
00089
00090 status = m_tuple1->addItem ("trigindex", m_trig_index, 0, 48);
00091 status = m_tuple1->addIndexedItem("trigcond", m_trig_index, m_trig_cond);
00092 status = m_tuple1->addIndexedItem("trigchan", m_trig_index, m_trig_chan);
00093 status = m_tuple1->addItem ("timewindow", m_trig_timewindow);
00094 status = m_tuple1->addItem ("timetype", m_trig_timetype);
00095 status = m_tuple1->addItem ("hlttype", m_hlt_type );
00096
00097 status = m_tuple1->addItem ("estime", m_est_start);
00098 status = m_tuple1->addItem ("status", m_est_status);
00099 status = m_tuple1->addItem ("quality", m_est_quality);
00100
00101 status = m_tuple1->addItem ("dang", m_dang);
00102 status = m_tuple1->addItem ("index", m_index,0,2);
00103 status = m_tuple1->addIndexedItem ("theta", m_index, m_theta);
00104 status = m_tuple1->addIndexedItem ("phi", m_index, m_phi);
00105 status = m_tuple1->addIndexedItem ("ene", m_index, m_ene);
00106 } else {
00107 log << MSG::ERROR << "Cannot book N-tuple theone"<<endmsg;
00108 return StatusCode::FAILURE;
00109 }
00110 }
00111
00112 log << MSG::INFO << "end initialize()" << endmsg;
00113
00114 }
00115
00116 return StatusCode::SUCCESS;
00117
00118 }
00119
00120
00121 StatusCode BbEmc::execute(){
00122
00123 MsgStream log(msgSvc(),name());
00124 log<<MSG::INFO<<"in execute()"<<endreq;
00125 m_cout_all++;
00126
00127
00128 setFilterPassed(false);
00129
00130 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00131 int runNum= eventHeader->runNumber();
00132 int eventNum= eventHeader->eventNumber();
00133 if(m_cout_all %1000 ==0) {
00134 std::cout<<name()<<"::"<< m_cout_all <<" events executed"<<std::endl;
00135 }
00136
00137 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), "/Event/EvtRec/EvtRecEvent");
00138 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00139 << evtRecEvent->totalCharged() << " , "
00140 << evtRecEvent->totalNeutral() << " , "
00141 << evtRecEvent->totalTracks() <<endreq;
00142
00143 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(),"/Event/EvtRec/EvtRecTrackCol");
00144 if((evtRecEvent->totalTracks()<2) || (evtRecEvent->totalTracks()>30)){
00145 return StatusCode::SUCCESS;
00146 }
00147 m_cout_tracks ++;
00148
00149
00150 int hlt_type=0;
00151 if(m_hlt_flag){
00152 SmartDataPtr<HltInf> hlt(eventSvc(),"/Event/Hlt/HltInf");
00153 DstHltInf *aHlt;
00154 if(!hlt) {
00155 log << MSG::WARNING << "Could not find Event HltInf, try DstHltInf now" << endreq;
00156 SmartDataPtr<DstHltInf> hltDst(eventSvc(),"/Event/Hlt/DstHltInf");
00157 if(!hltDst){
00158 log << MSG::ERROR << "Could not find Event DstHltInf too, please re-generate data" << endreq;
00159 return StatusCode::FAILURE;
00160 } else {
00161 aHlt=hltDst;
00162 }
00163 } else {
00164 aHlt=hlt;
00165 }
00166 uint32_t temp_type(aHlt->getEventType());
00167 int mask(1);
00168 for(int i=0; i<32; i++){
00169 if(temp_type & mask){
00170 hlt_type = i;
00171 break;
00172 }
00173 temp_type >>= 1;
00174 }
00175 }
00176
00177
00178 SmartDataPtr<TrigData> trigData(eventSvc(),"/Event/Trig/TrigData");
00179 if(m_trigger_flag & runNum>0){
00180 if (!trigData) {
00181 log << MSG::FATAL << "Could not find TrigData in TDS" << endreq;
00182 return StatusCode::FAILURE;
00183 }
00184
00185 log << MSG::DEBUG << "Trigger conditions: " << endreq;
00186 for(int trig_id=0; trig_id < 48; trig_id++){
00187 log << MSG::DEBUG << "i:" << trig_id
00188 << " name:" << trigData->getTrigCondName(trig_id)
00189 << " cond:" << trigData->getTrigCondition(trig_id) << endreq;
00190 }
00191 }
00192
00193
00194 SmartDataPtr<RecEsTimeCol> recEstimeCol(eventSvc(), "/Event/Recon/RecEsTimeCol");
00195 if(m_est_flag & runNum>0){
00196 if(!recEstimeCol){
00197 log << MSG::WARNING << "Can not get RecEsTimeCol" << endreq;
00198 return StatusCode::SUCCESS;
00199 }
00200 log << MSG::DEBUG <<"size of EsTime: " << recEstimeCol->size() << endreq;
00201 }
00202
00203
00204 Vp4 vGood;
00205 HepLorentzVector m_lv_ele;
00206 for(int i=0; i<evtRecEvent->totalTracks(); i++){
00207 if(i>=evtRecTrkCol->size()) break;
00208
00209 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00210 if(!(*itTrk)->isEmcShowerValid()) continue;
00211
00212 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00213 if(emcTrk->energy()<m_energy_cut) continue;
00214
00215 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
00216 m_lv_ele.setVect(emcpos);
00217 m_lv_ele.setE(emcTrk->energy());
00218 vGood.push_back(m_lv_ele);
00219 }
00220
00221
00222 if(vGood.size() - 2) return SUCCESS;
00223 m_cout_good ++;
00224
00225
00226 double dang = vGood[0].vect().angle(vGood[1].vect());
00227 if(dang < m_dangCut) return SUCCESS;
00228 m_cout_dang ++;
00229
00230 if(vGood[0].e() > vGood[1].e()) swap(vGood[0],vGood[1]);
00231
00232
00233
00234 double d0,z0,cosTheta,phi,mom;
00235 Vint iGood; iGood.clear();
00236 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00237 double m_vz,m_vr;
00238 EvtRecTrackIterator itTrk;
00239 if (m_kalTrk_flag){
00240 itTrk=evtRecTrkCol->begin() + i;
00241 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00242 RecMdcKalTrack *mdcTrk = (*itTrk)->mdcKalTrack();
00243 m_vz = mdcTrk->z();
00244 m_vr = mdcTrk->r();
00245 }else{
00246 if(!(*itTrk)->isMdcTrackValid()) continue;
00247 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00248 m_vz = mdcTrk->z();
00249 m_vr = mdcTrk->r();
00250 }
00251 if(fabs(m_vz) >= m_vz0cut) continue;
00252 if(m_vr >= m_vr0cut) continue;
00253 iGood.push_back(i);
00254 }
00255
00256 if (m_hist){
00257 m_run = runNum;
00258 m_event = eventNum;
00259
00260 m_num_Ctrk = evtRecEvent->totalCharged();
00261 m_num_Ntrk = evtRecEvent->totalNeutral();
00262 m_num_Gtrk = iGood.size();
00263 if (trigData){
00264 for(int trig_id=0; trig_id < 48; trig_id++){
00265 m_trig_index = trig_id;
00266 if(m_trig_index<16){
00267 m_trig_chan[m_trig_index] = trigData->getTrigChannel(m_trig_index);
00268 }
00269 m_trig_cond[m_trig_index] = trigData->getTrigCondition(m_trig_index);
00270 }
00271 m_trig_timewindow = trigData->getTimeWindow();
00272 m_trig_timetype = trigData->getTimingType();
00273 }
00274
00275 m_hlt_type = hlt_type;
00276
00277 if (recEstimeCol){
00278 m_est_start = (*(recEstimeCol->begin()))->getTest();
00279 m_est_status = (*(recEstimeCol->begin()))->getStat();
00280 m_est_quality = (*(recEstimeCol->begin()))->getQuality();
00281 }
00282
00283 m_dang = dang;
00284 for(int i=0; i<2; i++){
00285 m_index = i;
00286 m_theta[m_index] = vGood[m_index].vect().theta();
00287 m_phi[m_index] = vGood[m_index].vect().phi();
00288 m_ene[m_index] = vGood[m_index].e();
00289 }
00290
00291 m_tuple1->write();
00292 }
00293
00294 setFilterPassed(true);
00295
00296 return StatusCode::SUCCESS;
00297 }
00298
00299 StatusCode BbEmc::finalize() {
00300 MsgStream log(msgSvc(),name());
00301 log<<MSG::INFO<<"in finalize()"<<endreq;
00302
00303 std::cout <<name()<<" total event: "<< m_cout_all << std::endl;
00304 std::cout <<name()<<" total tracks >= 2, <= 30: "<<m_cout_tracks << std::endl;
00305 std::cout <<name()<<" good showers = 2: "<< m_cout_good << std::endl;
00306 std::cout <<name()<<" angle between two showers: "<< m_cout_dang << std::endl;
00307
00308 return StatusCode::SUCCESS;
00309 }
00310