EventAssemblyAlg Class Reference

#include <EventAssemblyAlg.h>

List of all members.

Public Member Functions

 EventAssemblyAlg (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode initialize ()
StatusCode execute ()
StatusCode finalize ()

Private Member Functions

double thetaCut (double p, double nSigma, int parId)
double phiCut (double p, double nSigma, int parId)

Private Attributes

bool myMsgFlag
bool myActiveCutFlag
double myEmcThetaCut
double myEmcPhiCut
double myEmcThetaNSigmaCut
double myEmcPhiNSigmaCut
int myParticleId
int m_Output
NTuple::Tuple * myNTuple1
NTuple::Item< double > myExted
NTuple::Item< double > myMatched
NTuple::Item< double > myEnergyMatched


Detailed Description

Definition at line 19 of file EventAssemblyAlg.h.


Constructor & Destructor Documentation

EventAssemblyAlg::EventAssemblyAlg ( const std::string name,
ISvcLocator *  pSvcLocator 
)

Definition at line 23 of file EventAssemblyAlg.cxx.

References m_Output, myActiveCutFlag, myEmcPhiCut, myEmcPhiNSigmaCut, myEmcThetaCut, myEmcThetaNSigmaCut, myMsgFlag, and myParticleId.

00023                                                                                  :Algorithm(name, pSvcLocator)
00024 {
00025   //Declare the properties
00026   declareProperty("MsgFlag",myMsgFlag=false);
00027   declareProperty("ActiveCutFlag",myActiveCutFlag=true);
00028   declareProperty("EmcThetaMatchCut",myEmcThetaCut=-1);
00029   declareProperty("EmcPhiMatchCut",myEmcPhiCut=-1);
00030   declareProperty("EmcThetaNSigmaCut",myEmcThetaNSigmaCut=5);
00031   declareProperty("EmcPhiNSigmaCut",myEmcPhiNSigmaCut=5);
00032   declareProperty("ParticleId",myParticleId);
00033   declareProperty("Output", m_Output = 0);
00034 }


Member Function Documentation

StatusCode EventAssemblyAlg::execute (  ) 

Definition at line 68 of file EventAssemblyAlg.cxx.

References abs, TofID::barrel_ec(), EventModel::EvtRec::Event, EventModel::EvtRec::EvtRecEvent, EventModel::EvtRec::EvtRecTrackCol, Bes_Common::FATAL, genRecEmupikp::i, Bes_Common::INFO, TofID::is_mrpc(), ganga-rec::j, TofID::layer(), m_Output, M_PI, TofID::module(), msgSvc(), myActiveCutFlag, myEmcPhiCut, myEmcPhiNSigmaCut, myEmcThetaCut, myEmcThetaNSigmaCut, myEnergyMatched, myExted, myMatched, myMsgFlag, myNTuple1, myParticleId, TofID::phi_module(), phiCut(), EventModel::Recon::RecEmcShowerCol, EventModel::Recon::RecExtTrackCol, EventModel::Recon::RecMdcDedxCol, EventModel::Recon::RecMdcKalTrackCol, EventModel::Recon::RecMdcTrackCol, EventModel::Recon::RecMucTrackCol, EventModel::Recon::RecTofTrackCol, EvtRecEvent::setTotalCharged(), EvtRecEvent::setTotalNeutral(), EvtRecEvent::setTotalTracks(), thetaCut(), and x.

00069 {
00070   MsgStream log(msgSvc(), name());
00071   log << MSG::INFO << "in execute()" << endreq;
00072 
00073   int numOfChargedTrks = 0;
00074   int numOfKalTrks = 0;
00075   int numOfDedx = 0;
00076   int numOfExt = 0;
00077   int numOfTof = 0;
00078   int numOfShower = 0;
00079   int numOfMatchedShower = 0;
00080   int numOfMucTrks = 0;
00081         
00082   // Part 1: Get the event header, print out event and run number
00083   SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00084   if (!eventHeader)
00085     {
00086       log << MSG::FATAL << "Could not find Event Header" << endreq;
00087       return( StatusCode::FAILURE);
00088     }
00089   log << MSG::INFO << ": retrieved event: " << eventHeader->eventNumber()  << "  run: " <<  eventHeader->runNumber() << endreq;
00090 
00091   //Part 2: Retrieve Rec data
00092   SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00093   RecMdcTrackCol::iterator beginOfMdcTrkCol = RecMdcTrackCol::iterator(NULL);
00094   if(!recMdcTrackCol)
00095     {
00096       log<<MSG::INFO<<"Could not find RecMdcTrackCol!"<<endreq;
00097     }
00098   else 
00099     {
00100       beginOfMdcTrkCol = recMdcTrackCol->begin();
00101       numOfChargedTrks = recMdcTrackCol->size();
00102     }
00103   
00104   SmartDataPtr<RecMdcKalTrackCol> mdcKalTrackCol(eventSvc(), EventModel::Recon::RecMdcKalTrackCol);
00105   RecMdcKalTrackCol::iterator beginOfMdcKalTrackCol = RecMdcKalTrackCol::iterator(NULL);
00106   if(!mdcKalTrackCol)
00107     {
00108       log<<MSG::INFO<<"Could not find MdcKalTrackCol!"<<endreq;
00109     }
00110   else 
00111     {
00112       numOfKalTrks = mdcKalTrackCol->size();
00113       beginOfMdcKalTrackCol = mdcKalTrackCol->begin();
00114     }
00115                 
00116   
00117   SmartDataPtr<RecMdcDedxCol> mdcDedxCol(eventSvc(), EventModel::Recon::RecMdcDedxCol);
00118   RecMdcDedxCol::iterator beginOfDedxCol = RecMdcDedxCol::iterator(NULL);
00119   if(!mdcDedxCol)
00120     {
00121       log<<MSG::INFO<<"Could not find MdcDedxCol!"<<endreq;
00122     }
00123     else {
00124       numOfDedx = mdcDedxCol->size();
00125       beginOfDedxCol = mdcDedxCol->begin();
00126     }
00127   
00128 
00129   SmartDataPtr<RecExtTrackCol> extTrackCol(eventSvc(), EventModel::Recon::RecExtTrackCol);
00130   RecExtTrackCol::iterator beginOfExtTrackCol = RecExtTrackCol::iterator(NULL);
00131   if(!extTrackCol)
00132     {
00133       log<<MSG::INFO<<"Could not find RecExtTrackCol!"<<endreq;
00134     }
00135   else
00136     {
00137       beginOfExtTrackCol = extTrackCol->begin();
00138       numOfExt = extTrackCol->size();
00139     }
00140   
00141   SmartDataPtr<RecTofTrackCol> TofTrackCol(eventSvc(), EventModel::Recon::RecTofTrackCol);
00142   RecTofTrackCol::iterator beginOfTofTrackCol = RecTofTrackCol::iterator(NULL);
00143   if(!TofTrackCol)
00144     {
00145       log<<MSG::INFO<<"Could not find RecTofTrackCol!"<<endreq;
00146     }
00147   else
00148     {
00149       beginOfTofTrackCol = TofTrackCol->begin();
00150       numOfTof = TofTrackCol->size();
00151     }
00152         
00153   SmartDataPtr<RecEmcShowerCol> emcRecShowerCol(eventSvc(), EventModel::Recon::RecEmcShowerCol);
00154   RecEmcShowerCol::iterator beginOfEmcTrackCol = RecEmcShowerCol::iterator(NULL);
00155   if(!emcRecShowerCol)
00156     {
00157       log<<MSG::INFO<<"Could not find RecEmcShowerCol!"<<endreq;
00158     }
00159   else
00160     {
00161       beginOfEmcTrackCol = emcRecShowerCol->begin();
00162       numOfShower = emcRecShowerCol->size();
00163     }
00164         
00165   SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc(), EventModel::Recon::RecMucTrackCol);
00166   RecMucTrackCol::iterator beginOfMucTrackCol = RecMucTrackCol::iterator(NULL);
00167   if(!mucTrackCol)
00168     {
00169       log<<MSG::INFO<<"Could not find RecMucTrackCol!"<<endreq;
00170     }
00171   else
00172     {
00173       beginOfMucTrackCol = mucTrackCol->begin();
00174       numOfMucTrks = mucTrackCol->size();
00175     }
00176 
00177     //Some match flag
00178     
00179     int MaxHit = numOfChargedTrks+numOfShower+numOfTof;
00180 
00181     vector<bool> dEdxMatched;
00182     vector<bool> kalTrkMatched;
00183     vector<bool> extMatched;
00184     vector<bool> TofMatched;
00185     vector<bool> emcMatched;
00186     vector<bool> mucMatched;
00187 
00188     for(int i=0;i<MaxHit;i++)
00189     {
00190       dEdxMatched.push_back(false);
00191       kalTrkMatched.push_back(false);
00192       extMatched.push_back(false);
00193       TofMatched.push_back(false);
00194       emcMatched.push_back(false);
00195       mucMatched.push_back(false);
00196     }
00197 
00198 
00199   //Part 3: Match tracks
00200   //  RecTrackListCol* aNewRecTrkListCol = new RecTrackListCol;
00201   EvtRecTrackCol* aNewEvtRecTrackCol = new EvtRecTrackCol;
00202   if(numOfChargedTrks+numOfShower)//if mdc trackCol and emc trackCol neither exists,there is no match job.
00203     {
00204       for(int i=0;i<numOfChargedTrks;i++)//match charged tracks
00205         {
00206           EvtRecTrack* aNewEvtRecTrack = new EvtRecTrack;
00207           int mdcTrkID = (*(beginOfMdcTrkCol+i))->trackId();
00208           aNewEvtRecTrack->setTrackId(mdcTrkID);
00209           aNewEvtRecTrack->setMdcTrack(*(beginOfMdcTrkCol+i));
00210           //      aNewEvtRecTrack->setMdcTrkIdx(i);
00211           for(int j=0;j<numOfKalTrks;j++)//match MdcKalTrack
00212             if( !kalTrkMatched[j] && mdcTrkID==(*(beginOfMdcKalTrackCol+j))->trackId() )
00213               {
00214                 aNewEvtRecTrack->setMdcKalTrack(*(beginOfMdcKalTrackCol+j));
00215                 //              aNewRecTrkList->setMdcKalTrkIdx(j);
00216                 kalTrkMatched[j]=true;
00217               }
00218           
00219           for(int j=0;j<numOfDedx;j++)//match dEdx
00220             if( !dEdxMatched[j] && mdcTrkID==(*(beginOfDedxCol+j))->trackId() )
00221               {
00222                 aNewEvtRecTrack->setMdcDedx(*(beginOfDedxCol+j));
00223                 //              aNewRecTrkList->setDedxTrkIdx(j);
00224                 dEdxMatched[j]=true;
00225               }
00226             
00227           for(int j=0;j<numOfExt;j++)//match ExtTrack
00228             if( !extMatched[j] && mdcTrkID==(*(beginOfExtTrackCol+j))->trackId() )
00229               {
00230                 aNewEvtRecTrack->setExtTrack(*(beginOfExtTrackCol+j));
00231                 //              aNewRecTrkList->setExtTrkIdx(j);
00232                 extMatched[j]=true;
00233               }
00234           for(int j=0;j<numOfTof;j++)//match Tof
00235             if( !TofMatched[j] && mdcTrkID==(*(beginOfTofTrackCol+j))->trackID() )
00236               {
00237                 aNewEvtRecTrack->addTofTrack(*(beginOfTofTrackCol+j));
00238                 //              aNewRecTrkList->setTofTrkIdx(j);
00239                 TofMatched[j]=true;
00240               }
00241           for(int j=0;j<numOfMucTrks;j++)//match MucTrack
00242             if( !mucMatched[j] && mdcTrkID==(*(beginOfMucTrackCol+j))->trackId() )
00243               {
00244                 aNewEvtRecTrack->setMucTrack(*(beginOfMucTrackCol+j));
00245                 //              aNewRecTrkList->setMucTrkIdx(j);
00246                 mucMatched[j]=true;
00247               }
00248 
00249           if(m_Output==1){
00250             myExted = -1.;
00251             myMatched = -1.;
00252             myEnergyMatched = -1.;
00253           }
00254           if(aNewEvtRecTrack->isExtTrackValid()&&(aNewEvtRecTrack->extTrack())->emcVolumeNumber()!=-1)//start to match Emc showers
00255             {
00256               if(m_Output==1)myExted = 1.;
00257               //resieve ext information @ Emc
00258               Hep3Vector extPos = (aNewEvtRecTrack->extTrack())->emcPosition();
00259               double extTheta = extPos.theta();
00260               double extPhi = extPos.phi();
00261               
00262               double dAngleMin = 9999.;
00263               int indexOfEmcNearest = -1;
00264               
00265               for(int j=0;j<numOfShower;j++)//start to find nearest emc shower
00266                 {
00267                   if(!emcMatched[j])
00268                     {
00269                       HepPoint3D emcPot = (*(beginOfEmcTrackCol+j))->position();
00270                       Hep3Vector emcPos(emcPot.x(),emcPot.y(),emcPot.z());
00271                       double dAngle = extPos.angle(emcPos);
00272                       if(dAngle<dAngleMin)
00273                         {
00274                           dAngleMin = dAngle;
00275                           indexOfEmcNearest = j;
00276                         }
00277                     }
00278                 }
00279               if(indexOfEmcNearest>=0)//if nearest Emc shower is found
00280                 {
00281                   HepPoint3D emcPot = (*(beginOfEmcTrackCol+indexOfEmcNearest))->position();
00282                   double theta = emcPot.theta();
00283                   double phi = emcPot.phi();
00284                   double deltaTheta = (extTheta-theta);
00285                   double deltaPhi = fmod(extPhi-phi+5*M_PI,2*M_PI)-M_PI;
00286                   if(myActiveCutFlag)
00287                     {
00288                       double tanLamda = (*(beginOfMdcTrkCol+i))->helix(4);
00289                       double kappa = (*(beginOfMdcTrkCol+i))->helix(2);
00290                       double p = sqrt(1+tanLamda*tanLamda)/fabs(kappa);
00291                       double aThetaCut = thetaCut(p,myEmcThetaNSigmaCut,myParticleId);
00292                       double aPhiCUt = phiCut(p,myEmcPhiNSigmaCut,myParticleId);
00293                       if(fabs(deltaTheta)<aThetaCut  && fabs(deltaPhi)<aPhiCUt)
00294                         {
00295                           aNewEvtRecTrack->setEmcShower(*(beginOfEmcTrackCol+indexOfEmcNearest));
00296                           //                      aNewRecTrkList->setEmcTrkIdx(indexOfEmcNearest);
00297                           emcMatched[indexOfEmcNearest]=true;
00298                           numOfMatchedShower++;
00299                           if(m_Output==1)myMatched = 1.;
00300                         }
00301                     }
00302                   else if(fabs(deltaTheta)<myEmcThetaCut && fabs(deltaPhi)<myEmcPhiCut)//if the nearest shower is in the fixed match window,it will be matched.
00303                     {
00304                       aNewEvtRecTrack->setEmcShower(*(beginOfEmcTrackCol+indexOfEmcNearest));
00305                       //                      aNewRecTrkList->setEmcTrkIdx(indexOfEmcNearest);
00306                       emcMatched[indexOfEmcNearest]=true;
00307                       numOfMatchedShower++;
00308                       if(m_Output==1)myMatched = 1.;
00309                     }
00310                 }
00311             }                                   
00312           aNewEvtRecTrackCol->push_back(aNewEvtRecTrack);
00313           if(m_Output==1){
00314             myEnergyMatched = 0.;
00315             myNTuple1->write();
00316           }
00317         }//end charged track match
00318                 
00319       int id=0;
00320       for(int i=0;i<numOfShower;i++)//If Emc shower isn't macthed, it will be added to the end of RecTrackListCol
00321         {
00322           if(!emcMatched[i])
00323             {
00324               EvtRecTrack* aNewEvtRecTrack = new EvtRecTrack;
00325               aNewEvtRecTrack->setTrackId(id+numOfChargedTrks);
00326               aNewEvtRecTrack->setEmcShower(*(beginOfEmcTrackCol+i));
00327               //              aNewRecTrkList->setEmcTrkIdx(i);
00328 
00329   //EMC and TOF match
00330         HepPoint3D emcPos((*(beginOfEmcTrackCol+i))->position());
00331         double emcPhi = emcPos.phi();
00332         emcPhi = emcPhi < 0 ? emcPhi+CLHEP::twopi : emcPhi;
00333         
00334         int module=9999, nphi1=9999, nphi2=9999, nphi1_mrpc=9999;
00335         module = (*(beginOfEmcTrackCol+i))->module();
00336         if(module==1) {    //barrel
00337           nphi1 = (int)(emcPhi*88./CLHEP::twopi);      //tof layer1 number
00338           nphi2 = (int)(emcPhi*88./CLHEP::twopi+0.5);  //tof layer2 number
00339           nphi2 += 88;
00340         } else {    //endcap
00341           nphi1 = (int)(emcPhi*48./CLHEP::twopi+0.5);
00342           
00343 
00344       //new mrpc detector 
00345           double fi_endtof_mrpc = atan2(emcPos.y(),emcPos.x())-3.141592654/2.; 
00346           if(fi_endtof_mrpc<0) fi_endtof_mrpc += 2*3.141592654;
00347           nphi1_mrpc=((int)(fi_endtof_mrpc/(3.141593/18)))+1; //We have 36 modules in each endcap (divided again into two layers)! The numbering starts with 1.
00348 
00349           if(nphi1_mrpc%2==0){nphi1_mrpc=nphi1_mrpc/2;}  //Calculate the correct number as in each endcap we have two layers of modules furnished with numbers 1 to 18.
00350           else {nphi1_mrpc=(nphi1_mrpc+1)/2;}
00351 
00352           if(emcPos.z()>0) {(nphi1_mrpc -= 19); nphi1_mrpc=nphi1_mrpc*(-1);} //Consider that for the east endcap the module at x=0;y=max has the number 18 (the numbering is != to west endcap)
00353           // nphi1_mrpc is now the number of the nearest ToF module!
00354          
00355          }
00356 
00357         int dPhiMin = 9999;
00358         int partid_dPhiMin = 7;
00359         bool mrpcmatch = false;
00360 
00361         int indexOfTofNearest = -1;
00362         for(int j=0;j<numOfTof;j++) { //start to find nearest Tof shower
00363           if( !TofMatched[j] ) {
00364 
00365             Identifier id((*(beginOfTofTrackCol+j))->tofID());
00366             bool is_mrpc = TofID::is_mrpc(id);
00367 
00368             if(!is_mrpc)
00369               {
00370                 int barrel_ec = TofID::barrel_ec(id);
00371                 int layer = TofID::layer(id);
00372                 int im = TofID::phi_module(id);
00373                 im += layer * 88;
00374 
00375                 if(module!=barrel_ec) continue;
00376             
00377                 int dPhi;
00378                 if(layer == 0) {
00379                   dPhi = abs(im - nphi1);
00380                 } else {
00381                   dPhi = abs(im - nphi2);
00382                 }
00383                 if(module==1) {
00384                   dPhi = dPhi>=44 ? 88-dPhi : dPhi;
00385                 } else {
00386                   dPhi = dPhi>=24 ? 48-dPhi : dPhi;
00387                 } 
00388 
00389                 if(dPhi < dPhiMin) {
00390                   dPhiMin = dPhi;
00391                   indexOfTofNearest = j;
00392                 }
00393               }//close if(!is_mrpc)
00394             else//is_mrpc
00395               {
00396                 int barrel_ec = TofID::barrel_ec(id);
00397                 int module_mrpc = TofID::module(id);
00398                 
00399                 int dphi;
00400                 dphi = abs(module_mrpc-nphi1_mrpc);
00401                 
00402                 if(dphi < dPhiMin)
00403                   {
00404                     dPhiMin=dphi;
00405                     partid_dPhiMin=barrel_ec;
00406                     indexOfTofNearest = j;
00407                     mrpcmatch = true;
00408                   }
00409               }//close else is_mrpc
00410 
00411           }//close if( !TofMatched[j] ) 
00412         }//close for(int j=0;j<numOfTof;j++)
00413 
00414         /*
00415         if(0)
00416           {
00417         cout << "Event Assembly Algorithm" <<endl;
00418         cout << "EMC Shower number   =  " << i << endl;
00419         cout << "mrpcmatch ?            " << mrpcmatch << endl;
00420         cout << "indexOfTofNearest      " << indexOfTofNearest <<endl;
00421         cout << "dphiMin                " << dPhiMin << endl;
00422         cout << "partid_mrpc            " << partid_dPhiMin <<endl;
00423         cout << "partid_emc             " << module  <<endl <<endl;
00424           }
00425         */
00426 
00427         //match in z-direction
00428         if(indexOfTofNearest>=0) {
00429           HepPoint3D tofPos(0,0.867,(*(beginOfTofTrackCol+indexOfTofNearest))->zrhit());
00430           double matWindow = 0;                    
00431           double eShower = (*(beginOfEmcTrackCol+i))->e5x5();//energy of a 5x5 cluster
00432           if(eShower>0) {
00433             matWindow = 0.01277+0.004268/sqrt(eShower);
00434           }
00435         }
00436 
00437 
00438         if(indexOfTofNearest>=0 && dPhiMin<=1) {  //if nearest Tof shower is found
00439           // calculate matching window for barrel
00440           double eShower = (*(beginOfEmcTrackCol+i))->e5x5();
00441           double matWindow = 0;
00442           if(eShower>0) {
00443             matWindow = 0.01277+0.004268/sqrt(eShower);
00444           }
00445           matWindow = 0.2;  //temp code
00446 
00447           HepPoint3D tofPos(0,0.867,(*(beginOfTofTrackCol+indexOfTofNearest))->zrhit());
00448           if(fabs(tofPos.theta()-emcPos.theta())<matWindow || module!=1) {
00449             aNewEvtRecTrack->addTofTrack(*(beginOfTofTrackCol+indexOfTofNearest));
00450             TofMatched[indexOfTofNearest]=true;
00451             //cout << "EventAssembly :  New neutral Toftrack added  " << endl;
00452 
00453           }
00454         }//close if(indexOfTofNearest>=0 && dPhiMin<=1)
00455 
00456               aNewEvtRecTrackCol->push_back(aNewEvtRecTrack);
00457               id++;
00458             }//close if(!)
00459         }//close for(int i=0;i<numOfShower;i++) <---This are the EMC showers
00460     }//end match
00461 
00462   //match statistics
00463   if(myMsgFlag) cout<<"Charged tracks: "<<numOfChargedTrks<<"   Ext track: "<<numOfExt  <<"  Shower:"<<numOfShower<<"  Matched shower:"<<numOfMatchedShower<<endl;
00464 
00465   //Part 4: register RecTrackListCol and EventList in TDS
00466   DataObject *aDataObject1;
00467   eventSvc()->findObject(EventModel::EvtRec::Event, aDataObject1);
00468   if ( aDataObject1 == NULL ) {
00469      EvtRecObject* aRecObject = new EvtRecObject; 
00470      StatusCode sc = eventSvc()->registerObject(EventModel::EvtRec::Event, aRecObject);
00471      if ( sc.isFailure() ) {
00472         log << MSG::FATAL << "Could not register EvtRecObject in TDS!" << endreq;
00473         return StatusCode::FAILURE;
00474      }
00475   }
00476 
00477   DataObject* aDataObject2;
00478   eventSvc()->findObject(EventModel::EvtRec::EvtRecEvent, aDataObject2);
00479   if ( aDataObject2 == NULL ) {
00480      EvtRecEvent* aRecEvent = new EvtRecEvent;
00481      aRecEvent->setTotalTracks(numOfChargedTrks+numOfShower-numOfMatchedShower);
00482      aRecEvent->setTotalCharged(numOfChargedTrks);
00483      aRecEvent->setTotalNeutral(numOfShower-numOfMatchedShower);
00484      StatusCode sc = eventSvc()->registerObject(EventModel::EvtRec::EvtRecEvent, aRecEvent);
00485      if ( sc.isFailure() ) {
00486         log << MSG::FATAL << "Could not register EvtRecEvent in TDS!" << endreq;
00487         return( StatusCode::FAILURE);
00488      }
00489   }
00490   else {
00491      EvtRecEvent* aRecEvent = dynamic_cast<EvtRecEvent*>(aDataObject2);
00492      aRecEvent->setTotalTracks(numOfChargedTrks+numOfShower-numOfMatchedShower);
00493      aRecEvent->setTotalCharged(numOfChargedTrks);
00494      aRecEvent->setTotalNeutral(numOfShower-numOfMatchedShower);
00495   }
00496 
00497   DataObject* aDataObject3;
00498   eventSvc()->findObject(EventModel::EvtRec::EvtRecTrackCol, aDataObject3);
00499   if ( aDataObject3 != NULL ) {
00500     //IDataManagerSvc* dataMgrSvc = dynamic_cast<IDataManagerSvc*> (eventSvc());
00501     SmartIF<IDataManagerSvc> dataMgrSvc(eventSvc());
00502      dataMgrSvc->clearSubTree(EventModel::EvtRec::EvtRecTrackCol);
00503      eventSvc()->unregisterObject(EventModel::EvtRec::EvtRecTrackCol);
00504   }
00505   
00506   StatusCode sc = eventSvc()->registerObject(EventModel::EvtRec::EvtRecTrackCol, aNewEvtRecTrackCol);
00507   if ( sc.isFailure() ) {
00508      log << MSG::FATAL << "Could not register EvtRecTrackCol in TDS!" << endreq;
00509      return( StatusCode::FAILURE);
00510   }
00511         
00512   //Part 5: check RecTrackListCol in TDS
00513   SmartDataPtr<EvtRecTrackCol> evtRecTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00514   if(!evtRecTrackCol)
00515     {
00516       log<<MSG::FATAL<<"Could not find RecTrackListCol in TDS!"<<endreq;
00517     }
00518   
00519   EvtRecTrackCol::iterator itOfRecTrkListCol = evtRecTrackCol->begin();
00520   if(myMsgFlag)
00521     for(int i=0;itOfRecTrkListCol!=evtRecTrackCol->end();itOfRecTrkListCol++,i++)
00522       {
00523         cout<<"******************** Track "<<i<<": *********************"<<endl;
00524         cout<<"Track ID: "<<(*itOfRecTrkListCol)->trackId()<<endl;
00525         if((*itOfRecTrkListCol)->isMdcTrackValid())
00526           {
00527             RecMdcTrack* mdcTrk = (*itOfRecTrkListCol)->mdcTrack();
00528             double kappa = mdcTrk->helix(2);
00529             double tanl = mdcTrk->helix(4);
00530             cout<<"Mdc kappa, tanl: "<<kappa<<", "<<tanl<<endl;
00531           }
00532         if((*itOfRecTrkListCol)->isExtTrackValid())
00533           {
00534             RecExtTrack* extTrack = (*itOfRecTrkListCol)->extTrack();
00535             Hep3Vector emcPos = extTrack->emcPosition();
00536             cout<<"Ext Emc pos:"<<emcPos.x()<<","<<emcPos.y()<<","<<emcPos.z()<<endl;
00537           }
00538         if((*itOfRecTrkListCol)->isEmcShowerValid())
00539           {
00540             RecEmcShower* emcTrack = (*itOfRecTrkListCol)->emcShower();
00541             HepPoint3D emcPos = emcTrack->position();
00542             double x = emcPos.x();
00543             double y = emcPos.y();
00544             double z = emcPos.z();
00545             cout<<"Emc rec pos:"<<x<<","<<y<<","<<z<<endl;
00546           }
00547       }
00548   
00549   return StatusCode::SUCCESS;
00550 }

