BOSS_7.0.2: /home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/DQA/DQA_MDC/DQA_MDC-00-00-03/src/DQA

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 #include "GaudiKernel/Algorithm.h"
00008 
00009 #include "GaudiKernel/INTupleSvc.h"
00010 #include "GaudiKernel/NTuple.h"
00011 #include "GaudiKernel/ITHistSvc.h"
00012 
00013 #include "CLHEP/Vector/ThreeVector.h"
00014 #include "CLHEP/Vector/LorentzVector.h"
00015 
00016 #include "EventModel/EventModel.h"
00017 #include "EventModel/Event.h"
00018 #include "RawEvent/RawDataUtil.h"
00019 
00020 #include "EvtRecEvent/EvtRecEvent.h"
00021 #include "EvtRecEvent/EvtRecTrack.h"
00022 #include "DstEvent/TofHitStatus.h"
00023 #include "EventModel/EventHeader.h"
00024 
00025 #include "MdcRecEvent/RecMdcTrack.h"
00026 #include "MdcRecEvent/RecMdcKalTrack.h"
00027 #include "MdcRecEvent/RecMdcHit.h"
00028 #include "ReconEvent/ReconEvent.h"
00029 
00030 #include "DQAEvent/DQAEvent.h"
00031 #include "DQA_MDC/DQA_MDC.h"
00032 
00033 using CLHEP::Hep3Vector;
00035 
00036 DQA_MDC::DQA_MDC(const std::string& name, ISvcLocator* pSvcLocator) :
00037     Algorithm(name, pSvcLocator) {
00038 
00039         //Declare the properties  
00040     }
00041 
00042 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
00043 StatusCode DQA_MDC::initialize(){
00044     MsgStream log(msgSvc(), name());
00045 
00046     log << MSG::INFO << "in initialize()" << endmsg;
00047     StatusCode status;
00048 
00049     // Call Histogram service
00050     if(service("THistSvc", m_thsvc).isFailure()) {
00051         log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
00052         return StatusCode::FAILURE;
00053     }
00054 
00055     // Resolution histograms
00056     m_hresAllIncBb = new TH1F("HResAllIncBb", "", 200, -1.0, 1.0);
00057     if(m_thsvc->regHist("/DQAHist/MDC/hresAllIncBb", m_hresAllIncBb).isFailure()) 
00058     { log << MSG::ERROR << "Couldn't register HResAllIncBb" << endreq; }
00059 
00060     m_hresAllExcBb = new TH1F("HResAllExcBb", "", 200, -1.0, 1.0);
00061     if(m_thsvc->regHist("/DQAHist/MDC/hresAllExcBb", m_hresAllExcBb).isFailure()) 
00062     { log << MSG::ERROR << "Couldn't register HResAllExcBb" << endreq; }
00063  
00064     m_hresAllEvaBb = new TH1F("HResAllEvaBb", "", 200, -1.0, 1.0);
00065     if(m_thsvc->regHist("/DQAHist/MDC/hresAllEvaBb", m_hresAllEvaBb).isFailure()) 
00066     { log << MSG::ERROR << "Couldn't register HResAllEvaBb" << endreq; }
00067  
00068     m_ppLabBb = new TH1F("PpLabBb", "", 800, 0, 3);
00069     if(m_thsvc->regHist("/DQAHist/MDC/hppLabBb", m_ppLabBb).isFailure()) 
00070     { log << MSG::ERROR << "Couldn't register PpLabBb" << endreq; }
00071     m_pmLabBb = new TH1F("PmLabBb", "", 800, 0, 3);
00072     if(m_thsvc->regHist("/DQAHist/MDC/hpmLabBb", m_pmLabBb).isFailure()) 
00073     { log << MSG::ERROR << "Couldn't register PmLabBb" << endreq; }
00074 
00075     m_ppCmsBb = new TH1F("PpCmsBb", "", 800, 0, 3);
00076     if(m_thsvc->regHist("/DQAHist/MDC/hppCmsBb", m_ppCmsBb).isFailure()) 
00077     { log << MSG::ERROR << "Couldn't register PpCmsBb" << endreq; }
00078     m_pmCmsBb = new TH1F("PmCmsBb", "", 800, 0, 3);
00079     if(m_thsvc->regHist("/DQAHist/MDC/hpmCmsBb", m_pmCmsBb).isFailure()) 
00080     { log << MSG::ERROR << "Couldn't register PmCmsBb" << endreq; }
00081 
00082     m_pTotLabBb = new TH1F("PTotLabBb", "", 800, 0, 3);
00083     if(m_thsvc->regHist("/DQAHist/MDC/hpTotLabBb", m_pTotLabBb).isFailure()) 
00084     { log << MSG::ERROR << "Couldn't register PTotLabBb" << endreq; }
00085     m_pTotCmsBb = new TH1F("PTotCmsBb", "", 800, 0, 3);
00086     if(m_thsvc->regHist("/DQAHist/MDC/hpTotCmsBb", m_pTotCmsBb).isFailure()) 
00087     { log << MSG::ERROR << "Couldn't register PTotCmsBb" << endreq; }
00088 
00089     m_hresAllIncDimu = new TH1F("HResAllIncDimu", "", 200, -1.0, 1.0);
00090     if(m_thsvc->regHist("/DQAHist/MDC/hresAllIncDimu", m_hresAllIncDimu).isFailure()) 
00091     { log << MSG::ERROR << "Couldn't register HResAllIncDimu" << endreq; }
00092 
00093     m_hresAllExcDimu = new TH1F("HResAllExcDimu", "", 200, -1.0, 1.0);
00094     if(m_thsvc->regHist("/DQAHist/MDC/hresAllExcDimu", m_hresAllExcDimu).isFailure()) 
00095     { log << MSG::ERROR << "Couldn't register HResAllExcDimu" << endreq; }
00096  
00097     m_hresAllEvaDimu = new TH1F("HResAllEvaDimu", "", 200, -1.0, 1.0);
00098     if(m_thsvc->regHist("/DQAHist/MDC/hresAllEvaDimu", m_hresAllEvaDimu).isFailure()) 
00099     { log << MSG::ERROR << "Couldn't register HResAllEvaDimu" << endreq; }
00100  
00101     m_ppLabDimu = new TH1F("PpLabDimu", "", 800, 0, 3);
00102     if(m_thsvc->regHist("/DQAHist/MDC/hppLabDimu", m_ppLabDimu).isFailure()) 
00103     { log << MSG::ERROR << "Couldn't register PpLabDimu" << endreq; }
00104     m_pmLabDimu = new TH1F("PmLabDimu", "", 800, 0, 3);
00105     if(m_thsvc->regHist("/DQAHist/MDC/hpmLabDimu", m_pmLabDimu).isFailure()) 
00106     { log << MSG::ERROR << "Couldn't register PmLabDimu" << endreq; }
00107 
00108     m_ppCmsDimu = new TH1F("PpCmsDimu", "", 800, 0, 3);
00109     if(m_thsvc->regHist("/DQAHist/MDC/hppCmsDimu", m_ppCmsDimu).isFailure()) 
00110     { log << MSG::ERROR << "Couldn't register PpCmsDimu" << endreq; }
00111     m_pmCmsDimu = new TH1F("PmCmsDimu", "", 800, 0, 3);
00112     if(m_thsvc->regHist("/DQAHist/MDC/hpmCmsDimu", m_pmCmsDimu).isFailure()) 
00113     { log << MSG::ERROR << "Couldn't register PmCmsDimu" << endreq; }
00114 
00115     m_pTotLabDimu = new TH1F("PTotLabDimu", "", 800, 0, 3);
00116     if(m_thsvc->regHist("/DQAHist/MDC/hpTotLabDimu", m_pTotLabDimu).isFailure()) 
00117     { log << MSG::ERROR << "Couldn't register PTotLabDimu" << endreq; }
00118     m_pTotCmsDimu = new TH1F("PTotCmsDimu", "", 800, 0, 3);
00119     if(m_thsvc->regHist("/DQAHist/MDC/hpTotCmsDimu", m_pTotCmsDimu).isFailure()) 
00120     { log << MSG::ERROR << "Couldn't register PTotCmsDimu" << endreq; }
00121 
00122 
00123     log << MSG::INFO << "Initialize done!" <<endmsg;
00124     return StatusCode::SUCCESS;
00125 }
00126 
00127 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
00128 StatusCode DQA_MDC::execute() {
00129 
00130     MsgStream log(msgSvc(), name());
00131     log << MSG::INFO << "in execute()" << endreq;
00132 
00133     SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00134     m_run   = eventHeader->runNumber();
00135     m_event = eventHeader->eventNumber();
00136     log << MSG::DEBUG <<"Run " << m_run << "\tEvent "  << m_event <<endreq;
00137 
00138     SmartDataPtr<DQAEvent::DQAEvent> dqaevt(eventSvc(), "/Event/DQATag");
00139     if ( dqaevt ) {
00140         log << MSG::INFO << "success get DQAEvent" << endreq;
00141     } else {
00142         log << MSG::ERROR << "Error accessing DQAEvent" << endreq;
00143         return StatusCode::FAILURE;
00144     }
00145     log << MSG::DEBUG << "DQA event tag = " << dqaevt->EventTag() << endreq;
00146 
00147     double p;
00148     double p_cms;
00149     bool isDimu = false;
00150     if ( dqaevt->Dimu() )  isDimu = true;
00151     log << MSG::INFO << "DQADimuTag:\t" << dqaevt->Dimu() << endreq;
00152     //cout << "DQADimuTag:\t" << dqaevt->Dimu() << endl;
00153 
00154     bool isBb = false;
00155     if ( dqaevt->Bhabha() )  isBb= true;
00156     log << MSG::INFO << "DQABbTag:\t" << dqaevt->Bhabha() << endreq;
00157     //cout << "DQABbTag:\t" << dqaevt->Bhabha() << endl;
00158 
00159     SmartDataPtr<RecMdcKalTrackCol> kaltrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
00160     if (!kaltrkCol) {
00161         log << MSG::FATAL << "Could not find RecMdcKalTrackCol" << endreq;
00162         return StatusCode::FAILURE;
00163     }
00164     RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
00165     for(; iter_trk != kaltrkCol->end(); iter_trk++) {
00166         if (isBb) RecMdcKalTrack::setPidType(RecMdcKalTrack::electron);
00167         else if (isDimu) RecMdcKalTrack::setPidType(RecMdcKalTrack::muon);
00168         //else return StatusCode::SUCCESS;
00169 
00170         double dr = (*iter_trk)->dr();
00171         double phi0 = (*iter_trk)->fi0();
00172         double kappa = (*iter_trk)->kappa();
00173         double dz = (*iter_trk)->dz();
00174         double tanl = (*iter_trk)->tanl();
00175         double pt = 1.0 / kappa;
00176         HepLorentzVector p4;
00177         p4.setPx(- sin(phi0) / fabs(kappa));
00178         p4.setPy(cos(phi0) / fabs(kappa));
00179         p4.setPz(tanl / fabs(kappa));
00180 
00181         double p3 = p4.mag();
00182         double mass;
00183         if(isBb) mass = 0.000511;
00184         else if(isDimu) mass = 0.105658;
00185 
00186         p4.setE(sqrt(p3 * p3 + mass * mass));
00187         p = p4.rho();
00188 
00189         double ecm = 3.68632;
00190         HepLorentzVector psip(0.011 * ecm, 0, 0.0075, ecm);
00191         Hep3Vector boostv = psip.boostVector();
00192         p4.boost(- boostv);
00193         p_cms = p4.rho();
00194         TH1 *hmom(0);
00195         if (isBb){
00196             if (m_thsvc->getHist("/DQAHist/MDC/hpTotLabBb", hmom).isSuccess()) {
00197                 hmom->Fill(p);
00198             } else {
00199                 log << MSG::ERROR << "Couldn't retrieve hpTotLabBb" << endreq;
00200             }
00201             if (m_thsvc->getHist("/DQAHist/MDC/hpTotCmsBb", hmom).isSuccess()) {
00202                 hmom->Fill(p_cms);
00203             } else {
00204                 log << MSG::ERROR << "Couldn't retrieve hpTotCmsBb" << endreq;
00205             }
00206             if (pt > 0) {
00207                 if (m_thsvc->getHist("/DQAHist/MDC/hppLabBb", hmom).isSuccess()) {
00208                     hmom->Fill(p);
00209                 } else {
00210                     log << MSG::ERROR << "Couldn't retrieve hppLabBb" << endreq;
00211                 }
00212                 if (m_thsvc->getHist("/DQAHist/MDC/hppCmsBb", hmom).isSuccess()) {
00213                     hmom->Fill(p);
00214                 } else {
00215                     log << MSG::ERROR << "Couldn't retrieve hppCmsBb" << endreq;
00216                 }
00217             } else if (pt < 0) {
00218                 if (m_thsvc->getHist("/DQAHist/MDC/hpmLabBb", hmom).isSuccess()) {
00219                     hmom->Fill(p);
00220                 } else {
00221                     log << MSG::ERROR << "Couldn't retrieve hpmLabBb" << endreq;
00222                 }
00223                 if (m_thsvc->getHist("/DQAHist/MDC/hpmCmsBb", hmom).isSuccess()) {
00224                     hmom->Fill(p);
00225                 } else {
00226                     log << MSG::ERROR << "Couldn't retrieve hpmCmsBb" << endreq;
00227                 }
00228             }
00229         } else if (isDimu){
00230             if (m_thsvc->getHist("/DQAHist/MDC/hpTotLabDimu", hmom).isSuccess()) {
00231                 hmom->Fill(p);
00232             } else {
00233                 log << MSG::ERROR << "Couldn't retrieve hpTotLabDimu" << endreq;
00234             }
00235             if (m_thsvc->getHist("/DQAHist/MDC/hpTotCmsDimu", hmom).isSuccess()) {
00236                 hmom->Fill(p_cms);
00237             } else {
00238                 log << MSG::ERROR << "Couldn't retrieve hpTotCmsDimu" << endreq;
00239             }
00240             if (pt > 0) {
00241                 if (m_thsvc->getHist("/DQAHist/MDC/hppLabDimu", hmom).isSuccess()) {
00242                     hmom->Fill(p);
00243                 } else {
00244                     log << MSG::ERROR << "Couldn't retrieve hppLabDimu" << endreq;
00245                 }
00246                 if (m_thsvc->getHist("/DQAHist/MDC/hppCmsDimu", hmom).isSuccess()) {
00247                     hmom->Fill(p);
00248                 } else {
00249                     log << MSG::ERROR << "Couldn't retrieve hppCmsDimu" << endreq;
00250                 }
00251             } else if (pt < 0) {
00252                 if (m_thsvc->getHist("/DQAHist/MDC/hpmLabDimu", hmom).isSuccess()) {
00253                     hmom->Fill(p);
00254                 } else {
00255                     log << MSG::ERROR << "Couldn't retrieve hpmLabDimu" << endreq;
00256                 }
00257                 if (m_thsvc->getHist("/DQAHist/MDC/hpmCmsDimu", hmom).isSuccess()) {
00258                     hmom->Fill(p);
00259                 } else {
00260                     log << MSG::ERROR << "Couldn't retrieve hpmCmsDimu" << endreq;
00261                 }
00262             }
00263         }
00264 
00265     }
00266 
00267     SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
00268     if(!newtrkCol){
00269         log << MSG::ERROR << "Could not find RecMdcTrackCol" << endreq;
00270         return ( StatusCode::FAILURE );
00271     }
00272 
00273     RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
00274     for(; it_trk != newtrkCol->end(); it_trk++){
00275         HitRefVec gothits = (*it_trk) -> getVecHits();
00276         HitRefVec::iterator it_hit = gothits.begin();
00277         for(; it_hit != gothits.end(); it_hit++){
00278 
00279             //cout << __LINE__ << endl;
00280             //cout << __LINE__ << endl;
00281             double lr = (*it_hit) -> getFlagLR();
00282             if(-1 == lr) lr = 0;        // definition not same as MdcRecHit
00283             double resilrInc;
00284             double resilrExc;
00285             //for boss before 650
00286             double docaInc = (*it_hit) -> getDoca();
00287             double docaExc = docaInc;
00288             //boss650
00289             //double docaInc = (*it_hit) -> getDocaIncl();
00290             //double docaExc = (*it_hit) -> getDocaExcl();
00291             double dmeas;
00292             if(0 == lr){
00293                 dmeas = (*it_hit) -> getDriftDistLeft();
00294             }else {
00295                 dmeas = (*it_hit) -> getDriftDistRight();
00296             }
00297             // the following is for cm to mm
00298             double df = 10.0;
00299             docaInc *= df;
00300             docaExc *= df;
00301             dmeas *= df;
00302             double resiInc = fabs(dmeas) - fabs(docaInc);
00303             if( 0 == lr )     resilrInc = -1.0 * resiInc;
00304             else              resilrInc = resiInc;
00305 
00306             double resiExc = fabs(dmeas) - fabs(docaExc);
00307             if( 0 == lr )     resilrExc = -1.0 * resiExc;
00308             else              resilrExc = resiExc;
00309 
00310             double resilrEva = 0.5 * (resilrInc + resilrExc);
00311             //cout << "resilrInc = " << resilrInc << ", resilrExc = " << resilrExc << ", resilrEva = " << resilrEva << endl;
00312             TH1 *htmp(0);
00313             if (isBb){
00314                 if (m_thsvc->getHist("/DQAHist/MDC/hresAllIncBb", htmp).isSuccess()) {
00315                     htmp->Fill(resilrInc);
00316                 } else {
00317                     log << MSG::ERROR << "Couldn't retrieve hresAllIncBb" << endreq;
00318                 }
00319                 if (m_thsvc->getHist("/DQAHist/MDC/hresAllExcBb", htmp).isSuccess()) {
00320                     htmp->Fill(resilrExc);
00321                 } else {
00322                     log << MSG::ERROR << "Couldn't retrieve hresAllExcBb" << endreq;
00323                 }     
00324                 if (m_thsvc->getHist("/DQAHist/MDC/hresAllEvaBb", htmp).isSuccess()) {
00325                     htmp->Fill(resilrEva);
00326                 } else {
00327                     log << MSG::ERROR << "Couldn't retrieve hresAllEvaBb" << endreq;
00328                 }     
00329 
00330             } else if (isDimu){
00331                 if (m_thsvc->getHist("/DQAHist/MDC/hresAllIncDimu", htmp).isSuccess()) {
00332                     htmp->Fill(resilrInc);
00333                 } else {
00334                     log << MSG::ERROR << "Couldn't retrieve hresAllIncDimu" << endreq;
00335                 }
00336                 if (m_thsvc->getHist("/DQAHist/MDC/hresAllExcDimu", htmp).isSuccess()) {
00337                     htmp->Fill(resilrExc);
00338                 } else {
00339                     log << MSG::ERROR << "Couldn't retrieve hresAllExcDimu" << endreq;
00340                 }     
00341                 if (m_thsvc->getHist("/DQAHist/MDC/hresAllEvaDimu", htmp).isSuccess()) {
00342                     htmp->Fill(resilrEva);
00343                 } else {
00344                     log << MSG::ERROR << "Couldn't retrieve hresAllEvaDimu" << endreq;
00345                 }     
00346             }
00347         }
00348     }
00349 
00350     return StatusCode::SUCCESS;
00351 }
00352 
00353 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
00354 StatusCode DQA_MDC::finalize() {
00355 
00356     MsgStream log(msgSvc(), name());
00357     log << MSG::INFO << "in finalize()" << endmsg;
00358     log << MSG::INFO << "Finalize successfully! " << endmsg;
00359 
00360     return StatusCode::SUCCESS;
00361 }
00362 
00363
/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/DQA/DQA_MDC/DQA_MDC-00-00-03/src/DQA_MDC.cxx
