StRoot: StTofPool/StTofpMcAnalysisMaker/StTofpMcAnalysisMaker.cxx Source File

00001 /*************************************************
00002  *
00003  * $Id: StTofpMcAnalysisMaker.cxx,v 1.4 2011/04/03 15:52:57 fisyak Exp $
00004  * $Log: StTofpMcAnalysisMaker.cxx,v $
00005  * Revision 1.4  2011/04/03 15:52:57  fisyak
00006  * Fix effect of constness in StAssociationMaker
00007  *
00008  * Revision 1.3  2007/04/17 23:11:04  dongx
00009  * replaced with standard STAR Loggers
00010  *
00011  * Revision 1.2  2005/10/06 19:58:15  fisyak
00012  * Adjust for persistent StMcEvent
00013  *
00014  * Revision 1.1  2004/03/16 04:58:55  geurts
00015  * *** empty log message ***
00016  *
00017  *
00018  *************************************************/
00019 
00020 #include <Stiostream.h>
00021 #include "TH1.h"
00022 #include "TH2.h"
00023 #include "StThreeVectorF.hh"
00024 #include "StThreeVectorD.hh"
00025 #include "StEventTypes.h"
00026 #include "StMcEventTypes.hh"
00027 #include "StAssociationMaker/StAssociationMaker.h"
00028 #include "StAssociationMaker/StTrackPairInfo.hh"
00029 #include "StEventMaker/StEventMaker.h"
00030 #include "StTofUtil/StTofGeometry.h"
00031 #include "StTofUtil/tofPathLength.hh"
00032 #include "StMessMgr.h"
00033 #include "StTofpMcAnalysisMaker.h"
00034 
00035 ClassImp(StTofpMcAnalysisMaker)
00036 
00037 
00038 const Int_t   StTofpMcAnalysisMaker::mPtBin = 50;
00039 const Int_t   StTofpMcAnalysisMaker::mYBin  = 14;
00040 const Float_t StTofpMcAnalysisMaker::mPtMin = 0.;
00041 const Float_t StTofpMcAnalysisMaker::mPtMax = 5.;
00042 const Float_t StTofpMcAnalysisMaker::mYMin  = -1.2;
00043 const Float_t StTofpMcAnalysisMaker::mYMax  =  0.2;
00044 
00045 
00046 
00047 //----------------------------------------------------------------
00048 const StGlobalTrack*  partnerTrack(mcTrackMapType* map, StMcTrack* mT) {
00049   mcTrackMapIter i = map->find(mT);
00050   const StGlobalTrack* rT = 0;
00051   if (i != map->end()) rT = (*i).second->partnerTrack();
00052   return rT;
00053 }
00054 
00055 
00056 
00057 //_________________________________________________
00058 StTofpMcAnalysisMaker::StTofpMcAnalysisMaker(const char *name, const char *title):StMaker(name,title)
00059 {
00060   //  StTofpMcAnalysisMaker Constructor
00061 
00062   mOuterTrackGeometry = true;
00063   mMinHitsPerTrack = 0;
00064   
00065 
00066   // - zero all pointers defined in the header file
00067   hMcPionPlus      = 0;
00068   hMcPionMin       = 0;
00069   hMcKaonPlus      = 0;
00070   hMcKaonMin       = 0;
00071   hMcProton        = 0;
00072   hMcAntiProton    = 0;
00073   hMcElectron      = 0;
00074   hRcPionPlus      = 0;
00075   hRcPionMin       = 0;
00076   hRcKaonPlus      = 0;
00077   hRcKaonMin       = 0;
00078   hRcProton        = 0;
00079   hRcAntiProton    = 0;
00080   hRcElectron      = 0;  
00081   hMatchPionPlus   = 0;
00082   hMatchPionMin    = 0;
00083   hMatchKaonPlus   = 0;
00084   hMatchKaonMin    = 0;
00085   hMatchProton     = 0;
00086   hMatchAntiProton = 0;
00087   hMatchElectron   = 0;
00088   //-
00089   hMomResPtPion = 0;
00090   hMultRef = 0;
00091   //-
00092   hTofpHitMap1  = 0;
00093   hTofpHitMap2  = 0;
00094   hTofpHitMap3  = 0;
00095   hTofpHitMap4  = 0;
00096   hTofpSlatIdA0 = 0;
00097   hTofpSlatIdA1 = 0;
00098   hTofpSlatIdB1 = 0;
00099   hTofpSlatIdD1 = 0;
00100   hTofpSlatIdD2 = 0;
00101   hTofpSlatIdE1 = 0; 
00102   hTofpSlatIdE2 = 0;
00103   hTofpSlatIdE3 = 0;
00104   hTofpSlatIdE4 = 0;
00105   hTofpSlatIdE5 = 0;
00106   hTofpSlatIdF1 = 0;
00107   hTofpSlatHitVecSize = 0;
00108 }
00109 
00110 //_________________________________________________
00111 StTofpMcAnalysisMaker::~StTofpMcAnalysisMaker()
00112 {
00113   //  StTofpMcAnalysisMaker Destructor
00114   if (Debug()) LOG_INFO << "Inside StTofpMcAnalysisMaker Destructor" << endm;
00115 }
00116 
00117 //_____________________________________________________________________________
00118 
00119 void StTofpMcAnalysisMaker::Clear(const char*)
00120 {
00121   StMaker::Clear();
00122 }
00123 
00124 //_________________________________________________
00125 Int_t StTofpMcAnalysisMaker::Finish()
00126 {
00127   return StMaker::Finish();
00128 }
00129 
00130 
00131 //_________________________________________________
00132 Int_t StTofpMcAnalysisMaker::Init()
00133 {
00134   // StTofpMcAnalysisMaker - Init
00135   
00136   hMomResPtPion = new TH2F("momResMomTPion","(|pmc| - |p|)/|pmc| as a Function of Pt Pion",mPtBin,mPtMin,mPtMax,200,-1.,1.);
00137   hMomResPtPion->SetXTitle("Pt (GeV/c)");
00138   hMomResPtPion->SetYTitle("(|pmc| - |p|)/|pmc| (GeV/c)");
00139 
00140   hMultRef = new TH1F("multRef","Mult Ref",100,0.,1000.);
00141   hMultRef->SetXTitle("uncorrected negative primary tracks");
00142 
00143   //-
00144   hMcElectron = new TH2F("mcElectron","Mc Electron",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00145   hMcElectron->SetXTitle("p_{T} (GeV/c)");
00146   hMcElectron->SetYTitle("pseudorapidity");
00147   hMcElectron->Sumw2();
00148   hRcElectron = new TH2F("rcElectron","RcElectron",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00149   hRcElectron->SetXTitle("p_{T} (GeV/c)");
00150   hRcElectron->SetYTitle("pseudorapidity");
00151   hRcElectron->Sumw2();
00152   hMatchElectron = new TH2F("mElectron","TOF matched electron",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00153   hMatchElectron->Sumw2();
00154   hMatchElectron->SetXTitle("p_{T} (GeV/c)");
00155   hMatchElectron->SetYTitle("pseudorapidity");
00156 
00157   //-
00158   hMcPionPlus = new TH2F("mcPionPlus","Mc Pion Plus",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00159   hMcPionPlus->SetXTitle("p_{T} (GeV/c)");
00160   hMcPionPlus->SetYTitle("pseudorapidity");
00161   hMcPionPlus->Sumw2();
00162   hRcPionPlus = new TH2F("rcPionPlus","Rc Pion Plus",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00163   hRcPionPlus->SetXTitle("p_{T} (GeV/c)");
00164   hRcPionPlus->SetYTitle("pseudorapidity");
00165   hRcPionPlus->Sumw2();
00166   hMatchPionPlus = new TH2F("mPionPlus","TOF matched pi+",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00167   hMatchPionPlus->Sumw2();
00168   hMatchPionPlus->SetXTitle("p_{T} (GeV/c)");
00169   hMatchPionPlus->SetYTitle("pseudorapidity");
00170 
00171   //-
00172   hMcPionMin = new TH2F("mcPionMin","Mc Pion Min",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00173   hMcPionMin->SetXTitle("p_{T} (GeV/c)");
00174   hMcPionMin->SetYTitle("pseudorapidity");
00175   hMcPionMin->Sumw2();
00176   hRcPionMin = new TH2F("rcPionMin","Rc Pion Min",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00177   hRcPionMin->SetXTitle("p_{T} (GeV/c)");
00178   hRcPionMin->SetYTitle("pseudorapidity");
00179   hRcPionMin->Sumw2();
00180   hMatchPionMin = new TH2F("mPionMin","TOF matched pi-",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00181   hMatchPionMin->Sumw2();
00182   hMatchPionMin->SetXTitle("p_{T} (GeV/c)");
00183   hMatchPionMin->SetYTitle("pseudorapidity");
00184 
00185   //-
00186   hMcKaonPlus = new TH2F("mcKaonPlus","Mc Kaon Plus",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00187   hMcKaonPlus->SetXTitle("p_{T} (GeV/c)");
00188   hMcKaonPlus->SetYTitle("pseudorapidity");
00189   hMcKaonPlus->Sumw2();
00190   hRcKaonPlus = new TH2F("rcKaonPlus","Rc Kaon Plus",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00191   hRcKaonPlus->SetXTitle("p_{T} (GeV/c)");
00192   hRcKaonPlus->SetYTitle("pseudorapidity");
00193   hRcKaonPlus->Sumw2();
00194   hMatchKaonPlus = new TH2F("mKaonPlus","TOF matched K+",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00195   hMatchKaonPlus->Sumw2();
00196   hMatchKaonPlus->SetXTitle("p_{T} (GeV/c)");
00197   hMatchKaonPlus->SetYTitle("pseudorapidity");
00198 
00199   //-
00200   hMcKaonMin = new TH2F("mcKaonMin","Mc Kaon Min",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00201   hMcKaonMin->SetXTitle("p_{T} (GeV/c)");
00202   hMcKaonMin->SetYTitle("pseudorapidity");
00203   hMcKaonMin->Sumw2();
00204   hRcKaonMin = new TH2F("rcKaonMin","Rc Kaon Min",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00205   hRcKaonMin->SetXTitle("p_{T} (GeV/c)");
00206   hRcKaonMin->SetYTitle("pseudorapidity");
00207   hRcKaonMin->Sumw2();
00208   hMatchKaonMin = new TH2F("mKaonMin","TOF matched K-",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00209   hMatchKaonMin->Sumw2();
00210   hMatchKaonMin->SetXTitle("p_{T} (GeV/c)");
00211   hMatchKaonMin->SetYTitle("pseudorapidity");
00212 
00213   //-
00214   hMcProton = new TH2F("mcProton","Mc Proton",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00215   hMcProton->SetXTitle("p_{T} (GeV/c)");
00216   hMcProton->SetYTitle("pseudorapidity");
00217   hMcProton->Sumw2();
00218   hRcProton = new TH2F("rcProton","RcProton",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00219   hRcProton->SetXTitle("p_{T} (GeV/c)");
00220   hRcProton->SetYTitle("pseudorapidity");
00221   hRcProton->Sumw2();
00222   hMatchProton = new TH2F("mProton","TOF matched proton",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00223   hMatchProton->Sumw2();
00224   hMatchProton->SetXTitle("p_{T} (GeV/c)");
00225   hMatchProton->SetYTitle("pseudorapidity");
00226 
00227   //-
00228   hMcAntiProton = new TH2F("mcAntiProton","Mc Antiproton",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00229   hMcAntiProton->SetXTitle("p_{T} (GeV/c)");
00230   hMcAntiProton->SetYTitle("pseudorapidity");
00231   hMcAntiProton->Sumw2();
00232   hRcAntiProton = new TH2F("rcAntiProton","Rc Antiproton",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00233   hRcAntiProton->SetXTitle("p_{T} (GeV/c)");
00234   hRcAntiProton->SetYTitle("pseudorapidity");
00235   hRcAntiProton->Sumw2();
00236   hMatchAntiProton = new TH2F("mAntiProton","TOF matched pbar",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
00237   hMatchAntiProton->Sumw2();
00238   hMatchAntiProton->SetXTitle("p_{T} (GeV/c)");
00239   hMatchAntiProton->SetYTitle("pseudorapidity");
00240 
00241   // Match QA Histograms
00242   hTofpSlatIdA0 = new TH1D("SlatIdA0","events per slat",41,0.5,41.5);
00243   hTofpSlatIdA1 = new TH1D("SlatIdA1","valid slat",41,0.5,41.5);
00244   hTofpSlatIdB1 = new TH1D("SlatIdB1","#tracks match  valid slat",41,0.5,41.5);
00245   hTofpSlatIdD1 = new TH1D("SlatIdD1","track match per valid slat",41,0.5,41.5);
00246   hTofpSlatIdD2 = new TH1D("SlatIdD2","single track match per slat",41,0.5,41.5);
00247   hTofpSlatIdE1 = new TH1D("SlatIdE1","one slat for one track match",41,0.5,41.5);
00248   hTofpSlatIdE2 = new TH1D("SlatIdE2","recovered from hitprof-weight",41,0.5,41.5);
00249   hTofpSlatIdE3 = new TH1D("SlatIdE3","recovered from ss",41,0.5,41.5);
00250   hTofpSlatIdE4 = new TH1D("SlatIdE4","recovered from closest hitplane",41,0.5,41.5);
00251   hTofpSlatIdE5 = new TH1D("SlatIdE5","total recovered slat per track match",41,0.5,41.5);
00252   hTofpSlatIdF1 = new TH1D("SlatIdF1","primary track match per slat",41,0.5,41.5);
00253   hTofpSlatHitVecSize = new TH1D("SlatMult","Slat Mult per Track",10,0,10);
00254   hTofpHitMap1 = new TH2D("tofpHitMap1","valid hit positions", 500,-250,0, 120,-1.23, -1.08);
00255   hTofpHitMap2 = new TH2D("tofpHitMap2","valid hit positions", 500,-250,0, 120,-1.23, -1.08);
00256   hTofpHitMap3 = new TH2D("tofpHitMap3","valid hit positions", 500,-250,0, 120,-1.23, -1.08);
00257   hTofpHitMap4 = new TH2D("tofpHitMap4","valid hit positions", 500,-250,0, 120,-1.23, -1.08);
00258 
00259   return StMaker::Init();
00260 }
00261 
00262 
00263 Int_t StTofpMcAnalysisMaker::InitRun(int runnumber){
00264   gMessMgr->Info("StTofpMcAnalysisMaker -- reinitializing TofGeometry","OS");
00265   mTofGeom = new StTofGeometry();
00266   mTofGeom->init(this);
00267   return kStOk;
00268 }
00269 
00270 Int_t StTofpMcAnalysisMaker::FinishRun(int runnumber){
00271   gMessMgr->Info("StTofpMcAnalysisMaker -- cleaning up geometry","OS" );
00272   if (mTofGeom) delete mTofGeom;
00273   mTofGeom=0;
00274   return kStOk;
00275 }
00276 
00277 //_________________________________________________
00278 Int_t StTofpMcAnalysisMaker::Make(){
00279 
00280   // Get the pointers we need, we have to use the titles we gave them in the
00281   // macro.  I just used the defaults.
00282 
00283   // StEvent
00284   StEvent* rEvent = 0;
00285   rEvent = (StEvent*) GetInputDS("StEvent");
00286 
00287   // StMcEvent
00288   StMcEvent* mEvent = (StMcEvent*) GetDataSet("StMcEvent");
00289 
00290   // StAssociationMaker
00291   StAssociationMaker* assoc = 0;
00292   assoc = (StAssociationMaker*) GetMaker("StAssociationMaker");
00293 
00294   // the Multimaps...
00295   rcTpcHitMapType* theHitMap = 0;
00296   theHitMap = assoc->rcTpcHitMap();
00297   rcTrackMapType* theTrackMap = 0;
00298   theTrackMap = assoc->rcTrackMap();
00299   mcV0MapType* theMcV0Map = 0;
00300   theMcV0Map = assoc->mcV0Map(); 
00301 
00302   if (!theHitMap) {
00303     gMessMgr->Warning() << "----------WARNING----------\n"
00304                         << "No Hit Map found for this event!" << endm;
00305     return kStWarn;
00306   }
00307 
00308   mcTrackMapType* mcTrackMap = 0;
00309   mcTrackMap = assoc->mcTrackMap();
00310   if (!mcTrackMap){
00311     LOG_INFO << "No mcTrackMap!!! " << endm;
00312     return 0;
00313   }
00314 
00315 
00316   // Example: look at the position of the primary vertex
00317   //          Map is not needed for this, but it's a good check,
00318   //          tracking will not be good if primary vertex was not well placed.
00319 
00320 
00321   // First check whether the Primary Vertex is there at all.
00322   StThreeVectorD VertexPos(0,0,0);
00323   if (rEvent->primaryVertex()) {
00324     VertexPos = rEvent->primaryVertex()->position();
00325     LOG_INFO << "Position of Primary Vertex from StEvent:" << endm;
00326     LOG_INFO << VertexPos << endm;
00327   }
00328   else {
00329     LOG_INFO << "----------- WARNING ------------" << endm;
00330     LOG_INFO << "No primary vertex from StEvent!!" << endm;
00331     LOG_INFO << "Assume it is at (0,0,0)         " << endm;
00332   }
00333   LOG_INFO << "Position of Primary Vertex from StMcEvent:" << endm;
00334   LOG_INFO << mEvent->primaryVertex()->position() << endm;
00335     
00336   if (!theTrackMap) {
00337     gMessMgr->Warning() << "----------WARNING----------\n"
00338                         << "No Track Map found for this event!" << endm;
00339     return kStWarn;
00340   }
00341 
00342   // -----------------------------------------
00343   // --- Stage I: MATCH RC-TRACKS to SLATS ---
00344   // -----------------------------------------
00345 
00346   // Note: why not replaced Stage I w/ the actual TOFpMatchMaker and
00347   //       use the StTofSlat in StEvent as the reference for Stage IIb?
00348 
00349 
00350   StEvent *event = rEvent;
00351   bool mHisto(true);
00352   int NTOFP(41);
00353   //.........................................................................
00354   // A. build vector of candidate slats with valid ADC signals 
00355   idVector validSlatIdVec;
00356   for (int i=0;i<NTOFP;i++){
00357     unsigned short slatId = mTofGeom->daqToSlatId(i);
00358     hTofpSlatIdA0->Fill(i+1);
00359 
00360     validSlatIdVec.push_back(slatId);
00361     hTofpSlatIdA1->Fill(i+1);
00362   }
00363   gMessMgr->Info("","OST") << "A: #valid slats: " << validSlatIdVec.size() << endm;
00364   // end of Sect.A
00365 
00366 
00367   //.........................................................................
00368   // B. loop over global tracks and determine all slat-track matches
00369   //
00370   tofSlatHitVector allSlatsHitVec;
00371   allSlatsHitVec.clear();
00372   StSPtrVecTrackNode& nodes = event->trackNodes();
00373   int nAllTracks=0;
00374   for (unsigned int iNode=0; iNode<nodes.size(); iNode++){
00375     tofSlatHitVector slatHitVec;
00376     slatHitVec.clear();
00377     StTrack *theTrack = nodes[iNode]->track(global);
00378 
00379     // make sure we have a track, a miniDST might have removed it...
00380     if (validTrack(theTrack)){
00381       nAllTracks++;
00382       StPhysicalHelixD theHelix = trackGeometry(theTrack)->helix();
00383       slatHitVec = mTofGeom->tofHelixToArray(theHelix, validSlatIdVec);
00384       if (slatHitVec.size()>0 && mHisto) hTofpSlatHitVecSize->Fill(slatHitVec.size());
00385 
00386       for (size_t ii = 0; ii < slatHitVec.size(); ii++) {
00387         slatHitVec[ii].trackIdVec.push_back(iNode);
00388         allSlatsHitVec.push_back(slatHitVec[ii]);
00389         if (mHisto){
00390           int id = mTofGeom->slatIdToDaq(slatHitVec[ii].slatIndex);
00391           float xhit = slatHitVec[ii].hitPosition.x();
00392           float yhit = slatHitVec[ii].hitPosition.y();
00393           float zhit = slatHitVec[ii].hitPosition.z();
00394           float phihit = atan2(yhit,xhit);
00395           hTofpSlatIdB1->Fill(id);
00396           hTofpHitMap1->Fill(zhit,phihit);
00397         }
00398 
00399         if (Debug()){ 
00400           LOG_INFO << "B: trackid=";
00401           idVectorIter ij = slatHitVec[ii].trackIdVec.begin();
00402           while (ij != slatHitVec[ii].trackIdVec.end()) {LOG_INFO << " " << *ij; ij++;}
00403           LOG_INFO << "\tind=" << mTofGeom->slatIdToDaq(slatHitVec[ii].slatIndex)
00404                << "\thitprof="<< slatHitVec[ii].hitProfile 
00405                << "\ts="<<slatHitVec[ii].s << "\tthxy="<<slatHitVec[ii].theta_xy 
00406                << "\tthzr="<<slatHitVec[ii].theta_zr;
00407           if (slatHitVec.size()>1) LOG_INFO << " M" << endm;
00408           else LOG_INFO << endm;
00409         }
00410       }
00411     } // existing global track
00412   } // loop over nodes
00413 
00414   gMessMgr->Info("","OST") << "B: #matched/#avail/#total tracknodes: " 
00415                            <<allSlatsHitVec.size() << "/" <<nAllTracks 
00416                            << "/" << nodes.size() << endm;
00417   // end of Sect.B
00418 
00419 
00420   //.........................................................................
00421   // C Neighbours -- identify crosstalk, geometry shifts (only fill histograms)
00422   //
00423   // [ ... ]
00424   // end of Sect.C
00425 
00426 
00427   //.........................................................................
00428   // D. sort hit vectors  and deal with (discard) slats matched by multiple tracks
00429   //
00430   int nSingleHitSlats(0);
00431   tofSlatHitVector singleHitSlatsVec;
00432   StructSlatHit slatHit;
00433 
00434   tofSlatHitVector tempVec = allSlatsHitVec;
00435   tofSlatHitVector erasedVec = tempVec;
00436   while (tempVec.size() != 0) {
00437     int nTracks = 0;
00438     vector<StThreeVectorD> vPosition;
00439     vector<vector<StThreeVectorD> > vLayerHitPositions;
00440     vector<Int_t> vHitProfile;
00441     vector<Float_t> vS, vTheta_xy, vTheta_zr;
00442     idVector trackIdVec;
00443 
00444     tofSlatHitVectorIter tempIter=tempVec.begin();
00445     tofSlatHitVectorIter erasedIter=erasedVec.begin();
00446     while(erasedIter!= erasedVec.end()) {
00447       if(tempIter->slatIndex == erasedIter->slatIndex) {
00448         nTracks++;
00449         // save all hit data in temporary vectors
00450         trackIdVec.push_back(erasedIter->trackIdVec.back());
00451         vPosition.push_back(erasedIter->hitPosition);
00452         vLayerHitPositions.push_back(erasedIter->layerHitPositions);
00453         vHitProfile.push_back(erasedIter->hitProfile);
00454         vS.push_back(erasedIter->s);
00455         vTheta_xy.push_back(erasedIter->theta_xy);
00456         vTheta_zr.push_back(erasedIter->theta_zr);
00457 
00458         if (mHisto){
00459          float xhit = erasedIter->hitPosition.x();
00460          float yhit = erasedIter->hitPosition.y();
00461          float zhit = erasedIter->hitPosition.z();
00462          float phihit = atan2(yhit,xhit);
00463          hTofpHitMap2->Fill(zhit,phihit);
00464         }
00465 
00466         erasedVec.erase(erasedIter);
00467         erasedIter--;
00468       }
00469       erasedIter++;
00470     }
00471 
00472     if (mHisto){
00473       int daqId = mTofGeom->slatIdToDaq(tempIter->slatIndex);
00474       hTofpSlatIdD1->Fill(daqId);
00475     }
00476 
00477     if (nTracks==1){
00478       nSingleHitSlats++;
00479       // for singly hit slat, copy data in singleHitSlatsVec
00480       slatHit.slatIndex = tempIter->slatIndex;
00481       slatHit.hitPosition = vPosition[0];
00482       slatHit.layerHitPositions = vLayerHitPositions[0];
00483       slatHit.trackIdVec = trackIdVec;
00484       slatHit.hitProfile =  vHitProfile[0];
00485       slatHit.s    =  vS[0];
00486       slatHit.theta_xy =  vTheta_xy[0];
00487       slatHit.theta_zr =  vTheta_zr[0];
00488 
00489       singleHitSlatsVec.push_back(slatHit);
00490 
00491       if (mHisto){
00492         int daqId = mTofGeom->slatIdToDaq(tempIter->slatIndex);
00493         float xhit = slatHit.hitPosition.x();
00494         float yhit = slatHit.hitPosition.y();
00495         float zhit = slatHit.hitPosition.z();
00496         float phihit = atan2(yhit,xhit);
00497         hTofpSlatIdD2->Fill(daqId);
00498         hTofpHitMap3->Fill(zhit,phihit);
00499       }
00500 
00501       // debugging output
00502       if (Debug()) {
00503         LOG_INFO << "D: ind=" << mTofGeom->slatIdToDaq(slatHit.slatIndex)
00504              << "\thitprof="<< slatHit.hitProfile << "\ts="<<slatHit.s
00505              << "\tthxy="<<slatHit.theta_xy << "\tthzr="<<slatHit.theta_zr << "\ttrackid:";
00506         idVectorIter ij=trackIdVec.begin();
00507         while (ij != trackIdVec.end()) { LOG_INFO << " " << *ij; ij++; }
00508         LOG_INFO <<endm;
00509       }
00510     }
00511     else if (nTracks>1){
00512       // for multiple hit slats either discard (yes) or
00513       // find the most likely candidate.
00514     } else
00515        gMessMgr->Warning("","OST")  << "D: no tracks extrapolate to matched slat ... should not happen!" << endm;
00516 
00517     tempVec = erasedVec;
00518   }
00519    gMessMgr->Info("","OST") << "D: #before/#after: " << allSlatsHitVec.size() 
00520                             << "/" << singleHitSlatsVec.size() << endm;
00521   //end of Sect.D
00522 
00523 
00524 
00525   //.........................................................................
00526   // E. sort and deal singleHitSlatsVector for multiple slats associated to single tracks
00527   //
00528   tofSlatHitVector allMatchedSlatsVec;
00529   tempVec = singleHitSlatsVec;
00530   erasedVec = tempVec;
00531   while (tempVec.size() != 0) {
00532     StructSlatHit slatHit;
00533     int nSlats = 0;
00534     vector<StThreeVectorD> vPosition;
00535     vector< vector<StThreeVectorD> > vLayerHitPositions;
00536     vector<Int_t> vHitProfile;
00537     vector<Float_t> vS, vTheta_xy, vTheta_zr;
00538     idVector vTrackId;
00539     vector<Int_t> slatIndex;
00540 
00541     tofSlatHitVectorIter tempIter=tempVec.begin();
00542     tofSlatHitVectorIter erasedIter=erasedVec.begin();
00543     while(erasedIter!= erasedVec.end()) {
00544       if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
00545         nSlats++;
00546         // save all hit data in temporary vectors
00547         slatIndex.push_back(erasedIter->slatIndex);
00548         vTrackId.push_back(erasedIter->trackIdVec.back());
00549         vPosition.push_back(erasedIter->hitPosition);
00550         vLayerHitPositions.push_back(erasedIter->layerHitPositions);
00551         vHitProfile.push_back(erasedIter->hitProfile);
00552         vS.push_back(erasedIter->s);
00553         vTheta_xy.push_back(erasedIter->theta_xy);
00554         vTheta_zr.push_back(erasedIter->theta_zr);
00555 
00556         erasedVec.erase(erasedIter);
00557         erasedIter--;
00558       }
00559       erasedIter++;
00560     }
00561 
00562     
00563     if (nSlats==1){
00564       // for singly hit slat, copy data in singleHitSlatsVec
00565       slatHit.slatIndex = slatIndex[0];
00566       slatHit.hitPosition = vPosition[0];
00567       slatHit.layerHitPositions = vLayerHitPositions[0];
00568       slatHit.trackIdVec.push_back(vTrackId[0]);
00569       slatHit.hitProfile =  vHitProfile[0];
00570       slatHit.s    =  vS[0];
00571       slatHit.theta_xy =  vTheta_xy[0];
00572       slatHit.theta_zr =  vTheta_zr[0];
00573       slatHit.matchFlag = 0; 
00574 
00575       allMatchedSlatsVec.push_back(slatHit);
00576 
00577       if (mHisto){
00578         int daqId = mTofGeom->slatIdToDaq(slatIndex[0]);
00579         hTofpSlatIdE1->Fill(daqId);
00580       }
00581 
00582       // debugging output
00583       if (Debug()) {
00584         LOG_INFO << "E: ind=" << mTofGeom->slatIdToDaq(slatHit.slatIndex)
00585              << "\thitprof="<< slatHit.hitProfile << "\ts="<<slatHit.s
00586              << "\tthxy="<<slatHit.theta_xy << "\tthzr="<<slatHit.theta_zr << "\ttrackid:";
00587         idVectorIter ij=vTrackId.begin();
00588         while (ij != vTrackId.end()) { LOG_INFO << " " << *ij; ij++; }
00589         LOG_INFO <<endm;
00590       }
00591     }
00592     else if (nSlats>1){   // for multiple hit slats  find the most likely candidate.
00593       int thiscandidate(-99);
00594       int thisMatchFlag(0);
00595 
00596       // 1. sort on hitprofile weight
00597       int weight(0);
00598       vector<int> weightCandidates;
00599       thisMatchFlag = 1;
00600       if (Debug()) LOG_INFO << "E: find ... weight ";
00601       for (int i=0;i<nSlats;i++){
00602         int hitWeight = vLayerHitPositions[i].size();
00603         if (Debug()) LOG_INFO << mTofGeom->slatIdToDaq(slatIndex[i]) << "("<<hitWeight<<")"<<" ";
00604         if (hitWeight>weight) {
00605           weight=hitWeight;
00606           weightCandidates.clear();
00607           weightCandidates.push_back(i);
00608         } else if (hitWeight == weight)
00609           weightCandidates.push_back(i);
00610       }
00611       if (weightCandidates.size()==1){
00612         thiscandidate = weightCandidates[0];
00613         int daqId = mTofGeom->slatIdToDaq(slatIndex[thiscandidate]);
00614         if (mHisto) hTofpSlatIdE2->Fill(daqId);
00615         if (Debug()) LOG_INFO << "candidate =" << daqId << endm;
00616       }
00617 
00618       // 2. if still undecided check on ss
00619       if (weightCandidates.size()>1){
00620         Float_t ss(0);
00621         vector<int> ssCandidates;
00622         thisMatchFlag = 2;
00623         if (Debug()) LOG_INFO << " ss ";
00624         for (unsigned int i=0;i<weightCandidates.size();i++){
00625           int ii=weightCandidates[i];     
00626           if (Debug()) LOG_INFO << mTofGeom->slatIdToDaq(slatIndex[ii]) << " ";
00627           if (vS[ii]>ss){
00628             ss = vS[ii];
00629             ssCandidates.clear();
00630             ssCandidates.push_back(ii);
00631           }else if  (vS[ii]==ss)
00632             ssCandidates.push_back(ii);   
00633         }
00634         if (ssCandidates.size()==1){
00635           thiscandidate = ssCandidates[0];
00636           int daqId = mTofGeom->slatIdToDaq(slatIndex[thiscandidate]);
00637           if (mHisto) hTofpSlatIdE3->Fill(daqId);
00638           if (Debug()) LOG_INFO << "candidate =" << daqId << endm;
00639         }
00640 
00641         // 3. if still undecided go for closest/first hit
00642         if (ssCandidates.size()>1){
00643           Int_t hitprof(0);
00644           vector<int> profileCandidates;
00645           thisMatchFlag = 3;
00646           if (Debug()) LOG_INFO << " hprof ";
00647           for (unsigned int i=0;i<ssCandidates.size();i++){
00648             int ii=ssCandidates[i];
00649             if (Debug()) LOG_INFO << mTofGeom->slatIdToDaq(slatIndex[ii]) << " ";
00650             if (vHitProfile[ii]>hitprof){
00651               hitprof = vHitProfile[ii];
00652               profileCandidates.clear();
00653               profileCandidates.push_back(ii);
00654           }else if  (vHitProfile[ii]==hitprof)
00655             profileCandidates.push_back(ii);      
00656           } 
00657           if (profileCandidates.size()==1){
00658             thiscandidate = profileCandidates[0];
00659             int daqId = mTofGeom->slatIdToDaq(slatIndex[thiscandidate]);
00660             if (mHisto) hTofpSlatIdE4->Fill(daqId);
00661             if (Debug()) LOG_INFO << "candidate =" << daqId << endm;
00662           }
00663           else
00664             if (Debug()) LOG_INFO << "none" << endm;
00665         }
00666 
00667 
00668         // forget it, and let user know of the non-decision
00669         if (thiscandidate == -99 && Debug()){
00670           LOG_INFO << "E: ind=";
00671           for (unsigned int ii=0;ii<slatIndex.size();ii++) 
00672             LOG_INFO << mTofGeom->slatIdToDaq(slatIndex[ii]) << " ";
00673           LOG_INFO << "\ttrkid:" << vTrackId[0] << " Unable to decide. "; 
00674           LOG_INFO << "(hitprofs:";
00675           for (unsigned int ii=0;ii<slatIndex.size();ii++) 
00676             LOG_INFO << vHitProfile[ii] << " ";
00677           LOG_INFO << " ss:";
00678           for (unsigned int ii=0;ii<slatIndex.size();ii++) 
00679             LOG_INFO << vS[ii] << " ";
00680           LOG_INFO << ")" << endm;
00681         }
00682 
00683       }
00684 
00685 
00686       if (thiscandidate>=0){
00687         slatHit.slatIndex = slatIndex[thiscandidate];
00688         slatHit.hitPosition = vPosition[thiscandidate];
00689         slatHit.layerHitPositions = vLayerHitPositions[thiscandidate];
00690         //slatHit.trackIdVec.clear();
00691         slatHit.trackIdVec.push_back(vTrackId[thiscandidate]);
00692         slatHit.hitProfile =  vHitProfile[thiscandidate];
00693         slatHit.s    =  vS[thiscandidate];
00694         slatHit.theta_xy =  vTheta_xy[thiscandidate];
00695         slatHit.theta_zr =  vTheta_zr[thiscandidate];
00696         slatHit.matchFlag = thisMatchFlag;
00697         
00698         allMatchedSlatsVec.push_back(slatHit);
00699 
00700         if (mHisto){
00701           int daqId = mTofGeom->slatIdToDaq(slatIndex[thiscandidate]);
00702           hTofpSlatIdE5->Fill(daqId);
00703         }
00704         
00705         // debugging output
00706         if (Debug()) {
00707           LOG_INFO << "E: ind=" << mTofGeom->slatIdToDaq(slatHit.slatIndex)
00708                << "\thitprof="<< slatHit.hitProfile << "\ts="<<slatHit.s
00709                << "\tthxy="<<slatHit.theta_xy << "\tthzr="<<slatHit.theta_zr << "\ttrackid:"
00710                << vTrackId[thiscandidate] << endm;
00711         }
00712       }
00713     } else
00714        gMessMgr->Warning("","OS")  << "E: no slats belong to this track ... should not happen!" << endm;
00715 
00716     tempVec = erasedVec;
00717   }
00718   gMessMgr->Info("","OST") << "E: #before/#after: " << singleHitSlatsVec.size()
00719                             << "/" << allMatchedSlatsVec.size() << endm;
00720   // end of Sect.E
00721   
00722 
00723   // --------------------------------------------------------------------//
00724   // Note: section F is not required: the primary-track requirement will //
00725   //       be enforced lateron in the MC-RC track association.           //
00726   // --------------------------------------------------------------------//
00727 
00728 
00729   //.........................................................................
00730   // F. perform further selection and
00731   //    fill valid track histograms and SlatCollection
00732   //
00733   //StTofSlatCollection *mSlatCollection =  new StTofSlatCollection;
00734   int nValidSinglePrimHitSlats(0); //nValidSingleHitSlats(0)
00735 
00736   for (size_t ii=0; ii < allMatchedSlatsVec.size(); ii++){
00737     int daqId = mTofGeom->slatIdToDaq(allMatchedSlatsVec[ii].slatIndex);
00738     //int jj = daqId-1;
00739 
00740     if (allMatchedSlatsVec[ii].trackIdVec.size()!=1)
00741        gMessMgr->Warning("","OST") << "F: WHAT!?!  mult.matched slat in single slat list " << daqId << endm;
00742 
00743     // 1. fill valid single track AND valid tdc histograms
00744     //if (validTdc(mTofpTdc[jj])) nValidSingleHitSlats++;
00745 
00746     // get track-id from slat hit vector
00747     unsigned int trackNode = allMatchedSlatsVec[ii].trackIdVec[0];
00748     StTrack *theTrack = nodes[trackNode]->track(primary);
00749 
00750     // 2. continue only if the (primary) track exists
00751     if (validTofTrack(theTrack)){
00752       nValidSinglePrimHitSlats++;
00753       if (mHisto){
00754         float xhit = allMatchedSlatsVec[ii].hitPosition.x();
00755         float yhit = allMatchedSlatsVec[ii].hitPosition.y();
00756         float zhit = allMatchedSlatsVec[ii].hitPosition.z();
00757         float phihit = atan2(yhit,xhit);
00758         hTofpHitMap4->Fill(zhit,phihit);
00759         hTofpSlatIdF1->Fill(daqId);
00760       }
00761           
00762       //--- store number of hits per track
00763       int nHitsPerTrack = theTrack->topologyMap().numberOfHits(kTpcId);
00764       //if (mHisto) hTofpNumberOfTrackHits->Fill(nHitsPerTrack);
00765           
00766       // select the apropriate track geometry
00767       const StTrackGeometry *theTrackGeometry = trackGeometry(theTrack);
00768 
00769       //--- get momentum from track
00770       const StThreeVectorF momentum = theTrackGeometry->momentum();
00771       //if (mHisto) hTofpPtTrack->Fill(momentum.perp());
00772             
00773       //--- calculate flight path
00774       double pathLength = tofPathLength(&event->primaryVertex()->position(), 
00775                                         &allMatchedSlatsVec[ii].hitPosition,
00776                                         theTrackGeometry->helix().curvature());
00777             
00778 
00779       //--- calculate local hit position on slat based first, last and middle plane
00780       //    (middle plane is the average of first and last plane, which is mathematically
00781       //     the same as SlatHitVec.hitPosition ... )
00782       StThreeVectorD *pInnerLayer, *pOuterLayer;
00783       pInnerLayer =  &(*(allMatchedSlatsVec[ii].layerHitPositions.begin()));
00784       pOuterLayer =  &(*(allMatchedSlatsVec[ii].layerHitPositions.end() - 1));
00785 
00786       //--- dig out from the dedx and rich pid traits (only for Debug mode)
00787       float dedx(0.), cherang(0);
00788       int dedx_np(0), cherang_nph(0);
00789       if (Debug()){
00790         StSPtrVecTrackPidTraits& traits = theTrack->pidTraits();
00791         for (unsigned int it=0;it<traits.size();it++){
00792           if (traits[it]->detector() == kTpcId){
00793             StDedxPidTraits *pid = dynamic_cast<StDedxPidTraits*>(traits[it]);
00794             if (pid && pid->method() ==kTruncatedMeanId){
00795               dedx    = pid->mean();
00796               dedx_np =  pid->numberOfPoints();
00797             }
00798           } else if  (traits[it]->detector() == kRichId){
00799             StRichPidTraits *pid = dynamic_cast<StRichPidTraits*>(traits[it]);
00800             if (pid){ 
00801               StRichSpectra* pidinfo = pid->getRichSpectra();
00802               if (pidinfo && pidinfo->getCherenkovPhotons()>2){
00803                 cherang     = pidinfo->getCherenkovAngle();
00804                 cherang_nph = pidinfo->getCherenkovPhotons();
00805               }
00806             }
00807           }
00808         }
00809       }
00810 
00811       //--- calculate local hit position on slat based on average hitposition
00812       //float localHitPos = mTofGeom->slatHitPosition(&allMatchedSlatsVec[ii].hitPosition);
00813 
00814       // Fill TOF Slat Collection
00815       //StTofSlat *tofSlat = new StTofSlat(jj+1,(int)mTofpAdc[jj],(int)mTofpTdc[jj],theTrack);
00816       //StTofSlat *tofSlat = new StTofSlat(jj+1,(int)mTofpAdc[jj],(int)mTofpTdc[jj],theTrack,
00817       //                             localHitPos, allMatchedSlatsVec[ii].hitProfile,
00818       //                             allMatchedSlatsVec[ii].matchFlag);
00819       //tofSlat->setPosition(allMatchedSlatsVec[ii].hitPosition);
00820       //mSlatCollection->push_back(tofSlat);
00821 
00822       // dump debug data
00823       if (Debug()){
00824         LOG_INFO << "F: ind=" << mTofGeom->slatIdToDaq(allMatchedSlatsVec[ii].slatIndex) 
00825              << "\ttrackid:";
00826         idVectorIter ij=allMatchedSlatsVec[ii].trackIdVec.begin();
00827         while (ij != allMatchedSlatsVec[ii].trackIdVec.end()) { LOG_INFO << " " << *ij; ij++; }
00828         LOG_INFO << "\tR=" << 1/(theTrackGeometry->helix().curvature())
00829              << "\tpT=" << momentum.perp() << "\tp=" << momentum.mag()
00830              << "\thits="<< nHitsPerTrack << "\ts="<< pathLength
00831              << "\t#fitp=" <<theTrack->fitTraits().numberOfFitPoints(kTpcId)
00832              << "\t#trkp=" <<theTrack->detectorInfo()->numberOfPoints(kTpcId)
00833              << " \tdedx=" << dedx 
00834              << " \tdca="<< theTrack->geometry()->helix().distance(event->primaryVertex()->position());
00835         if (cherang!=0) LOG_INFO  << " \trich="<< cherang << " (" << cherang_nph << ")";
00836         LOG_INFO << endm;
00837       }
00838 
00839     } // track exists 
00840   }
00841 
00842 
00843     
00844   // ---------------------------------------------
00845   // --- Stage II: FIND RC-MC pairs (TPC eff.) ---
00846   // ---------------------------------------------
00847    
00848   // ----
00849   // IIa. loop over mcTracks and check for associated rcTrack (TPC tracking eff.)
00850   // IIb. check matching results for rcTrack (TOFp matching eff.)
00851   // ----
00852 
00853 
00854   const StSPtrVecMcTrack& mcTracks = mEvent->primaryVertex()->daughters();
00855   if (Debug()) LOG_INFO << "mcTrack Loop: "<< mcTracks.size() << " entries" << endm;
00856   StMcTrack* mcTrack;
00857   const StTrack* rcTrack = 0;
00858 
00859   int nRecon(0), nMatch(0);
00860 
00861   //--- Loop over MC tracks
00862   for (StMcTrackConstIterator mcItr = mcTracks.begin(); mcItr  != mcTracks.end(); mcItr++) {
00863     mcTrack = *mcItr;
00864     if (!mcTrack) {
00865       gMessMgr->Warning() << "No McTrack?!? Should not happen!"<<endm;
00866       continue;
00867     }
00868 
00869 
00870     //-- Fill histograms with McTrack data 
00871     StThreeVectorF mcP = mcTrack->momentum();
00872     if      (mcTrack->particleDefinition()->pdgEncoding()== 11)   hMcElectron->Fill(mcP.perp(),mcP.pseudoRapidity());
00873     else if (mcTrack->particleDefinition()->pdgEncoding()== 211)  hMcPionPlus->Fill(mcP.perp(),mcP.pseudoRapidity());
00874     else if (mcTrack->particleDefinition()->pdgEncoding()==-211)  hMcPionMin->Fill(mcP.perp(),mcP.pseudoRapidity());
00875     else if (mcTrack->particleDefinition()->pdgEncoding()== 321)  hMcKaonPlus->Fill(mcP.perp(),mcP.pseudoRapidity());
00876     else if (mcTrack->particleDefinition()->pdgEncoding()==-321)  hMcKaonMin->Fill(mcP.perp(),mcP.pseudoRapidity());
00877     else if (mcTrack->particleDefinition()->pdgEncoding()== 2212) hMcProton->Fill(mcP.perp(),mcP.pseudoRapidity());
00878     else if (mcTrack->particleDefinition()->pdgEncoding()==-2212) hMcAntiProton->Fill(mcP.perp(),mcP.pseudoRapidity());
00879     else LOG_INFO << "mcTracks has wrong pdgEncoding: "<< mcTrack->particleDefinition()->pdgEncoding() << endm; 
00880 
00881     //-- get associated reconstructed rctrack
00882     rcTrack = partnerTrack(mcTrackMap,mcTrack);
00883     if (!rcTrack) continue;
00884     if (rcTrack->flag()<0) continue;
00885 
00886     //-- require primary track
00887     rcTrack = rcTrack->node()->track(primary);
00888     if (!rcTrack) continue;
00889 
00890     nRecon++;
00891 
00892     //-- Fill histograms for the TPC-reconstructed tracks
00893     //             (select the approriate trackgometry)
00894 
00895     //StThreeVectorF rcP = rcTrack->geometry()->momentum();
00896     StThreeVectorF rcP = trackGeometry(rcTrack)->momentum();
00897 
00898     if      (mcTrack->particleDefinition()->pdgEncoding()==   11) hRcElectron->Fill(rcP.perp(),rcP.pseudoRapidity());
00899     else if (mcTrack->particleDefinition()->pdgEncoding()==  211) hRcPionPlus->Fill(rcP.perp(),rcP.pseudoRapidity());
00900     else if (mcTrack->particleDefinition()->pdgEncoding()== -211) hRcPionMin->Fill(rcP.perp(),rcP.pseudoRapidity());
00901     else if (mcTrack->particleDefinition()->pdgEncoding()==  321) hRcKaonPlus->Fill(rcP.perp(),rcP.pseudoRapidity());
00902     else if (mcTrack->particleDefinition()->pdgEncoding()== -321) hRcKaonMin->Fill(rcP.perp(),rcP.pseudoRapidity());
00903     else if (mcTrack->particleDefinition()->pdgEncoding()== 2212) hRcProton->Fill(rcP.perp(),rcP.pseudoRapidity());
00904     else if (mcTrack->particleDefinition()->pdgEncoding()==-2212) hRcAntiProton->Fill(rcP.perp(),rcP.pseudoRapidity());
00905     else LOG_INFO << "cannot fill rc histo" << endm;
00906     float diff = (mcP.mag()-rcP.mag()) / mcP.mag();           
00907     if (mcTrack->particleDefinition()->pdgEncoding()==-211) hMomResPtPion->Fill(mcP.perp(),diff);
00908 
00909 
00910 
00911 
00912     //-- Loop over matched slats from Stage 1 (allMatchedSlatsVec) and look for similar trackKey
00913     //
00914 
00915     if (Debug()) LOG_INFO  << "RC/MC track Id="<< rcTrack->key() << "/" << mcTrack->key() << endm;
00916 
00917     //    Note on trackId vs. StTrack->key():
00918     //      these are *not* the same. trackId (based on the StEvent->trackNodes
00919     //      array index) is used internally to quickly locate tracks in the trackNodes
00920     //      array. StTrack->key(), however, should *always* be used to identify tracks
00921     //      outside the scope of this local code.
00922     //      To get the proper Key() from the trackId use this:
00923     //      trackKey = StEvent->trackNodes[trackId]->track()->key();
00924 
00925     for (tofSlatHitVectorIter matchedSlat=allMatchedSlatsVec.begin();
00926          matchedSlat != allMatchedSlatsVec.end(); matchedSlat++){
00927 
00928       //- All matches from Stage I should have one and only one associated track
00929       if (matchedSlat->trackIdVec.size()!=1){
00930         gMessMgr->Warning() << "Slat with " << matchedSlat->trackIdVec.size() 
00931                             << "track matches? Should not happen!" << endm;
00932         continue;
00933       }
00934 
00935       //- Map trackId to trackKey
00936       unsigned int thisTrackId = matchedSlat->trackIdVec[0];
00937       unsigned int thisKey = nodes[thisTrackId]->track(global)->key();
00938 
00939       if (Debug()) LOG_INFO << " match trackKey (trackId): " << thisKey 
00940                         << " (" << thisTrackId << ")";
00941       
00942       if (thisKey == rcTrack->key()){
00943         if (Debug()) LOG_INFO << "RC-MC-TOFp match";
00944 
00945         nMatch++;
00946 
00947         //- Fill histograms for the reconstructed and matched tracks
00948         if      (mcTrack->particleDefinition()->pdgEncoding()==   11) hMatchElectron->Fill(rcP.perp(),rcP.pseudoRapidity());
00949         else if (mcTrack->particleDefinition()->pdgEncoding()==  211) hMatchPionPlus->Fill(rcP.perp(),rcP.pseudoRapidity());
00950         else if (mcTrack->particleDefinition()->pdgEncoding()== -211) hMatchPionMin->Fill(rcP.perp(),rcP.pseudoRapidity());
00951         else if (mcTrack->particleDefinition()->pdgEncoding()==  321) hMatchKaonPlus->Fill(rcP.perp(),rcP.pseudoRapidity());
00952         else if (mcTrack->particleDefinition()->pdgEncoding()== -321) hMatchKaonMin->Fill(rcP.perp(),rcP.pseudoRapidity());
00953         else if (mcTrack->particleDefinition()->pdgEncoding()== 2212) hMatchProton->Fill(rcP.perp(),rcP.pseudoRapidity());
00954         else if (mcTrack->particleDefinition()->pdgEncoding()==-2212) hMatchAntiProton->Fill(rcP.perp(),rcP.pseudoRapidity());
00955         else LOG_INFO << "cannot fill match histo" << endm;
00956       }
00957       if (Debug()) LOG_INFO << endm;
00958 
00959     }
00960 
00961 
00962 
00963     // Alternative check ... reextrapolate the track and compare slatIds
00964 
00965 //    // A. make sure the track makes it to the ctb
00966 //    StPhysicalHelixD thisHelix = rcTrack->geometry()->helix();
00967 //    //pairD pairL = thisHelix.pathLength(mTofGeom->tofParam().r + mTofGeom->tofParam().counter_thickness);
00968 //    //if(fabs(pairL.first)>=500.0 && fabs(pairL.second)>=500.0) continue;
00969 //          
00970 //    // B. only consider slat extrapolations (no Nslat or Ntrack/slat cuts)
00971 //    tofSlatHitVector hitVec = mTofGeom->tofHelixToArray(thisHelix, validSlatIdVec);
00972 //    if (hitVec.size()==0) continue;
00973 //
00974 //    {
00975 //      LOG_INFO << "TOF HitVec ("<< hitVec.size() << ") slatIn reads:";
00976 //      for (tofSlatHitVectorIter ii=hitVec.begin();ii!=hitVec.end();ii++)
00977 //        LOG_INFO << " " << ii->slatIndex;
00978 //      LOG_INFO << endm;
00979 //    }
00980 //
00981 //
00982 //
00983 //    // C  i. loop over AllHits and select tracks w/  one or more single hit slats
00984 //    //   ii. select on hit profile
00985 //    bool singleHit(false),profiledSingleHit(false);
00986 //    tofSlatHitVectorIter thisHit=hitVec.begin();
00987 //    while (thisHit!=hitVec.end()){
00988 //      tofSlatHitVectorIter allHits = allMatchedSlatsVec.begin();
00989 //      while (allHits!=allMatchedSlatsVec.end()){
00990 //        if (thisHit->slatIndex == allHits->slatIndex){
00991 //          LOG_INFO << "TOF found it in allHits "
00992 //               <<  allHits->trackIdVec.size() << endm;
00993 //          if (allHits->trackIdVec.size()==1){
00994 //            singleHit=true;
00995 //            if (thisHit->hitProfile == 31) profiledSingleHit=true;
00996 //          }
00997 //        }
00998 //        allHits++;
00999 //      }
01000 //      thisHit++;
01001 //    }
01002 //    if (!singleHit) {LOG_INFO << "TOF multiple hit"<<endm;continue;}
01003 //    LOG_INFO << "TOF single hit"<< endm;
01004 //    //if (!profiledSingleHit) continue;
01005 //    //LOG_INFO << "TOF profiled single hit"<< endm;
01006 
01007           
01008   }
01009   gMessMgr->Info("","OST") << "Embedding: #mcTrack/#recon/#match: " <<  mcTracks.size()
01010                            <<"/"<< nRecon << "/" << nMatch << endm;
01011 
01012   return kStOK;
01013 }
01014 
01015 
01016 //---------------------------------------------------------------------------
01017 bool StTofpMcAnalysisMaker::validTrack(StTrack *track){
01018   // 1. no track, no go.
01019   if (!track) return false;
01020 
01021   // 2. track quality flag, should be >0
01022   if (track->flag()<=0) return false;
01023 
01024   // 3. minimum #hits per track
01025   if (track->topologyMap().numberOfHits(kTpcId) < mMinHitsPerTrack) return false;
01026 
01027   // 4. rlast (not implemented, yet)
01028 
01029   return true;
01030 }
01031 
01032 
01033 bool StTofpMcAnalysisMaker::validTofTrack(StTrack *track){
01034   // select valid tracks for time-of-flight calculations
01035 
01036   // 1. track must exist
01037   if (!track) return false;
01038 
01039   // 2. track must be a primary track
01040   if (!dynamic_cast<StPrimaryTrack*>(track)) return false;
01041 
01042   // 3. DCA cut (obsolete?)
01043   double DCA= track->impactParameter();
01044   //int charge = track->geometry()->charge();
01045   //if (mHisto) hTofpDCATrackprimVertex->Fill(DCA*charge);
01046   double mMaxDCA = 999.;
01047   if (DCA > mMaxDCA) {
01048     gMessMgr->Info("","OST") << "dca>max:" << DCA<< endm;
01049     return false;
01050   }
01051 
01052   return true;
01053 }  
01054 
01055 
01056 const StTrackGeometry* StTofpMcAnalysisMaker::trackGeometry(const StTrack* track){
01057   // returns apropriate StTrackGeometry (standard or outerGeometry)
01058   if (!track) return 0;
01059   const StTrackGeometry *thisTrackGeometry;
01060   if (mOuterTrackGeometry)
01061     thisTrackGeometry = track->outerGeometry();
01062   else
01063     thisTrackGeometry = track->geometry();
01064   return thisTrackGeometry;
01065 }
01066 
01067 //--------------------------------------------------------------------------