BOSS_7.0.2: /home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/Pi0EtaToGGRecAlg/Pi0EtaToGGRecAlg-00-00-11/src/Pi0EtaToGGRecAlg.cxx Source File

00001 #include "GaudiKernel/SmartIF.h"
00002 #include "GaudiKernel/MsgStream.h"
00003 #include "GaudiKernel/SmartDataPtr.h"
00004 #include "GaudiKernel/IDataManagerSvc.h"
00005 
00006 #include "EventModel/EventHeader.h"
00007 #include "EvtRecEvent/EvtRecEvent.h"
00008 #include "EvtRecEvent/EvtRecTrack.h"
00009 #include "EvtRecEvent/EvtRecPi0.h"
00010 #include "EvtRecEvent/EvtRecEtaToGG.h"
00011 #include "VertexFit/KalmanKinematicFit.h"
00012 
00013 #include "Pi0EtaToGGRecAlg/Pi0EtaToGGRecAlg.h"
00014 
00015 const double xmpi0 =  0.134976;  // BOSS 6.4.1 pdt.table value 
00016 const double xmeta =  0.54784;   // PDG08 = 0.547853, BOSS 6.4.1 pdt.table = 0.54784
00017 
00018 //*************************************************************************************
00019 Pi0EtaToGGRecAlg::Pi0EtaToGGRecAlg(const std::string& name, ISvcLocator* pSvcLocator) :
00020    Algorithm(name, pSvcLocator)
00021 {
00022    //Declare the properties
00023 
00024    declareProperty("RejectBothInEndcap",  m_rejectEndEnd = false);
00025 
00027    declareProperty("Pi0MinMass",  m_pi0minmass_cut = 0.08);
00028    declareProperty("Pi0MaxMass",  m_pi0maxmass_cut = 0.18);
00029    declareProperty("Pi0ChisqCut", m_pi0chisq_cut   = 2500);
00030 
00032    declareProperty("EtaMinMass",  m_etaminmass_cut = 0.4);
00033    declareProperty("EtaMaxMass",  m_etamaxmass_cut = 0.7);
00034    declareProperty("EtaChisqCut", m_etachisq_cut   = 2500);
00035 
00036 
00038    declareProperty("PhotonMinEnergy", m_minEnergy = 0.025);
00039 
00040    declareProperty("PhotonInBarrelOrEndcap",  m_useBarrelEndcap   = true);
00041    declareProperty("PhotonMaxCosThetaBarrel", m_maxCosThetaBarrel = 0.8);
00042    declareProperty("PhotonMinCosThetaEndcap", m_minCosThetaEndcap = 0.86);
00043    declareProperty("PhotonMaxCosThetaEndcap", m_maxCosThetaEndcap = 0.92);
00044    declareProperty("PhotonMinEndcapEnergy",   m_minEndcapEnergy   = 0.050);
00045 
00046    declareProperty("PhotonApplyTimeCut", m_applyTimeCut = true);
00047    declareProperty("PhotonMinTime",      m_minTime      = 0.);
00048    declareProperty("PhotonMaxTime",      m_maxTime      = 14.);
00049    declareProperty("PhotonDeltaTime",    m_deltaTime    = 10.);
00050 
00051    declareProperty("PhotonApplyDangCut", m_applyDangCut = false);
00052    declareProperty("PhotonMinDang",      m_minDang      = 20.0);
00053 
00055    declareProperty("MaxNGoodPhoton",     m_maxNGoodPhoton  =  50);
00056    declareProperty("MaxNPi0",            m_maxNPi0         =  100);
00057    //declareProperty("MaxNEta",            m_maxNEta         =  100);
00058 
00059    declareProperty("DoClear",            m_doClear         =  false);
00060 }
00061 
00062 //******************************************************************************************
00063 StatusCode Pi0EtaToGGRecAlg::initialize() {
00064 
00065    MsgStream log(msgSvc(), name());
00066    log << MSG::INFO << "in initialize()" <<endreq;
00067 
00068    return StatusCode::SUCCESS;
00069 }
00070 
00071 
00072 //*********************************************************************************************
00073 StatusCode Pi0EtaToGGRecAlg::execute() {
00074 
00075    MsgStream log(msgSvc(), name());
00076    log << MSG::INFO << "in execute()" <<endreq;
00077 
00078    // get event header, eventlist and tracklist from TDS
00079    SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00080    SmartDataPtr<EvtRecEvent> recEvt(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00081    SmartDataPtr<EvtRecTrackCol> recTrkCol(eventSvc(),  EventModel::EvtRec::EvtRecTrackCol);
00082 
00083    SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), EventModel::EvtRec::EvtRecPi0Col);
00084    if ( !recPi0Col ) {
00085       recPi0Col = new EvtRecPi0Col;
00086       StatusCode sc = eventSvc()->registerObject(EventModel::EvtRec::EvtRecPi0Col, recPi0Col);
00087       if ( sc.isFailure() ) {
00088          log << MSG::ERROR << "could not register EvtRecPi0Col in TDS" <<endreq;
00089          return StatusCode::FAILURE;
00090       }
00091    }
00092 
00093    SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(), EventModel::EvtRec::EvtRecEtaToGGCol);
00094    if ( !recEtaToGGCol ) {
00095       recEtaToGGCol = new EvtRecEtaToGGCol;
00096       StatusCode sc = eventSvc()->registerObject(EventModel::EvtRec::EvtRecEtaToGGCol, recEtaToGGCol);
00097       if ( sc.isFailure() ) {
00098          log << MSG::ERROR << "could not register EvtRecEtaToGGCol in TDS" <<endreq;
00099          return StatusCode::FAILURE;
00100       }
00101    }
00102 
00103    if ( m_doClear ) {
00104        SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
00105        dataManSvc->clearSubTree(EventModel::EvtRec::EvtRecPi0Col);
00106        dataManSvc->clearSubTree(EventModel::EvtRec::EvtRecEtaToGGCol);
00107    }
00108 
00109    // select good photons
00110    std::vector<EvtRecTrack*> vGoodPhotons;
00111    for ( int i = recEvt->totalCharged(); i < recEvt->totalTracks(); i++ ) {
00112       EvtRecTrack* gTrk = *(recTrkCol->begin() + i);
00113       if ( !validPhoton(gTrk,recEvt,recTrkCol) ) continue;
00114       vGoodPhotons.push_back(gTrk);
00115    }
00116 
00117    int nGoodPhoton = vGoodPhotons.size();
00118    if ( nGoodPhoton > m_maxNGoodPhoton ) return StatusCode::SUCCESS;
00119 
00121    KalmanKinematicFit * kmfit = KalmanKinematicFit::instance();
00122 
00123    //
00124    //  Form two-photon combinations
00125    //
00126    for(int i1 = 0; i1 < (nGoodPhoton-1); i1++) {
00127       EvtRecTrack* g1Trk = vGoodPhotons[i1];
00128       RecEmcShower* g1Shower = g1Trk->emcShower();
00129       HepLorentzVector g1P4 = getP4(g1Shower);
00130 
00131       for(int i2 = i1+1; i2 < nGoodPhoton; i2++) {
00132          EvtRecTrack* g2Trk = vGoodPhotons[i2];
00133          RecEmcShower* g2Shower = g2Trk->emcShower();
00134          HepLorentzVector g2P4 = getP4(g2Shower);
00135 
00136          HepLorentzVector p2g = g1P4 + g2P4;
00137 
00139          if(m_rejectEndEnd&&bothInEndcap(g1P4,g2P4)) continue;
00140 
00142          if ((p2g.m() > m_pi0minmass_cut) && (p2g.m() < m_pi0maxmass_cut )) 
00143          {
00144 
00146            kmfit->init();
00147            kmfit->setIterNumber(5);
00148            kmfit->AddTrack(0, 0.0, g1Shower);
00149            kmfit->AddTrack(1, 0.0, g2Shower);
00150            kmfit->AddResonance(0, xmpi0, 0, 1);
00151 
00152            kmfit->Fit(0);  // Perform fit
00153            kmfit->BuildVirtualParticle(0);
00154 
00155            double pi0_chisq = kmfit->chisq(0);
00156            if ( pi0_chisq >= m_pi0chisq_cut ) continue;
00157 
00158 
00160            EvtRecPi0* recPi0 = new EvtRecPi0();
00161            recPi0->setUnconMass(p2g.restMass());
00162            recPi0->setChisq(pi0_chisq);
00163 
00164            if ( g1P4.e() >= g2P4.e() ) {
00165              recPi0->setHiPfit(kmfit->pfit(0));
00166              recPi0->setLoPfit(kmfit->pfit(1));
00167              recPi0->setHiEnGamma(g1Trk);
00168              recPi0->setLoEnGamma(g2Trk);
00169            }
00170            else {
00171              recPi0->setHiPfit(kmfit->pfit(1));
00172              recPi0->setLoPfit(kmfit->pfit(0));
00173              recPi0->setHiEnGamma(g2Trk);
00174              recPi0->setLoEnGamma(g1Trk);
00175            }
00176            recPi0Col->push_back(recPi0);
00177 
00178          }  // End of pi0 mass window IF
00179 
00180 
00182          if ((p2g.m() > m_etaminmass_cut) && (p2g.m() < m_etamaxmass_cut )) 
00183          {
00184 
00186            kmfit->init();
00187            kmfit->setIterNumber(5);
00188            kmfit->AddTrack(0, 0.0, g1Shower);
00189            kmfit->AddTrack(1, 0.0, g2Shower);
00190            kmfit->AddResonance(0, xmeta, 0, 1);
00191 
00192            kmfit->Fit(0);  // Perform fit
00193            kmfit->BuildVirtualParticle(0);
00194 
00195            double eta_chisq = kmfit->chisq(0);
00196            if ( eta_chisq >= m_etachisq_cut ) continue;
00197 
00198 
00200            EvtRecEtaToGG* recEtaToGG = new EvtRecEtaToGG();
00201            recEtaToGG->setUnconMass(p2g.restMass());
00202            recEtaToGG->setChisq(eta_chisq);
00203 
00204            if ( g1P4.e() >= g2P4.e() ) {
00205              recEtaToGG->setHiPfit(kmfit->pfit(0));
00206              recEtaToGG->setLoPfit(kmfit->pfit(1));
00207              recEtaToGG->setHiEnGamma(g1Trk);
00208              recEtaToGG->setLoEnGamma(g2Trk);
00209            }
00210            else {
00211              recEtaToGG->setHiPfit(kmfit->pfit(1));
00212              recEtaToGG->setLoPfit(kmfit->pfit(0));
00213              recEtaToGG->setHiEnGamma(g2Trk);
00214              recEtaToGG->setLoEnGamma(g1Trk);
00215            }
00216            recEtaToGGCol->push_back(recEtaToGG);
00217 
00218          }  // End of etaToGG mass window IF
00219 
00220       }  // End of g2Trk evtRec FOR 
00221    }  // End of g1Trk evtRec FOR 
00222 
00223    // If there are too many pi0s, it's possibly a bad event and all pi0s are abandoned
00224    if ( recPi0Col->size() > m_maxNPi0 ) {
00225       recPi0Col->clear();
00226    }
00227 
00228    return StatusCode::SUCCESS;
00229 }
00230 
00231 
00232 //********************************************************************************************
00233 StatusCode Pi0EtaToGGRecAlg::finalize() {
00234 
00235    MsgStream log(msgSvc(), name());
00236    log << MSG::INFO << "in finalize()" << endreq;
00237 
00238    return StatusCode::SUCCESS;
00239 }
00240 
00241 
00242 
00244 HepLorentzVector Pi0EtaToGGRecAlg::getP4(RecEmcShower* gTrk){
00245 
00246    double eraw = gTrk->energy();
00247    double phi =  gTrk->phi();
00248    double the =  gTrk->theta();
00249 
00250    return HepLorentzVector( eraw * sin(the) * cos(phi),
00251                             eraw * sin(the) * sin(phi),
00252                             eraw * cos(the),
00253                             eraw );
00254 }
00255 
00256 
00257 bool Pi0EtaToGGRecAlg::validPhoton(EvtRecTrack* recTrk, 
00258                                    SmartDataPtr<EvtRecEvent> recEvt, 
00259                                    SmartDataPtr<EvtRecTrackCol> recTrkCol) 
00260 {
00261    if ( !recTrk->isEmcShowerValid() ) return false;
00262 
00263    RecEmcShower *emcTrk = recTrk->emcShower();
00264 
00265    HepLorentzVector shP4 = getP4(emcTrk);
00266    double cosThetaSh = shP4.vect().cosTheta();
00267 
00268 
00270    if (shP4.e() <= m_minEnergy) return false;
00271 
00272 
00274    if ( m_useBarrelEndcap ) {
00275       bool inBarrelEndcap = false;
00276 
00277       if(fabs(cosThetaSh) < m_maxCosThetaBarrel) inBarrelEndcap = true;
00278 
00279       if((fabs(cosThetaSh) > m_minCosThetaEndcap)
00280             &&(fabs(cosThetaSh) < m_maxCosThetaEndcap)
00281             &&(shP4.e() > m_minEndcapEnergy)) inBarrelEndcap = true;
00282 
00283       if ( !inBarrelEndcap ) return false;
00284    }
00285 
00286 
00288    if ( m_applyTimeCut ) {
00289       double time = emcTrk->time();
00290       if ( recEvt->totalCharged() > 0 ) {
00291          if ( time < m_minTime || time > m_maxTime ) return false;
00292       }
00293       else {
00294          RecEmcShower* firstG = (*(recTrkCol->begin()))->emcShower();
00295          double deltaTime = fabs(time - firstG->time());
00296          if ( deltaTime > 10 ) return false;
00297       }
00298    }
00299 
00300 
00302    if (m_applyDangCut && recEvt->totalCharged() > 0) 
00303    {      
00304       Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
00305       double dang = 200.; 
00306 
00307       for (int j = 0; j < recEvt->totalCharged(); j++) {
00308          EvtRecTrackIterator jtTrk = recTrkCol->begin() + j;
00309          if ( !(*jtTrk)->isExtTrackValid() ) continue;
00310          RecExtTrack* extTrk = (*jtTrk)->extTrack();
00311          if ( extTrk->emcVolumeNumber() == -1 ) continue;
00312          Hep3Vector extpos = extTrk->emcPosition();
00313          double  angd1 = extpos.angle(emcpos);
00314          if ( angd1 < dang ) dang = angd1;
00315       }
00316 
00317       if ( dang < 200 ) {
00318          dang = dang * 180 / (CLHEP::pi);
00319          if ( dang <= m_minDang ) return false;
00320       }
00321    }  // End of "recEvt->totalCharged() > 0" IF
00322 
00323    return true;
00324 }
00325 
00326 
00327 bool Pi0EtaToGGRecAlg::bothInEndcap(HepLorentzVector g1_P4, HepLorentzVector g2_P4){
00328 
00329   double g1_CosTheta = g1_P4.vect().cosTheta();
00330   double g2_CosTheta = g2_P4.vect().cosTheta();
00331 
00332   bool g1InEndcap = false;
00333   bool g2InEndcap = false;
00334 
00335   if((fabs(g1_CosTheta) > m_minCosThetaEndcap)
00336      &&(fabs(g1_CosTheta) < m_maxCosThetaEndcap)) g1InEndcap = true;
00337 
00338   if((fabs(g2_CosTheta) > m_minCosThetaEndcap)
00339      &&(fabs(g2_CosTheta) < m_maxCosThetaEndcap)) g2InEndcap = true;
00340 
00341   if(g1InEndcap&&g2InEndcap) return( true );
00342 
00343   return( false );
00344 }