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 #include "EventModel/EventModel.h"
00008 #include "EventModel/Event.h"
00009
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 using CLHEP::Hep3Vector;
00022 using CLHEP::Hep2Vector;
00023 using CLHEP::HepLorentzVector;
00024 #include "CLHEP/Geometry/Point3D.h"
00025 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00026 typedef HepGeom::Point3D<double> HepPoint3D;
00027 #endif
00028 #include "McCor/McCor.h"
00029
00030 #include "Identifier/Identifier.h"
00031
00032 #include "TGraphErrors.h"
00033 #include "TGraph2D.h"
00034 #include "TGraph2DErrors.h"
00035 #include <vector>
00036 #include <string>
00037 #include <fstream>
00038 #include <strstream>
00039
00040
00041
00042
00044
00045 TGraph2DErrors *dt = new TGraph2DErrors();
00046
00047 McCor::McCor(const std::string& name, ISvcLocator* pSvcLocator) :
00048 Algorithm(name, pSvcLocator) {
00049
00050
00051 ntOut = true;
00052 declareProperty("NTupleOut",ntOut);
00053
00054
00055
00056
00057
00058
00059
00060 }
00061
00062
00063 StatusCode McCor::initialize(){
00064 MsgStream log(msgSvc(), name());
00065 log << MSG::INFO << "in initialize()" << endmsg;
00066 StatusCode status;
00067 if(ntOut == true){
00068 NTuplePtr nt1(ntupleSvc(), "FILE1/ec");
00069 if ( nt1 ) m_tuple1 = nt1;
00070 else {
00071 m_tuple1 = ntupleSvc()->book ("FILE1/ec", CLID_ColumnWiseTuple, "ks N-Tuple example");
00072 if ( m_tuple1 ) {
00073 status = m_tuple1->addItem ("ef", m_ef);
00074 status = m_tuple1->addItem ("e5", m_e5);
00075 status = m_tuple1->addItem ("ec", m_ec);
00076 status = m_tuple1->addItem ("ct", m_ct);
00077 }
00078 else {
00079 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00080 return StatusCode::FAILURE;
00081 }
00082 }
00083 }
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094 double energy,thetaid,peak,peakerr,res,reserr;
00095 string DataPath;
00096 DataPath=getenv("MCCORROOT");
00097 DataPath += "/share/evset.txt";
00098 ifstream in1;
00099 in1.open(DataPath.c_str(),ios::in);
00100
00101 double ep[18]={0.03,0.04,0.05,0.075,0.1,0.125,0.15,0.2,0.25,0.3,0.4,0.5,0.75,1.0,1.25,1.5,1.75,2.0};
00102 for(int i=0;i<504;i++){
00103 in1>>energy;
00104 in1>>thetaid;
00105 in1>>peak;
00106 in1>>peakerr;
00107 in1>>res;
00108 in1>>reserr;
00109 int j = i/28;
00110 dt->SetPoint(i,energy,thetaid,peak);
00111 dt->SetPointError(i,0,0,peakerr);
00112 }
00113 in1.close();
00114 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00115 return StatusCode::SUCCESS;
00116
00117 }
00118
00119
00120 StatusCode McCor::execute() {
00121
00122
00123 MsgStream log(msgSvc(), name());
00124 log << MSG::INFO << "in execute()" << endreq;
00125
00126 SmartDataPtr<Event::EventH> eventHeader(eventSvc(),"/Event");
00127 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00128
00129 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00130 << evtRecEvent->totalCharged() << " , "
00131 << evtRecEvent->totalNeutral() << " , "
00132 << evtRecEvent->totalTracks() <<endreq;
00133
00134 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00135
00136
00137 for(int i = 0; i< evtRecEvent->totalTracks(); i++) {
00138
00139 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00140 if(!(*itTrk)->isEmcShowerValid()) continue;
00141 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00142
00143 int intid = emcTrk->cellId();
00144 Identifier id=Identifier(intid);
00145 int getthetaid = EmcID::theta_module(id);
00146 int getmodule = EmcID::barrel_ec(id);
00147 if(getthetaid>21)getthetaid=43-getthetaid;
00148 if(getmodule==1)getthetaid=getthetaid+6;
00149 double energyF = emcTrk->energy();
00150 double e5x5 = emcTrk->e5x5();
00151 double costheta = cos(emcTrk->theta());
00152 double dthetaid=double(getthetaid);
00153 double energyC = McCor::corEnergyMc(e5x5,dthetaid);
00154
00155 if(ntOut == true){
00156 m_ct = costheta;
00157 m_ef = energyF;
00158 m_e5 = e5x5;
00159 m_ec = energyC;
00160 m_tuple1->write();
00161 }
00162 emcTrk->setEnergy(energyC);
00163
00164 }
00165 }
00166
00167
00168 StatusCode McCor::finalize() {
00169
00170 MsgStream log(msgSvc(), name());
00171 log << MSG::INFO << "in finalize()" << endmsg;
00172 return StatusCode::SUCCESS;
00173 }
00174
00175
00176 double McCor::corEnergyMc(double eg,double theid){
00177
00178 double Energy5x5=eg;
00179 if(eg<0.029)eg=0.029;
00180 if(eg>1.76)eg=1.76;
00181 if(theid<=0)theid=0.001;
00182 if(theid>=27)theid=26.999;
00183 Float_t einter = eg + 0.00001;
00184 Float_t tinter = theid+0.0001;
00185 double ecor=dt->Interpolate(einter,tinter);
00186 if(!(ecor))return Energy5x5;
00187 if(ecor<0.5)return Energy5x5;
00188 double EnergyCor=Energy5x5/ecor;
00189 return EnergyCor;
00190 }
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203