00001 #include "BbEmcAlg/MdcBbEmcEff.h"
00002 #include "GaudiKernel/MsgStream.h"
00003 #include "GaudiKernel/SmartDataPtr.h"
00004 #include "EventModel/EventHeader.h"
00005 #include "RawEvent/RawDataUtil.h"
00006 #include "Identifier/MdcID.h"
00007 #include "MdcRawEvent/MdcDigi.h"
00008 #include "EvTimeEvent/RecEsTime.h"
00009 #include "EvtRecEvent/EvtRecEvent.h"
00010 #include "EvtRecEvent/EvtRecTrack.h"
00011 #include "CLHEP/Units/PhysicalConstants.h"
00012
00013
00014 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00015
00016 typedef HepGeom::Point3D<double> HepPoint3D;
00017 #endif
00018
00019
00020 #include "MdcRecEvent/RecMdcTrack.h"
00021 #include "MdcRecEvent/RecMdcHit.h"
00022
00023
00024 #include "GaudiKernel/INTupleSvc.h"
00025 #include "GaudiKernel/NTuple.h"
00026
00027 typedef std::vector<HepLorentzVector> Vp4;
00028 using namespace std;
00029 using namespace Event;
00030
00031 MdcBbEmcEff::MdcBbEmcEff(const std::string& name, ISvcLocator* pSvcLocator) :
00032 Algorithm(name, pSvcLocator) {
00033 declareProperty("hist", m_hist = false);
00034 declareProperty("debug", m_debug= 0);
00035
00036 declareProperty("emcEneCutLow", m_emcEneCutLow=1.44);
00037 declareProperty("emcEneCutHigh", m_emcEneCutHigh=1.88);
00038 declareProperty("emcEneCutTot", m_emcEneCutTot=3.16);
00039 declareProperty("emcDangCutLow", m_emcDangCutLow=2.94);
00040 declareProperty("emcDangCutHigh",m_emcDangCutHigh=3.08);
00041
00042 declareProperty("dPhi", m_dPhiCut= 0.2);
00043 declareProperty("dCosTheta", m_dCosThetaCut= 0.2);
00044 declareProperty("d0", m_d0Cut= 1.);
00045 declareProperty("z0", m_z0Cut= 5.);
00046
00047 declareProperty("barrelCut", m_barrelCut = 0.8);
00048 declareProperty("endcapCut", m_endcapCutLow = 0.83);
00049 declareProperty("endcapCutHigh", m_endcapCutHigh = 0.93);
00050 }
00051
00052
00053
00054 StatusCode MdcBbEmcEff::initialize(){
00055 MsgStream log(msgSvc(), name());
00056 StatusCode sc;
00057 m_evtIndex = 0;
00058
00059 m_effAllN1 = 0;
00060 m_effAllN2 = 0;
00061
00062 for (int i=0;i<3;i++){
00063 m_effN1[i] = 0;
00064 m_effN2[i] = 0;
00065 }
00066
00067 if(m_hist) {
00068 if(bookNTuple()<0) sc = StatusCode::FAILURE;
00069 }
00070
00071 return StatusCode::SUCCESS;
00072 }
00073
00074
00075 StatusCode MdcBbEmcEff::execute() {
00076 MsgStream log(msgSvc(), name());
00077 StatusCode sc = StatusCode::SUCCESS;
00078
00079 setFilterPassed(false);
00080
00081
00082 if(getEventInfo()<0) return StatusCode::FAILURE;
00083
00084
00085 if(selectBbByEmcShower()<0) return StatusCode::SUCCESS;
00086
00087
00088 if(bbEmcMdcTrackingEff()<0) return StatusCode::SUCCESS;
00089
00090 return StatusCode::SUCCESS;
00091 }
00092
00093
00094 StatusCode MdcBbEmcEff::finalize() {
00095 MsgStream log(msgSvc(), name());
00096 log << MSG::INFO << "in finalize()" << endreq;
00097 if((m_effAllN1+m_effAllN2)>0)std::cout<<" MdcBbEmcEff efficiency = N2/(N1+N2) = "
00098 << m_effAllN2<<"/("<<m_effAllN1<<"+"<<m_effAllN2<<") = "
00099 << (m_effAllN2)/((float)(m_effAllN1+m_effAllN2))<<std::endl;
00100 for(int i=0;i<3;i++){
00101 string pos;
00102 if (0==i) pos="barrel";
00103 if (1==i) pos="gap";
00104 if (2==i) pos="endcap";
00105 if((m_effN1[i]+m_effN2[i])>0){
00106 std::cout<<" MdcBbEmcEff of "<<pos <<" N2/(N1+N2) = "
00107 << m_effN2[i]<<"/("<<m_effN1[i]<<"+"<<m_effN2[i]<<") = "
00108 << (m_effN2[i])/((float)(m_effN1[i]+m_effN2[i]))<<std::endl;
00109 }
00110 }
00111 return StatusCode::SUCCESS;
00112 }
00113
00114
00115 int MdcBbEmcEff::bookNTuple(){
00116 MsgStream log(msgSvc(), name());
00117 StatusCode sc;
00118 NTuplePtr nt1(ntupleSvc(), "MdcBbEmcEff/evt");
00119 if ( nt1 ) { m_tuple1 = nt1;}
00120 else {
00121 m_tuple1 = ntupleSvc()->book ("MdcBbEmcEff/evt", CLID_ColumnWiseTuple, "event");
00122 if ( m_tuple1 ) {
00123
00124 sc = m_tuple1->addItem ("runNo", m_runNo);
00125 sc = m_tuple1->addItem ("evtNo", m_evtNo);
00126 sc = m_tuple1->addItem ("t0", m_t0);
00127 sc = m_tuple1->addItem ("t0Stat", m_t0Stat);
00128
00129
00130 sc = m_tuple1->addItem ("index", m_index,0,2);
00131 sc = m_tuple1->addIndexedItem ("emcTheta",m_index, m_emcTheta);
00132 sc = m_tuple1->addIndexedItem ("emcPhi", m_index, m_emcPhi);
00133 sc = m_tuple1->addIndexedItem ("emcEne", m_index, m_emcEne);
00134 sc = m_tuple1->addItem ("emcDang", m_emcDang);
00135
00136
00137 sc = m_tuple1->addItem ("dCosTheta", m_dCosTheta);
00138 sc = m_tuple1->addItem ("dPhi", m_dPhi);
00139 sc = m_tuple1->addItem ("nTk", m_nTk,0,2);
00140 sc = m_tuple1->addIndexedItem ("phi", m_nTk,m_phi);
00141 sc = m_tuple1->addIndexedItem ("cosTheta",m_nTk,m_cosTheta);
00142 sc = m_tuple1->addIndexedItem ("d0", m_nTk,m_d0);
00143 sc = m_tuple1->addIndexedItem ("z0", m_nTk,m_z0);
00144 sc = m_tuple1->addIndexedItem ("p", m_nTk,m_p);
00145 sc = m_tuple1->addIndexedItem ("pt", m_nTk,m_pt);
00146 } else {
00147 log << MSG::ERROR << "Cannot book tuple MdcBbEmcEff/evt" << endmsg;
00148 return -1;
00149 }
00150 }
00151 return 1;
00152 }
00153
00154 int MdcBbEmcEff::getEventInfo(){
00155 MsgStream log(msgSvc(), name());
00156
00157
00158 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
00159 if (evtHead) {
00160 t_runNo = evtHead->runNumber();
00161 t_evtNo = evtHead->eventNumber();
00162 }else{
00163 log << MSG::WARNING<< "Could not retreve event header" << endreq;
00164 }
00165
00166 std::cout<<m_evtIndex <<"------------evtNo:"<<t_evtNo << std::endl;
00167 m_evtIndex++;
00168
00169
00170 t_t0 = -1;
00171 t_t0Stat = -1;
00172 SmartDataPtr<RecEsTimeCol> aevtmeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00173 if ((!aevtmeCol)||(aevtmeCol->size()==0)) {
00174 log << MSG::WARNING << "Could not fnd RecEsTimeCol" << endreq;
00175 }else{
00176 RecEsTimeCol::iterator iter_evt = aevtmeCol->begin();
00177 t_t0 = (*iter_evt)->getTest();
00178 t_t0Stat = (*iter_evt)->getStat();
00179 }
00180 return 1;
00181 }
00182
00183 int MdcBbEmcEff::selectBbByEmcShower(){
00184 MsgStream log(msgSvc(), name());
00185
00186 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), "/Event/EvtRec/EvtRecEvent");
00187 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00188 << evtRecEvent->totalCharged() << " , "
00189 << evtRecEvent->totalNeutral() << " , "
00190 << evtRecEvent->totalTracks() <<endreq;
00191 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), "/Event/EvtRec/EvtRecTrackCol");
00192
00193 Vp4 vGood;
00194 HepLorentzVector m_lv_ele;
00195 for(int i=0; i<evtRecEvent->totalTracks(); i++){
00196 if(i>=evtRecTrkCol->size()) return -1;
00197
00198 EvtRecTrackIterator itTrk =evtRecTrkCol->begin() + i;
00199 if(!(*itTrk)->isEmcShowerValid()) {
00200 if(m_debug>1)std::cout<<"EmcShower NOT Valid "<<std::endl;
00201 continue;
00202 }
00203 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00204 if((emcTrk->energy()>m_emcEneCutHigh)||(emcTrk->energy()<m_emcEneCutLow)){
00205 if(m_debug>1)std::cout<<"Cut by EmcEnergy: "<<emcTrk->energy()<<std::endl;
00206 continue;
00207 }
00208 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
00209 m_lv_ele.setVect(emcpos);
00210 m_lv_ele.setE(emcTrk->energy());
00211
00212 vGood.push_back(m_lv_ele);
00213 }
00214
00215 if(m_debug>1)std::cout<<"vGood.size = "<<vGood.size()<<std::endl;
00216 if(vGood.size()!=2){
00217 if (m_debug>0) std::cout<<"Cut by vGood.size: "<<vGood.size()<<std::endl;
00218 return -2;
00219 }
00220
00221
00222
00223
00224 for(int i=0; i<2; i++){
00225 double cosEmcTheta = cos(vGood[i].vect().theta());
00226 if( fabs(cosEmcTheta) <= m_barrelCut ){
00227 m_posFlag = BARREL;
00228 }else if((cosEmcTheta>=m_endcapCutLow)&&(cosEmcTheta<=m_endcapCutHigh)){
00229 m_posFlag = ENDCAP;
00230 break;
00231 }else if((cosEmcTheta>m_barrelCut)&&(cosEmcTheta<m_endcapCutLow)){
00232 m_posFlag = GAP;
00233 }else{
00234 m_posFlag = OUT;
00235 }
00236 if(m_debug>1) std::cout<<"Emc cos(theta)="<<cosEmcTheta
00237 <<"Track in "<<m_posFlag<<std::endl;
00238 }
00239 if(m_posFlag == OUT) return -5;
00240
00241 double dang = vGood[0].vect().angle(vGood[1].vect());
00242 if(vGood[0].e() > vGood[1].e()) swap(vGood[0],vGood[1]);
00243 double emcTheta[2],emcEne[2];
00244 for(int index=0; index<2; index++){
00245 emcTheta[index] = vGood[index].vect().theta();
00246 t_emcPhi[index] = vGood[index].vect().phi();
00247 emcEne[index] = vGood[index].e();
00248 }
00249
00250 if(m_hist){
00251 m_index = 2;
00252 for (int index=0;index<2;index++){
00253 m_emcTheta[index] = emcTheta[index];
00254 m_emcPhi[index] = t_emcPhi[index];
00255 m_emcEne[index] = emcEne[index];
00256 }
00257 m_emcDang = dang;
00258 }
00259 if(m_debug>1) std::cout<<" dang "<<dang<<std::endl;
00260 if(m_debug>1) std::cout<<" energy "<<emcEne[0]<<" "<<emcEne[1]<<std::endl;
00261
00262 if(dang<=m_emcDangCutLow || dang>=m_emcDangCutHigh ) {
00263 if(m_debug>0) std::cout<<"Cut by emcDang "<<dang<<std::endl;
00264 return -3;
00265 }
00266
00267 if(emcEne[0]<=m_emcEneCutLow || emcEne[1]>=m_emcEneCutHigh){
00268 if(m_debug>0) std::cout<<"Cut by emc energy low or high "
00269 <<emcEne[0]<<" "<<emcEne[1]<<std::endl;
00270 return -4;
00271 }
00272 if( (emcEne[0] + emcEne[1])<=m_emcEneCutTot){
00273 if(m_debug>0) std::cout<<"Cut by emc total energy:"
00274 <<emcEne[0]<<" "<<emcEne[1]<<" tot:"<<emcEne[0]+emcEne[1]<<std::endl;
00275 return StatusCode::FAILURE;
00276 }
00277
00278 return 1;
00279 }
00280
00281 int MdcBbEmcEff::bbEmcMdcTrackingEff(){
00282 MsgStream log(msgSvc(), name());
00283
00284
00285 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
00286 if (!recMdcTrackCol){
00287 log << MSG::WARNING<< " Unable to retrieve recMdcTrackCol" << endreq;
00288 return -1;
00289 }
00290
00291
00292 t_nTk = recMdcTrackCol->size();
00293 if((t_nTk>2) || (0==t_nTk)) {
00294 if(m_debug>1)std::cout<<name()<<"Cut by nTk "<<t_nTk<<std::endl;
00295 return -2;
00296 }
00297
00298
00299 double d0[2],z0[2],cosTheta[2],phi[2],p[2],pt[2];
00300 RecMdcTrackCol::iterator it = recMdcTrackCol->begin();
00301 for(int iTk=0 ; it != recMdcTrackCol->end(); it++,iTk++ ) {
00302 RecMdcTrack* tk = (*it);
00303 d0[iTk] = tk->helix(0);
00304 z0[iTk] = tk->helix(3);
00305 if ((fabs(d0[iTk])>m_d0Cut) || (fabs(z0[iTk])>m_z0Cut)){
00306 if(m_debug>0) std::cout<<"Cut by d0 "
00307 <<d0[iTk]<<" z0 "<<z0[iTk]<<std::endl;
00308 return -3;
00309 }
00310 cosTheta[iTk] = cos(tk->theta());
00311 phi[iTk] = tk->phi();
00312 p[iTk] = tk->p();
00313 pt[iTk] = tk->pxy();
00314 }
00315
00316 double dCosTheta = cosTheta[0]+cosTheta[1];
00317 double dPhi = phi[0] - phi[1] + CLHEP::pi;
00318 if(dPhi > CLHEP::pi) dPhi-=CLHEP::twopi;
00319
00320
00321 if (m_hist){
00322 m_runNo = t_runNo;
00323 m_evtNo = t_evtNo;
00324 m_t0 = t_t0;
00325 m_t0Stat= t_t0Stat;
00326
00327 for(int i=0;i<2;i++){
00328 m_cosTheta[i]=cosTheta[i];
00329 m_phi[i]=phi[i];
00330 m_d0[i]=d0[i];
00331 m_z0[i]=z0[i];
00332 m_p[i]=p[i];
00333 m_pt[i]=pt[i];
00334 }
00335 m_nTk = t_nTk;
00336 m_dCosTheta= dCosTheta;
00337 m_dPhi = dPhi;
00338
00339 m_nTk = t_nTk;
00340 m_tuple1->write();
00341 }
00342
00343
00344 if(2==t_nTk){
00345 if((fabs(dCosTheta)>m_dCosThetaCut) || (fabs(dPhi)>m_dPhiCut)){
00346 if(m_debug>0) {
00347 if (fabs(dCosTheta)>m_dCosThetaCut){
00348 std::cout<<"Cut by dCosTheta "<<dCosTheta<<std::endl;
00349 }else{
00350 std::cout<<"Cut by dPhi "<<dPhi<<std::endl;
00351 }
00352 }
00353 return -4;
00354 }
00355 m_effAllN2++;
00356 m_effN2[m_posFlag]++;
00357 }else{
00358 m_effAllN1++;
00359 m_effN1[m_posFlag]++;
00360 }
00361
00362
00363 if ((1==t_nTk)&&(m_posFlag==BARREL)) setFilterPassed(true);
00364 return 1;
00365 }