StatusCode EventAssemblyAlg::finalize (  ) 

Definition at line 553 of file EventAssemblyAlg.cxx.

References Bes_Common::INFO, and msgSvc().

00554 {
00555   MsgStream log(msgSvc(), name());
00556   log << MSG::INFO << "in finalize()" << endreq;
00557   
00558   return StatusCode::SUCCESS;
00559 }

StatusCode EventAssemblyAlg::initialize (  ) 

Definition at line 37 of file EventAssemblyAlg.cxx.

References calibUtil::ERROR, Bes_Common::INFO, m_Output, msgSvc(), myEmcPhiCut, myEmcPhiNSigmaCut, myEmcThetaCut, myEmcThetaNSigmaCut, myEnergyMatched, myExted, myMatched, myNTuple1, and ntupleSvc().

00038 {
00039   MsgStream log(msgSvc(), name());
00040   log << MSG::INFO << "in initialize()" << endreq;
00041   
00042   if(myEmcThetaCut<0) myEmcThetaCut = 0.1745;
00043   if(myEmcPhiCut<0) myEmcPhiCut = 0.2618;
00044   
00045   log << MSG::INFO <<"EmcThetaCut,EmcPhiCut: "<<myEmcThetaCut<<", "<<myEmcPhiCut<<endreq;
00046   log << MSG::INFO <<"emcThetaNSigmaCut, emcPhiNSigmaCut: "<<myEmcThetaNSigmaCut<<", "<<myEmcPhiNSigmaCut<<endreq;
00047   
00048   if(m_Output == 1){
00049     NTuplePtr nt1(ntupleSvc(), "FILE998/match");
00050     if ( nt1 ) myNTuple1 = nt1;
00051     else {
00052       myNTuple1 = ntupleSvc()->book("FILE998/match", CLID_ColumnWiseTuple, "match");
00053       if ( myNTuple1 )    {
00054         myNTuple1->addItem("exted",myExted);
00055         myNTuple1->addItem("matched",myMatched);
00056         myNTuple1->addItem("Ematched",myEnergyMatched);
00057       }
00058       else    {
00059         log << MSG::ERROR << " Cannot book N-tuple:" << long(myNTuple1) << endmsg;
00060         //return StatusCode::FAILURE;
00061       }
00062     }
00063   }
00064         return StatusCode::SUCCESS;
00065 }

