3#include "GaudiKernel/MsgStream.h"
4#include "GaudiKernel/SmartDataPtr.h"
5#include "GaudiKernel/PropertyMgr.h"
6#include "GaudiKernel/Bootstrap.h"
8#include "GaudiKernel/INTupleSvc.h"
9#include "GaudiKernel/NTuple.h"
10#include "GaudiKernel/ITHistSvc.h"
12#include "CLHEP/Vector/ThreeVector.h"
13#include "CLHEP/Vector/LorentzVector.h"
15#include "EventModel/EventModel.h"
16#include "EventModel/Event.h"
18#include "EvtRecEvent/EvtRecEvent.h"
19#include "EvtRecEvent/EvtRecTrack.h"
20#include "DstEvent/TofHitStatus.h"
21#include "EventModel/EventHeader.h"
24#include "ParticleID/ParticleID.h"
26#include "DQAEvent/DQAEvent.h"
27#include "DQA_EMC/DQA_EMC.h"
29#ifndef ENABLE_BACKWARDS_COMPATIBILITY
32using CLHEP::HepLorentzVector;
38 Algorithm(name, pSvcLocator) {
41 declareProperty(
"NtupleOutput", m_NtupleOutput=0);
47 MsgStream log(
msgSvc(), name());
49 log << MSG::INFO <<
"in initialize()" << endmsg;
57 if ( nt ) m_tuple = nt;
59 m_tuple =
ntupleSvc()->book(
"DQAFILE/EMC", CLID_ColumnWiseTuple,
"EMC ntuple");
61 status = m_tuple->addItem(
"ixtal", m_ixtal);
62 status = m_tuple->addItem(
"npart", m_npart);
63 status = m_tuple->addItem(
"ntheta",m_ntheta);
64 status = m_tuple->addItem(
"nphi", m_nphi);
65 status = m_tuple->addItem(
"theta", m_theta);
66 status = m_tuple->addItem(
"phi", m_phi);
67 status = m_tuple->addItem(
"emcX", m_emcX);
68 status = m_tuple->addItem(
"emcY", m_emcY);
69 status = m_tuple->addItem(
"eSeed", m_eSeed);
70 status = m_tuple->addItem(
"e5x5", m_e5x5);
71 status = m_tuple->addItem(
"energy",m_energy);
72 status = m_tuple->addItem(
"time", m_time);
75 log << MSG::ERROR <<
"Can not book N-tuple:" << long(m_tuple) << endreq;
79 if(service(
"THistSvc", m_thistsvc).isFailure()) {
80 log << MSG::ERROR <<
"Couldn't get THistSvc" << endreq;
81 return StatusCode::FAILURE;
87 std::string HistName0, HistName;
88 HistName0=
"/DQAHist/EMC/";
90 sprintf( name,
"EMC_Bhabha_ShowerEneBarrelVsEvent");
93 m_HistEnergyB =
new TH1F(name,
title,200,0.,2.0);
94 m_HistEnergyB->GetXaxis()->SetTitle(
"shower energy(GeV)");
95 m_HistEnergyB->GetXaxis()->CenterTitle();
96 m_HistEnergyB->GetYaxis()->SetTitle(
"Number of event");
97 m_HistEnergyB->GetYaxis()->CenterTitle();
99 HistName=HistName0+name;
100 if(m_thistsvc->regHist(HistName, m_HistEnergyB).isFailure()){
101 log << MSG::ERROR <<
"Couldn't register " <<name<< endreq;
104 sprintf( name,
"EMC_Bhabha_ShowerEneEastVsEvent");
106 m_HistEnergyEast =
new TH1F(name,
title,200,0.,2.0);
107 m_HistEnergyEast->GetXaxis()->SetTitle(
"shower energy(GeV)");
108 m_HistEnergyEast->GetXaxis()->CenterTitle();
109 m_HistEnergyEast->GetYaxis()->SetTitle(
"Number of event");
110 m_HistEnergyEast->GetYaxis()->CenterTitle();
112 HistName=HistName0+name;
113 if(m_thistsvc->regHist(HistName, m_HistEnergyEast).isFailure()){
114 log << MSG::ERROR <<
"Couldn't register "<<name << endreq;
118 sprintf( name,
"EMC_Bhabha_ShowerEneWestVsEvent");
120 m_HistEnergyWest =
new TH1F(name,
title,200,0.,2.0);
121 m_HistEnergyWest->GetXaxis()->SetTitle(
"shower energy(GeV)");
122 m_HistEnergyWest->GetXaxis()->CenterTitle();
123 m_HistEnergyWest->GetYaxis()->SetTitle(
"Number of event");
124 m_HistEnergyWest->GetYaxis()->CenterTitle();
127 HistName=HistName0+name;
128 if( m_thistsvc->regHist(HistName, m_HistEnergyWest).isFailure()){
129 log << MSG::ERROR <<
"Couldn't register "<<name << endreq;
133 sprintf( name,
"EMC_Bhabha_ShowerThetaVsvent");
135 m_HistTheta =
new TH1F(name,
title,56, 0, 56);
136 m_HistTheta->GetXaxis()->SetTitle(
"shower ID(theta)");
137 m_HistTheta->GetXaxis()->CenterTitle();
138 m_HistTheta->GetYaxis()->SetTitle(
"Number of event");
139 m_HistTheta->GetYaxis()->CenterTitle();
141 HistName=HistName0+name;
142 if( m_thistsvc->regHist(HistName, m_HistTheta).isFailure()){
143 log << MSG::ERROR <<
"Couldn't register" <<name<< endreq;
146 sprintf(name,
"EMC_Bhabha_ShowerCosTheta");
148 m_HistCosTheta =
new TH1F(name,
title,200, -1.0, 1.0);
149 m_HistCosTheta->GetXaxis()->SetTitle(
"shower cos(theta)");
150 m_HistCosTheta->GetXaxis()->CenterTitle();
151 m_HistCosTheta->GetYaxis()->SetTitle(
"Number of event");
152 m_HistCosTheta->GetYaxis()->CenterTitle();
154 HistName=HistName0+name;
155 if(m_thistsvc->regHist(HistName, m_HistCosTheta).isFailure()){
156 log << MSG::ERROR <<
"Couldn't register "<<name << endreq;
160 sprintf( name,
"EMC_Bhabha_ShowerPhiBarrelVsEvent");
162 m_HistPhiB =
new TH1F(name,
title,120, 0, 120);
163 m_HistPhiB->GetXaxis()->SetTitle(
"shower ID(phi)");
164 m_HistPhiB->GetXaxis()->CenterTitle();
165 m_HistPhiB->GetYaxis()->SetTitle(
"Number of event");
166 m_HistPhiB->GetYaxis()->CenterTitle();
169 HistName=HistName0+name;
170 if( m_thistsvc->regHist(HistName, m_HistPhiB).isFailure()){
171 log << MSG::ERROR <<
"Couldn't register "<<name << endreq;
174 sprintf( name,
"EMC_Bhabha_ShowerPhiEastVsEvent");
176 m_HistPhiEast =
new TH1F(name,
title,256, -3.14, 3.14);
177 m_HistPhiEast->GetXaxis()->SetTitle(
"shower phi(radian)");
178 m_HistPhiEast->GetXaxis()->CenterTitle();
179 m_HistPhiEast->GetYaxis()->SetTitle(
"Number of event");
180 m_HistPhiEast->GetYaxis()->CenterTitle();
182 HistName=HistName0+name;
183 if(m_thistsvc->regHist(HistName, m_HistPhiEast).isFailure()){
184 log << MSG::ERROR <<
"Couldn't register"<<name << endreq;
187 sprintf( name,
"EMC_Bhabha_ShowerPhiWestVsEvent");
189 m_HistPhiWest =
new TH1F(name,
title,256, -3.14, 3.14);
190 m_HistPhiWest->GetXaxis()->SetTitle(
"shower phi(radian)");
191 m_HistPhiWest->GetXaxis()->CenterTitle();
192 m_HistPhiWest->GetYaxis()->SetTitle(
"Number of event");
193 m_HistPhiWest->GetYaxis()->CenterTitle();
195 HistName=HistName0+name;
196 if(m_thistsvc->regHist(HistName, m_HistPhiWest).isFailure()){
197 log << MSG::ERROR <<
"Couldn't register" <<name<<endreq;
200 sprintf( name,
"EMC_Bhabha_ShowerThetaPhi");
202 m_ThetaPhi =
new TH2F(name,
"Theta versus Phi",
203 2000, -3.15, 3.15, 2000, 0.1, 3.0);
204 m_ThetaPhi->GetXaxis()->SetTitle(
"shower phi(radian)");
205 m_ThetaPhi->GetXaxis()->CenterTitle();
206 m_ThetaPhi->GetYaxis()->SetTitle(
"shower theta(radian)");
207 m_ThetaPhi->GetYaxis()->CenterTitle();
209 HistName=HistName0+name;
210 if(m_thistsvc->regHist(HistName, m_ThetaPhi).isFailure()){
211 log << MSG::ERROR <<
"Couldn't register" <<name<< endreq;
216 sprintf( name,
"EMC_Bhabha_Time-T0");
218 m_HistTime =
new TH1F(name,
title,100, -40, 60);
219 m_HistTime->GetXaxis()->SetTitle(
"EmcTime-T0 (50ns)");
220 m_HistTime->GetXaxis()->CenterTitle();
221 m_HistTime->GetYaxis()->SetTitle(
"Number of event");
222 m_HistTime->GetYaxis()->CenterTitle();
224 HistName=HistName0+name;
225 if(m_thistsvc->regHist(HistName, m_HistTime).isFailure()){
226 log << MSG::ERROR <<
"Couldn't register" <<name<<endreq;
229 sprintf( name,
"EMC_Bhabha_ixtal");
231 m_HistHitMap =
new TH1F(name,
title,6240, 0, 6240);
232 m_HistHitMap->GetXaxis()->SetTitle(
"ixtalNumber");
233 m_HistHitMap->GetXaxis()->CenterTitle();
234 m_HistHitMap->GetYaxis()->SetTitle(
"Number of event");
235 m_HistHitMap->GetYaxis()->CenterTitle();
237 HistName=HistName0+name;
238 if(m_thistsvc->regHist(HistName, m_HistHitMap).isFailure()){
239 log << MSG::ERROR <<
"Couldn't register" <<name<<endreq;
242 sprintf( name,
"EMC_Bhabha_eSeedvsIxtal ");
244 m_eSeedIxtal =
new TH2F(name,
"eSeed:Ixtal",
245 6240, 0, 6240, 2000, 0, 1.8);
246 m_eSeedIxtal->GetXaxis()->SetTitle(
"Ixtal)");
247 m_eSeedIxtal->GetXaxis()->CenterTitle();
248 m_eSeedIxtal->GetYaxis()->SetTitle(
"eSeed(GeV)");
249 m_eSeedIxtal->GetYaxis()->CenterTitle();
251 HistName=HistName0+name;
252 if(m_thistsvc->regHist(HistName, m_eSeedIxtal).isFailure()){
253 log << MSG::ERROR <<
"Couldn't register" <<name<< endreq;
256 sprintf( name,
"EMC_Bhabha_emcX:emcYeast ");
257 sprintf(
title,
"EMC Bhabha emcX vs emcY of east endcap");
258 m_XYeast =
new TH2F(name,
"emcX:emcY",
259 2000, -100, 100, 2000, -100, 100);
260 m_XYeast->GetXaxis()->SetTitle(
"emcX");
261 m_XYeast->GetXaxis()->CenterTitle();
262 m_XYeast->GetYaxis()->SetTitle(
"emcY");
263 m_XYeast->GetYaxis()->CenterTitle();
265 HistName=HistName0+name;
266 if(m_thistsvc->regHist(HistName, m_XYeast).isFailure()){
267 log << MSG::ERROR <<
"Couldn't register" <<name<< endreq;
270 sprintf( name,
"EMC_Bhabha_emcX:emcYwest ");
271 sprintf(
title,
"EMC Bhabha emcX vs emcY of west endcap");
272 m_XYwest =
new TH2F(name,
"emcX:emcY",
273 2000, -100, 100, 2000, -100, 100);
274 m_XYwest->GetXaxis()->SetTitle(
"emcX");
275 m_XYwest->GetXaxis()->CenterTitle();
276 m_XYwest->GetYaxis()->SetTitle(
"emcY");
277 m_XYwest->GetYaxis()->CenterTitle();
279 HistName=HistName0+name;
280 if(m_thistsvc->regHist(HistName, m_XYwest).isFailure()){
281 log << MSG::ERROR <<
"Couldn't register" <<name<< endreq;
286 scCalib = Gaudi::svcLocator() -> service(
"EmcCalibConstSvc", m_emcCalibConstSvc);
287 if( scCalib != StatusCode::SUCCESS){
288 log << MSG::ERROR <<
"can not use EmcCalibConstSvc" << endreq;
292 log << MSG::INFO <<
"successfully return from initialize()" <<endmsg;
293 return StatusCode::SUCCESS;
301 MsgStream log(
msgSvc(), name());
302 log << MSG::INFO <<
"in execute()" << endreq;
304 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
305 m_runNo=eventHeader->runNumber();
306 m_event=eventHeader->eventNumber();
307 log << MSG::DEBUG <<
"run, evtnum = "
314 SmartDataPtr<DQAEvent::DQAEvent> dqaevt(eventSvc(),
"/Event/DQATag");
316 log << MSG::INFO <<
"success get DQAEvent" << endreq;
318 log << MSG::ERROR <<
"Error accessing DQAEvent" << endreq;
319 return StatusCode::FAILURE;
322 log << MSG::DEBUG <<
"event tag = " << dqaevt->EventTag() << endreq;
325 if ( dqaevt->Bhabha() ) {
326 log << MSG::INFO <<
"Bhabha event" << endreq;
330 double eneShower,theta,phi,
costheta;
335 std::string HistName;
337 for(
int i = 0; i < evtRecEvent->totalCharged(); i++){
340 log << MSG::DEBUG << i <<
" " << (*itTrk)->partId() <<
" "
341 << (*itTrk)->quality() << endreq;
345 if ( ! (*itTrk)->isElectron() )
continue;
357 if ( mdcTrk->
charge() > 0 ) {
358 log << MSG::DEBUG <<
"is electron" << endreq;
360 log << MSG::DEBUG <<
"is positron" << endreq;
363 if((*itTrk)->isEmcShowerValid()) {
366 eneShower=theShower->
energy();
367 theta = theShower->
theta();
368 phi= theShower->
phi();
381 log << MSG::DEBUG <<
"test program " << npart <<
" " << ntheta <<
" " << nphi << endmsg;
386 m_ixtal=m_emcCalibConstSvc->
getIndex(npart,ntheta,nphi);
389 m_emcX = theShower->
x();
390 m_emcY = theShower->
y();
391 m_eSeed = theShower->
eSeed();
392 m_e5x5 = theShower->
e5x5();
393 m_energy = theShower->
energy();
394 m_time = theShower->
time();
403 if (m_thistsvc->getHist(HistName, hmom).isSuccess()) {
406 log << MSG::ERROR <<
"Couldn't retrieve" <<HistName<< endreq;
411 if (m_thistsvc->getHist(HistName, hmom).isSuccess()) {
414 log << MSG::ERROR <<
"Couldn't retrieve" <<HistName<< endreq;
418 m_eSeedIxtal->Fill(m_ixtal,m_eSeed);
425 sprintf(
Name,
"/DQAHist/EMC/EMC_Bhabha_ShowerEneBarrelVsEvent");
427 if (m_thistsvc->getHist(HistName, hmom).isSuccess()) {
428 hmom->Fill(eneShower);
430 log << MSG::ERROR <<
"Couldn't retrieve" <<HistName<< endreq;
433 sprintf(
Name,
"/DQAHist/EMC/EMC_Bhabha_ShowerPhiBarrelVsEvent");
435 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
438 log << MSG::ERROR <<
"Couldn't retrieve" <<HistName<< endreq;
445 sprintf(
Name,
"/DQAHist/EMC/EMC_Bhabha_ShowerEneEastVsEvent");
447 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
448 hmom->Fill( eneShower);
450 log << MSG::ERROR <<
"Couldn't retrieve" <<HistName<< endreq;
453 sprintf(
Name,
"/DQAHist/EMC/EMC_Bhabha_ShowerPhiEastVsEvent");
455 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
458 log << MSG::ERROR <<
"Couldn't retrieve" <<HistName<< endreq;
461 m_XYeast->Fill(m_emcX,m_emcY);
466 sprintf(
Name,
"/DQAHist/EMC/EMC_Bhabha_ShowerEneWestVsEvent");
468 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
469 hmom->Fill( eneShower);
471 log << MSG::ERROR <<
"Couldn't retrieve" <<HistName<< endreq;
474 sprintf(
Name,
"/DQAHist/EMC/EMC_Bhabha_ShowerPhiWestVsEvent");
476 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
479 log << MSG::ERROR <<
"Couldn't retrieve" <<HistName<< endreq;
482 m_XYwest->Fill(m_emcX,m_emcY);
497 sprintf(
Name,
"/DQAHist/EMC/EMC_Bhabha_ShowerThetaVsvent");
499 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
502 log << MSG::ERROR <<
"Couldn't retrieve" <<HistName<< endreq;
505 sprintf(
Name,
"/DQAHist/EMC/EMC_Bhabha_ShowerCosTheta");
507 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
510 log << MSG::ERROR <<
"Couldn't retrieve" <<HistName<< endreq;
521 m_ThetaPhi->Fill(phi, theta);
523 if(m_NtupleOutput==1) m_tuple->write();
529 return StatusCode::SUCCESS;
536 MsgStream log(
msgSvc(), name());
537 log << MSG::INFO <<
"in finalize()" << endmsg;
538 return StatusCode::SUCCESS;
HepGeom::Point3D< double > HepPoint3D
EvtRecTrackCol::iterator EvtRecTrackIterator
double cos(const BesAngle a)
DQA_EMC(const std::string &name, ISvcLocator *pSvcLocator)
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
static unsigned int theta_module(const Identifier &id)
static unsigned int phi_module(const Identifier &id)
virtual int getIndex(unsigned int PartId, unsigned int ThetaIndex, unsigned int PhiIndex) const =0
sprintf(cut,"kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecTrackCol