00001 #include <vector>
00002
00003 #include "GaudiKernel/MsgStream.h"
00004 #include "GaudiKernel/SmartDataPtr.h"
00005 #include "GaudiKernel/PropertyMgr.h"
00006 #include "GaudiKernel/Bootstrap.h"
00007
00008 #include "GaudiKernel/INTupleSvc.h"
00009 #include "GaudiKernel/NTuple.h"
00010 #include "GaudiKernel/ITHistSvc.h"
00011
00012 #include "CLHEP/Vector/ThreeVector.h"
00013 #include "CLHEP/Vector/LorentzVector.h"
00014
00015 #include "EventModel/EventModel.h"
00016 #include "EventModel/Event.h"
00017
00018 #include "EvtRecEvent/EvtRecEvent.h"
00019 #include "EvtRecEvent/EvtRecTrack.h"
00020 #include "DstEvent/TofHitStatus.h"
00021 #include "EventModel/EventHeader.h"
00022
00023
00024 #include "ParticleID/ParticleID.h"
00025
00026 #include "DQAEvent/DQAEvent.h"
00027 #include "DQA_EMC/DQA_EMC.h"
00028
00029 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00030 typedef HepGeom::Point3D<double> HepPoint3D;
00031 #endif
00032 using CLHEP::HepLorentzVector;
00033
00034
00036
00037 DQA_EMC::DQA_EMC(const std::string& name, ISvcLocator* pSvcLocator) :
00038 Algorithm(name, pSvcLocator) {
00039
00040
00041 }
00042
00043
00044 StatusCode DQA_EMC::initialize(){
00045 MsgStream log(msgSvc(), name());
00046
00047 log << MSG::INFO << "in initialize()" << endmsg;
00048 StatusCode status;
00049
00050
00051
00052
00053
00054 NTuplePtr nt(ntupleSvc(), "DQAFILE/EMC");
00055 if ( nt ) m_tuple = nt;
00056 else {
00057 m_tuple = ntupleSvc()->book("DQAFILE/EMC", CLID_ColumnWiseTuple, "EMC ntuple");
00058 if( m_tuple ) {
00059 status = m_tuple->addItem("ixtal", m_ixtal);
00060 status = m_tuple->addItem("npart", m_npart);
00061 status = m_tuple->addItem("ntheta",m_ntheta);
00062 status = m_tuple->addItem("nphi", m_nphi);
00063 status = m_tuple->addItem("theta", m_theta);
00064 status = m_tuple->addItem("phi", m_phi);
00065 status = m_tuple->addItem("emcX", m_emcX);
00066 status = m_tuple->addItem("emcY", m_emcY);
00067 status = m_tuple->addItem("eSeed", m_eSeed);
00068 status = m_tuple->addItem("e5x5", m_e5x5);
00069 status = m_tuple->addItem("energy",m_energy);
00070 status = m_tuple->addItem("time", m_time);
00071
00072 } else {
00073 log << MSG::ERROR << "Can not book N-tuple:" << long(m_tuple) << endreq;
00074 }
00075 }
00076
00077
00078 if(service("THistSvc", m_thistsvc).isFailure()) {
00079 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
00080 return StatusCode::FAILURE;
00081 }
00082
00083 char name[60];
00084 char title[60];
00085
00086 std::string HistName0, HistName;
00087 HistName0="/DQAHist/EMC/";
00088
00089 sprintf( name, "EMC_Bhabha_ShowerEneBarrelVsEvent");
00090 sprintf( title, "EMC Bhabha ShowerEneBarrel vs event");
00091
00092 m_HistEnergyB =new TH1F(name,title,200,0.,2.0);
00093 m_HistEnergyB->GetXaxis()->SetTitle("shower energy(GeV)");
00094 m_HistEnergyB->GetXaxis()->CenterTitle();
00095 m_HistEnergyB->GetYaxis()->SetTitle("Number of event");
00096 m_HistEnergyB->GetYaxis()->CenterTitle();
00097
00098 HistName=HistName0+name;
00099 if(m_thistsvc->regHist(HistName, m_HistEnergyB).isFailure()){
00100 log << MSG::ERROR << "Couldn't register " <<name<< endreq;
00101 }
00102
00103 sprintf( name, "EMC_Bhabha_ShowerEneEastVsEvent");
00104 sprintf( title, "EMC Bhabha ShowerEneEast vs event");
00105 m_HistEnergyEast =new TH1F(name,title,200,0.,2.0);
00106 m_HistEnergyEast->GetXaxis()->SetTitle("shower energy(GeV)");
00107 m_HistEnergyEast->GetXaxis()->CenterTitle();
00108 m_HistEnergyEast->GetYaxis()->SetTitle("Number of event");
00109 m_HistEnergyEast->GetYaxis()->CenterTitle();
00110
00111 HistName=HistName0+name;
00112 if(m_thistsvc->regHist(HistName, m_HistEnergyEast).isFailure()){
00113 log << MSG::ERROR << "Couldn't register "<<name << endreq;
00114 }
00115
00116
00117 sprintf( name, "EMC_Bhabha_ShowerEneWestVsEvent");
00118 sprintf( title, "EMC Bhabha ShowerEneWest vs event");
00119 m_HistEnergyWest =new TH1F(name,title,200,0.,2.0);
00120 m_HistEnergyWest->GetXaxis()->SetTitle("shower energy(GeV)");
00121 m_HistEnergyWest->GetXaxis()->CenterTitle();
00122 m_HistEnergyWest->GetYaxis()->SetTitle("Number of event");
00123 m_HistEnergyWest->GetYaxis()->CenterTitle();
00124
00125
00126 HistName=HistName0+name;
00127 if( m_thistsvc->regHist(HistName, m_HistEnergyWest).isFailure()){
00128 log << MSG::ERROR << "Couldn't register "<<name << endreq;
00129 }
00130
00131
00132 sprintf( name, "EMC_Bhabha_ShowerThetaVsvent");
00133 sprintf( title, "EMC Bhabha ShowerTheta vs event");
00134 m_HistTheta =new TH1F(name,title,56, 0, 56);
00135 m_HistTheta->GetXaxis()->SetTitle("shower ID(theta)");
00136 m_HistTheta->GetXaxis()->CenterTitle();
00137 m_HistTheta->GetYaxis()->SetTitle("Number of event");
00138 m_HistTheta->GetYaxis()->CenterTitle();
00139
00140 HistName=HistName0+name;
00141 if( m_thistsvc->regHist(HistName, m_HistTheta).isFailure()){
00142 log << MSG::ERROR << "Couldn't register" <<name<< endreq;
00143 }
00144
00145 sprintf(name,"EMC_Bhabha_ShowerCosTheta");
00146 sprintf(title,"Emc Bhabha Costheta");
00147 m_HistCosTheta =new TH1F(name,title,200, -1.0, 1.0);
00148 m_HistCosTheta->GetXaxis()->SetTitle("shower cos(theta)");
00149 m_HistCosTheta->GetXaxis()->CenterTitle();
00150 m_HistCosTheta->GetYaxis()->SetTitle("Number of event");
00151 m_HistCosTheta->GetYaxis()->CenterTitle();
00152
00153 HistName=HistName0+name;
00154 if(m_thistsvc->regHist(HistName, m_HistCosTheta).isFailure()){
00155 log << MSG::ERROR << "Couldn't register "<<name << endreq;
00156 }
00157
00158
00159 sprintf( name, "EMC_Bhabha_ShowerPhiBarrelVsEvent");
00160 sprintf( title, "EMC Bhabha ShowerPhiBarrel vs event");
00161 m_HistPhiB =new TH1F(name,title,120, 0, 120);
00162 m_HistPhiB->GetXaxis()->SetTitle("shower ID(phi)");
00163 m_HistPhiB->GetXaxis()->CenterTitle();
00164 m_HistPhiB->GetYaxis()->SetTitle("Number of event");
00165 m_HistPhiB->GetYaxis()->CenterTitle();
00166
00167
00168 HistName=HistName0+name;
00169 if( m_thistsvc->regHist(HistName, m_HistPhiB).isFailure()){
00170 log << MSG::ERROR << "Couldn't register "<<name << endreq;
00171 }
00172
00173 sprintf( name, "EMC_Bhabha_ShowerPhiEastVsEvent");
00174 sprintf( title, "EMC Bhabha ShowerPhiEast vs event");
00175 m_HistPhiEast =new TH1F(name,title,256, -3.14, 3.14);
00176 m_HistPhiEast->GetXaxis()->SetTitle("shower phi(radian)");
00177 m_HistPhiEast->GetXaxis()->CenterTitle();
00178 m_HistPhiEast->GetYaxis()->SetTitle("Number of event");
00179 m_HistPhiEast->GetYaxis()->CenterTitle();
00180
00181 HistName=HistName0+name;
00182 if(m_thistsvc->regHist(HistName, m_HistPhiEast).isFailure()){
00183 log << MSG::ERROR << "Couldn't register"<<name << endreq;
00184 }
00185
00186 sprintf( name, "EMC_Bhabha_ShowerPhiWestVsEvent");
00187 sprintf( title, "EMC Bhabha ShowerPhiWest vs event");
00188 m_HistPhiWest =new TH1F(name,title,256, -3.14, 3.14);
00189 m_HistPhiWest->GetXaxis()->SetTitle("shower phi(radian)");
00190 m_HistPhiWest->GetXaxis()->CenterTitle();
00191 m_HistPhiWest->GetYaxis()->SetTitle("Number of event");
00192 m_HistPhiWest->GetYaxis()->CenterTitle();
00193
00194 HistName=HistName0+name;
00195 if(m_thistsvc->regHist(HistName, m_HistPhiWest).isFailure()){
00196 log << MSG::ERROR << "Couldn't register" <<name<<endreq;
00197 }
00198
00199 sprintf( name, "EMC_Bhabha_ShowerThetaPhi");
00200 sprintf( title, "EMC Bhabha ShowerThetaPhi");
00201 m_ThetaPhi =new TH2F(name,"Theta versus Phi",
00202 2000, -3.15, 3.15, 2000, 0.1, 3.0);
00203 m_ThetaPhi->GetXaxis()->SetTitle("shower phi(radian)");
00204 m_ThetaPhi->GetXaxis()->CenterTitle();
00205 m_ThetaPhi->GetYaxis()->SetTitle("shower theta(radian)");
00206 m_ThetaPhi->GetYaxis()->CenterTitle();
00207
00208 HistName=HistName0+name;
00209 if(m_thistsvc->regHist(HistName, m_ThetaPhi).isFailure()){
00210 log << MSG::ERROR << "Couldn't register" <<name<< endreq;
00211 }
00212
00213
00215 sprintf( name, "EMC_Bhabha_Time-T0");
00216 sprintf( title, "EMC Bhabha Time-T0 distribution");
00217 m_HistTime =new TH1F(name,title,100, -40, 60);
00218 m_HistTime->GetXaxis()->SetTitle("EmcTime-T0 (50ns)");
00219 m_HistTime->GetXaxis()->CenterTitle();
00220 m_HistTime->GetYaxis()->SetTitle("Number of event");
00221 m_HistTime->GetYaxis()->CenterTitle();
00222
00223 HistName=HistName0+name;
00224 if(m_thistsvc->regHist(HistName, m_HistTime).isFailure()){
00225 log << MSG::ERROR << "Couldn't register" <<name<<endreq;
00226 }
00227
00228 sprintf( name, "EMC_Bhabha_ixtal");
00229 sprintf( title, "EMC Bhabha ixtal distribution");
00230 m_HistHitMap =new TH1F(name,title,6240, 0, 6240);
00231 m_HistHitMap->GetXaxis()->SetTitle("ixtalNumber");
00232 m_HistHitMap->GetXaxis()->CenterTitle();
00233 m_HistHitMap->GetYaxis()->SetTitle("Number of event");
00234 m_HistHitMap->GetYaxis()->CenterTitle();
00235
00236 HistName=HistName0+name;
00237 if(m_thistsvc->regHist(HistName, m_HistHitMap).isFailure()){
00238 log << MSG::ERROR << "Couldn't register" <<name<<endreq;
00239 }
00240
00241 sprintf( name, "EMC_Bhabha_eSeedvsIxtal ");
00242 sprintf( title, "EMC Bhabha eSeed vs ixtal");
00243 m_eSeedIxtal =new TH2F(name,"eSeed:Ixtal",
00244 6240, 0, 6240, 2000, 0, 1.8);
00245 m_eSeedIxtal->GetXaxis()->SetTitle("Ixtal)");
00246 m_eSeedIxtal->GetXaxis()->CenterTitle();
00247 m_eSeedIxtal->GetYaxis()->SetTitle("eSeed(GeV)");
00248 m_eSeedIxtal->GetYaxis()->CenterTitle();
00249
00250 HistName=HistName0+name;
00251 if(m_thistsvc->regHist(HistName, m_eSeedIxtal).isFailure()){
00252 log << MSG::ERROR << "Couldn't register" <<name<< endreq;
00253 }
00254
00255 sprintf( name, "EMC_Bhabha_emcX:emcYeast ");
00256 sprintf( title, "EMC Bhabha emcX vs emcY of east endcap");
00257 m_XYeast =new TH2F(name,"emcX:emcY",
00258 2000, -100, 100, 2000, -100, 100);
00259 m_XYeast->GetXaxis()->SetTitle("emcX");
00260 m_XYeast->GetXaxis()->CenterTitle();
00261 m_XYeast->GetYaxis()->SetTitle("emcY");
00262 m_XYeast->GetYaxis()->CenterTitle();
00263
00264 HistName=HistName0+name;
00265 if(m_thistsvc->regHist(HistName, m_XYeast).isFailure()){
00266 log << MSG::ERROR << "Couldn't register" <<name<< endreq;
00267 }
00268
00269 sprintf( name, "EMC_Bhabha_emcX:emcYwest ");
00270 sprintf( title, "EMC Bhabha emcX vs emcY of west endcap");
00271 m_XYwest =new TH2F(name,"emcX:emcY",
00272 2000, -100, 100, 2000, -100, 100);
00273 m_XYwest->GetXaxis()->SetTitle("emcX");
00274 m_XYwest->GetXaxis()->CenterTitle();
00275 m_XYwest->GetYaxis()->SetTitle("emcY");
00276 m_XYwest->GetYaxis()->CenterTitle();
00277
00278 HistName=HistName0+name;
00279 if(m_thistsvc->regHist(HistName, m_XYwest).isFailure()){
00280 log << MSG::ERROR << "Couldn't register" <<name<< endreq;
00281 }
00282
00283
00284 StatusCode scCalib;
00285 scCalib = Gaudi::svcLocator() -> service("EmcCalibConstSvc", m_emcCalibConstSvc);
00286 if( scCalib != StatusCode::SUCCESS){
00287 log << MSG::ERROR << "can not use EmcCalibConstSvc" << endreq;
00288 }
00289
00290
00291 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00292 return StatusCode::SUCCESS;
00293
00294
00295 }
00296
00297
00298 StatusCode DQA_EMC::execute() {
00299
00300 MsgStream log(msgSvc(), name());
00301 log << MSG::INFO << "in execute()" << endreq;
00302
00303 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00304 m_runNo=eventHeader->runNumber();
00305 m_event=eventHeader->eventNumber();
00306 log << MSG::DEBUG <<"run, evtnum = "
00307 << m_runNo << " , "
00308 << m_event <<endreq;
00309
00310
00311 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00312
00313 SmartDataPtr<DQAEvent::DQAEvent> dqaevt(eventSvc(), "/Event/DQATag");
00314 if ( dqaevt ) {
00315 log << MSG::INFO << "success get DQAEvent" << endreq;
00316 } else {
00317 log << MSG::ERROR << "Error accessing DQAEvent" << endreq;
00318 return StatusCode::FAILURE;
00319 }
00320
00321 log << MSG::DEBUG << "event tag = " << dqaevt->EventTag() << endreq;
00322
00323
00324 if ( dqaevt->Bhabha() ) {
00325 log << MSG::INFO << "Bhabha event" << endreq;
00326 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00327
00328
00329 double eneShower,theta,phi,costheta;
00330 RecEmcID showerId;
00331 unsigned int npart;
00332 int ntheta,nphi;
00333 char Name[60];
00334 std::string HistName;
00335
00336 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00337 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00338
00339 log << MSG::DEBUG << i << " " << (*itTrk)->partId() << " "
00340 << (*itTrk)->quality() << endreq;
00341
00342
00343 if ( (*itTrk)->partId() != 1 ) continue;
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00355 if ( mdcTrk->charge() > 0 ) {
00356 log << MSG::DEBUG << "is electron" << endreq;
00357 } else {
00358 log << MSG::DEBUG << "is positron" << endreq;
00359 }
00360
00361 if((*itTrk)->isEmcShowerValid()) {
00362
00363 RecEmcShower *theShower = (*itTrk)->emcShower();
00364 eneShower=theShower->energy();
00365 theta = theShower->theta();
00366 phi= theShower->phi();
00367 costheta=cos(theta);
00368 showerId = theShower->getShowerId();
00369
00370 npart = EmcID::barrel_ec(showerId);
00371
00372
00373
00374
00375 ntheta=EmcID::theta_module(showerId);
00376
00377 nphi=EmcID::phi_module(showerId);
00378
00379 m_npart = npart;
00380 m_ntheta = ntheta;
00381 m_nphi = nphi;
00382 m_ixtal=m_emcCalibConstSvc->getIndex(npart,ntheta,nphi);
00383 m_theta = theta;
00384 m_phi = phi;
00385 m_emcX = theShower->x();
00386 m_emcY = theShower->y();
00387 m_eSeed = theShower->eSeed();
00388 m_e5x5 = theShower->e5x5();
00389 m_energy = theShower->energy();
00390 m_time = theShower->time();
00391
00392 TH1 *hmom(0);
00393
00394
00395
00396 sprintf( Name, "/DQAHist/EMC/EMC_Bhabha_Time-T0");
00397 HistName=Name;
00398 if (m_thistsvc->getHist(HistName, hmom).isSuccess()) {
00399 hmom->Fill(m_time);
00400 } else {
00401 log << MSG::ERROR << "Couldn't retrieve" <<HistName<< endreq;
00402 }
00403
00404 sprintf( Name, "/DQAHist/EMC/EMC_Bhabha_ixtal");
00405 HistName=Name;
00406 if (m_thistsvc->getHist(HistName, hmom).isSuccess()) {
00407 hmom->Fill(m_ixtal);
00408 } else {
00409 log << MSG::ERROR << "Couldn't retrieve" <<HistName<< endreq;
00410 }
00411
00412
00413 m_eSeedIxtal->Fill(m_ixtal,m_eSeed);
00414
00415
00416
00417
00418 if (npart==1){
00419
00420 sprintf( Name, "/DQAHist/EMC/EMC_Bhabha_ShowerEneBarrelVsEvent");
00421 HistName=Name;
00422 if (m_thistsvc->getHist(HistName, hmom).isSuccess()) {
00423 hmom->Fill(eneShower);
00424 } else {
00425 log << MSG::ERROR << "Couldn't retrieve" <<HistName<< endreq;
00426 }
00427
00428 sprintf( Name, "/DQAHist/EMC/EMC_Bhabha_ShowerPhiBarrelVsEvent");
00429 HistName=Name;
00430 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
00431 hmom->Fill(nphi);
00432 } else {
00433 log << MSG::ERROR << "Couldn't retrieve" <<HistName<< endreq;
00434 }
00435
00436 }
00437
00438 if (npart==0){
00439
00440 sprintf( Name, "/DQAHist/EMC/EMC_Bhabha_ShowerEneEastVsEvent");
00441 HistName=Name;
00442 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
00443 hmom->Fill( eneShower);
00444 } else {
00445 log << MSG::ERROR << "Couldn't retrieve" <<HistName<< endreq;
00446 }
00447
00448 sprintf( Name, "/DQAHist/EMC/EMC_Bhabha_ShowerPhiEastVsEvent");
00449 HistName=Name;
00450 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
00451 hmom->Fill(phi);
00452 } else {
00453 log << MSG::ERROR << "Couldn't retrieve" <<HistName<< endreq;
00454 }
00455
00456 m_XYeast->Fill(m_emcX,m_emcY);
00457
00458 }
00459 if (npart==2){
00460
00461 sprintf( Name, "/DQAHist/EMC/EMC_Bhabha_ShowerEneWestVsEvent");
00462 HistName=Name;
00463 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
00464 hmom->Fill( eneShower);
00465 } else {
00466 log << MSG::ERROR << "Couldn't retrieve" <<HistName<< endreq;
00467 }
00468
00469 sprintf( Name, "/DQAHist/EMC/EMC_Bhabha_ShowerPhiWestVsEvent");
00470 HistName=Name;
00471 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
00472 hmom->Fill(phi);
00473 } else {
00474 log << MSG::ERROR << "Couldn't retrieve" <<HistName<< endreq;
00475 }
00476
00477 m_XYwest->Fill(m_emcX,m_emcY);
00478
00479 }
00480
00481 int m_nthe=-9;
00482 if (npart==0){
00483 m_nthe = ntheta;
00484 }
00485 if (npart==2){
00486 m_nthe = 55-ntheta;
00487 }
00488 if (npart==1){
00489 m_nthe = ntheta+6;
00490 }
00491
00492 sprintf( Name, "/DQAHist/EMC/EMC_Bhabha_ShowerThetaVsvent");
00493 HistName=Name;
00494 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
00495 hmom->Fill(m_nthe);
00496 } else {
00497 log << MSG::ERROR << "Couldn't retrieve" <<HistName<< endreq;
00498 }
00499
00500 sprintf( Name,"/DQAHist/EMC/EMC_Bhabha_ShowerCosTheta");
00501 HistName=Name;
00502 if (m_thistsvc->getHist(HistName, hmom ).isSuccess()) {
00503 hmom->Fill(costheta);
00504 } else {
00505 log << MSG::ERROR << "Couldn't retrieve" <<HistName<< endreq;
00506 }
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516 m_ThetaPhi->Fill(phi, theta);
00517
00518 m_tuple->write();
00519
00520 }
00521 }
00522 }
00523
00524 return StatusCode::SUCCESS;
00525
00526 }
00527
00528
00529 StatusCode DQA_EMC::finalize() {
00530
00531 MsgStream log(msgSvc(), name());
00532 log << MSG::INFO << "in finalize()" << endmsg;
00533 return StatusCode::SUCCESS;
00534 }
00535
00536