double EventAssemblyAlg::phiCut ( double  p,
double  nSigma,
int  parId 
) [private]

Definition at line 581 of file EventAssemblyAlg.cxx.

References cut, exp(), and myEmcPhiCut.

Referenced by execute().

00582 {
00583   if(p<0.3) p=0.3;
00584   if(parId<1||parId>3) parId=3;
00585   if(nSigma<=0) nSigma = 5.;
00586   double cut = myEmcPhiCut;
00587   switch(parId)
00588     {
00589     case 1:cut = fabs(exp(-2.187*p-2.945)+0.03236)+(exp(-2.977*p-3.558)+0.005566)*nSigma; break;
00590     case 2:cut = fabs(exp(-2.216*p-2.13)+0.03314)+(exp(-2.436*p-3.392)+0.005049)*nSigma; break;
00591     case 3:cut = fabs(exp(-1.63*p-2.931)+0.03223)+(exp(-3.192*p-2.756)+0.01533)*nSigma; break;
00592     }
00593   
00594   return cut;
00595 }

double EventAssemblyAlg::thetaCut ( double  p,
double  nSigma,
int  parId 
) [private]

Definition at line 564 of file EventAssemblyAlg.cxx.

References cut, exp(), and myEmcThetaCut.

Referenced by execute().

00565 {
00566   if(p<0.3) p=0.3;
00567   if(parId<1||parId>3) parId=3;
00568   if(nSigma<=0) nSigma = 5.;
00569   double cut = myEmcThetaCut;
00570   switch(parId)
00571     {
00572     case 1:cut = fabs(-0.00159)+(exp(-3.566*p-4.14)+0.006804)*nSigma; break;
00573     case 2:cut = fabs(-0.00145)+(exp(-4.268*p-3.305)+0.01129)*nSigma; break;
00574     case 3:cut = fabs(-0.00161)+(exp(-5.955*p-3.052)+0.01751)*nSigma; break;
00575     }
00576 
00577   return cut;
00578 }


Member Data Documentation

int EventAssemblyAlg::m_Output [private]

Definition at line 40 of file EventAssemblyAlg.h.

Referenced by EventAssemblyAlg(), execute(), and initialize().

bool EventAssemblyAlg::myActiveCutFlag [private]

Definition at line 32 of file EventAssemblyAlg.h.

Referenced by EventAssemblyAlg(), and execute().

double EventAssemblyAlg::myEmcPhiCut [private]

Definition at line 35 of file EventAssemblyAlg.h.

Referenced by EventAssemblyAlg(), execute(), initialize(), and phiCut().

double EventAssemblyAlg::myEmcPhiNSigmaCut [private]

Definition at line 38 of file EventAssemblyAlg.h.

Referenced by EventAssemblyAlg(), execute(), and initialize().

double EventAssemblyAlg::myEmcThetaCut [private]

Definition at line 34 of file EventAssemblyAlg.h.

Referenced by EventAssemblyAlg(), execute(), initialize(), and thetaCut().

double EventAssemblyAlg::myEmcThetaNSigmaCut [private]

Definition at line 37 of file EventAssemblyAlg.h.

Referenced by EventAssemblyAlg(), execute(), and initialize().

NTuple::Item<double> EventAssemblyAlg::myEnergyMatched [private]

Definition at line 45 of file EventAssemblyAlg.h.

Referenced by execute(), and initialize().

NTuple::Item<double> EventAssemblyAlg::myExted [private]

Definition at line 43 of file EventAssemblyAlg.h.

Referenced by execute(), and initialize().

NTuple::Item<double> EventAssemblyAlg::myMatched [private]

Definition at line 44 of file EventAssemblyAlg.h.

Referenced by execute(), and initialize().

bool EventAssemblyAlg::myMsgFlag [private]

Definition at line 31 of file EventAssemblyAlg.h.

Referenced by EventAssemblyAlg(), and execute().

NTuple::Tuple* EventAssemblyAlg::myNTuple1 [private]

Definition at line 42 of file EventAssemblyAlg.h.

Referenced by execute(), and initialize().

int EventAssemblyAlg::myParticleId [private]

Definition at line 39 of file EventAssemblyAlg.h.

Referenced by EventAssemblyAlg(), and execute().


Generated on Tue Nov 29 23:18:49 2016 for BOSS_7.0.2 by  doxygen 1.4.7