StRoot: StBTofMatchMaker/StBTofMatchMaker.cxx Source File

00001 /*******************************************************************
00002  *
00003  * $Id: StBTofMatchMaker.cxx,v 1.19 2012/05/07 14:11:16 fisyak Exp $
00004  *
00005  * Author: Xin Dong
00006  *****************************************************************
00007  *
00008  * Description: BTof Match Maker to do the matching between the 
00009  *              fired celles and TPC tracks
00010  *
00011  *****************************************************************
00012  *
00013  * $Log: StBTofMatchMaker.cxx,v $
00014  * Revision 1.19  2012/05/07 14:11:16  fisyak
00015  * Keep btofGeometry in const area for future use
00016  *
00017  * Revision 1.18  2011/08/04 19:14:02  geurts
00018  * Bug fix: allow ideal geometry setting in the case an alignment file is used [Patrick]
00019  *
00020  * Revision 1.17  2011/07/27 16:13:58  geurts
00021  * Alignment calibration modifications [Patrick Huck]:
00022  *  -  modified to open the local Z window cut to determine the z offset
00023  *  -  variables mZLocalCut, mCalculateAlign and mAlignFileName added
00024  *  -  functions setCalculateAlign and setAlignFileName added
00025  *
00026  * Revision 1.16  2010/08/09 19:18:45  geurts
00027  * Include local theta calculation in CellHit structure. Pass LocalTheta info on to TOF PID traits. [Masa]
00028  *
00029  * Revision 1.15  2010/07/16 04:25:16  geurts
00030  * initialize mUseIdealGeometry to be kFALSE in ctor
00031  *
00032  * Revision 1.14  2010/07/14 20:35:21  geurts
00033  * introduce switch to enable ideal MC geometry, without alignment updates. Default: disabled
00034  *
00035  * Revision 1.13  2010/05/25 22:09:38  geurts
00036  * improved database handling and reduced log output
00037  *
00038  * Revision 1.12  2010/05/20 22:58:47  geurts
00039  * Keep BTofMatchMaker from crashing ungracefully when no mEvent or BTOF Collection is found
00040  *
00041  * Revision 1.11  2010/03/22 19:40:28  dongx
00042  * Fixed a bug in setting index2Primary in processMuDst
00043  *
00044  * Revision 1.10  2010/03/19 22:25:39  dongx
00045  * - Added getBTofGeom() function for outside use
00046  * - Remove AddConst(btofGeometry) to avoid crash due to duplication
00047  * - TOT selection window opened to 40 ns
00048  * - Added CPU timer printouts for processStEvent() funciton
00049  *
00050  * Revision 1.9  2010/03/04 21:59:27  dongx
00051  * Further addition in the initial clean up for primary tracks too
00052  *
00053  * Revision 1.8  2010/03/04 21:40:54  dongx
00054  * Added clean up in processMuDst to remove associations done before
00055  *
00056  * Revision 1.7  2010/03/04 03:54:28  dongx
00057  * Further improvement on the processMuDst speed in Step F.
00058  *
00059  * Revision 1.6  2010/03/04 00:08:47  dongx
00060  * Removed primary check for globals at projection in accessing MuDst function as it takes significant CPU time. Waiting for MuDst update
00061  * Added some more debugging output for CPU time usage
00062  *
00063  * Revision 1.5  2010/02/25 05:17:10  dongx
00064  * Geometry initalization moved from Init() to InitRun()
00065  *
00066  * Revision 1.4  2009/09/15 00:30:45  dongx
00067  * 1) Added the functionality to perform the matching with MuDst directly.
00068  * 2) Several updates on the track cuts used for matching
00069  *    - flag<1000 was added
00070  *    - nHits>15 cut was removed
00071  * 3) Created a new StBTofPidTraits for any primary track
00072  * 4) Local Z window cut set to symmetric (fabs(localz)<3.05)
00073  * 5) Some small changes in the LOGGER output.
00074  *
00075  * Revision 1.3  2009/08/26 20:33:56  dongx
00076  * Geometry init moved to Init() function, also allow reading in from others
00077  *
00078  * Revision 1.2  2009/07/24 18:52:53  dongx
00079  * - Local Z window restricted in the projection
00080  * - ToT selection is used firstly when more than one hits associated with a track
00081  * - matchFlag updated
00082  *    0:   no matching
00083  *    1:   1-1 matching
00084  *    2:   1-2 matching, pick up the one with higher ToT value (<25ns)
00085  *    3:   1-2 matching, pick up the one with closest projection posision along y
00086  *
00087  * Revision 1.1  2009/06/23 13:15:03  geurts
00088  * *** empty log message ***
00089  *
00090  *
00091  *******************************************************************/
00092 #include <iostream>
00093 #include "StEventTypes.h"
00094 #include "Stypes.h"
00095 #include "StThreeVectorD.hh"
00096 #include "StThreeVectorF.hh"
00097 #include "StHelix.hh"
00098 #include "StTrackGeometry.h"
00099 #include "StDcaGeometry.h"
00100 #include "StDedxPidTraits.h"
00101 #include "StTrackPidTraits.h"
00102 #include "StBTofPidTraits.h"
00103 #include "StarClassLibrary/StParticleTypes.hh"
00104 #include "StarClassLibrary/StParticleDefinition.hh"
00105 #include "StTpcDedxPidAlgorithm.h"
00106 #include "StEventUtilities/StuRefMult.hh"
00107 #include "PhysicalConstants.h"
00108 #include "StPhysicalHelixD.hh"
00109 #include "StHelix.hh"
00110 #include "StBTofCollection.h"
00111 #include "StBTofUtil/tofPathLength.hh"
00112 #include "StBTofUtil/StBTofGeometry.h"
00113 #include "StBTofUtil/StBTofHitCollection.h"
00114 #include "tables/St_tofConfig_Table.h"
00115 #include "tables/St_tofTrayConfig_Table.h"
00116 #include "TH1.h"
00117 #include "TH2.h"
00118 #include "TFile.h"
00119 #include "TTree.h"
00120 #include "StMessMgr.h"
00121 #include "StMemoryInfo.hh"
00122 #include "StTimer.hh"
00123 
00124 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00125 #include "StMuDSTMaker/COMMON/StMuDst.h"
00126 #include "StMuDSTMaker/COMMON/StMuBTofHit.h"
00127 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00128 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
00129 #include "StMuDSTMaker/COMMON/StMuBTofPidTraits.h"
00130 
00131 #include "StBTofMatchMaker.h"
00132 //#include "TMemStat.h"
00133 
00134 ClassImp(StBTofMatchMaker)
00135 
00136 //---------------------------------------------------------------------------
00137 StBTofMatchMaker::StBTofMatchMaker(const Char_t *name): StMaker(name){
00138   // set default values
00139   mEventCounter = 0;
00140   mAcceptedEventCounter = 0;
00141   mTofEventCounter = 0;
00142   mAcceptAndBeam = 0;
00143 
00144   mBTofGeom = 0;
00145 
00146   mWidthPad = 3.45;
00147   mZLocalCut = 3.05;
00148 
00149   setOuterTrackGeometry();
00150   setMinHitsPerTrack(15);
00151   setMinFitPointsPerTrack(15);
00152   setMinFitPointsOverMax(0.52);
00153   setMaxDCA(9999.);
00154 
00155   setCreateHistoFlag(kFALSE);
00156   setHistoFileName("tofana.root");
00157   setCreateTreeFlag(kFALSE);
00158   setSaveGeometry(kFALSE);
00159   mInitFromOther = kFALSE;
00160   mUseIdealGeometry = kFALSE;
00161   mCalculateAlign   = kFALSE;
00162   setAlignFileName("");
00163   doPrintMemoryInfo = kFALSE;
00164   doPrintCpuInfo    = kFALSE;
00165 
00166   mEvent = 0;
00167   mMuDst = 0;
00168   mMuDstIn = kFALSE;
00169 }
00170 
00171 StBTofMatchMaker::~StBTofMatchMaker(){ /* nope */}
00172 
00173 //void StBTofMatchMaker::Clear(Option_t *opt){StMaker::Clear();}
00174 
00175 //---------------------------------------------------------------------------
00176 Int_t StBTofMatchMaker::Init(){
00177   LOG_INFO << "Initializing match settings:" << endm;
00178   LOG_INFO << "  Minimum fitpoints per track: " << mMinFitPointsPerTrack << endm;
00179   LOG_INFO << "  Maximum DCA: " << mMaxDCA << endm;
00180   
00181   if (!mOuterTrackGeometry) {
00182     LOG_WARN << "using standard trackgeometry()" << endm;
00183   }
00184   
00185   // m_Mode can be set by SetMode() method
00186   if(m_Mode) {
00187 //    setHistoFileName("tofana.root");
00188   } else {
00189     setHistoFileName("");
00190   }
00191 
00192   if (mHisto){
00193     bookHistograms();
00194     LOG_INFO << "Histograms are booked" << endm;
00195     if (mHistoFileName!="") {
00196       LOG_INFO << "Histograms will be stored in " << mHistoFileName.c_str() << endm;
00197     }
00198   }
00199 
00200   // reset event counters
00201   mEventCounter = 0;
00202   mAcceptedEventCounter = 0;
00203   mTofEventCounter = 0;
00204   mAcceptAndBeam = 0;
00205 
00206   // for alignment calculate, we start from ideal geometry
00207   if(mCalculateAlign) {
00208     if (mAlignFileName=="") mUseIdealGeometry = kTRUE;
00209     mZLocalCut = 5.0;
00210   }
00211 
00212   return kStOK;
00213 }
00214 
00215 
00216 //---------------------------------------------------------------------------
00217 Int_t StBTofMatchMaker::InitRun(Int_t runnumber){
00218 
00219   // determine TOF configuration from run#
00220   LOG_INFO << "Initializing BTOF Geometry:" << endm;
00222   // TOF geometry initializtion -- from GEANT geometry directly
00223   //                               need St_geant_Maker be loaded before
00225   StTimer geomTimer;
00226   geomTimer.start();
00227   mBTofGeom = 0;
00228   if(TDataSet *geom = GetDataSet("btofGeometry")) {
00229     mBTofGeom = (StBTofGeometry *)geom->GetObject();
00230     LOG_INFO << " Found btofGeometry ... " << endm; 
00231     mInitFromOther = kTRUE;
00232   } else {
00233     mBTofGeom = new StBTofGeometry("btofGeometry","btofGeometry in MatchMaker");
00234     LOG_INFO << " Create a new btofGeometry ... " << endm;
00235     AddConst(new TObjectSet("btofGeometry",mBTofGeom));
00236   } 
00237   if(mBTofGeom && !mBTofGeom->IsInitDone()) {
00238     LOG_INFO << " BTofGeometry initialization ... " << endm;
00239     //fg if(runnumber<1000000) mBTofGeom->SetMCOn();
00240     if (mUseIdealGeometry) mBTofGeom->SetMCOn();
00241     else                   mBTofGeom->SetMCOff();
00242     LOG_INFO << " Alignment file: " << mAlignFileName.c_str() << endm;
00243     mBTofGeom->SetAlignFile(mAlignFileName.c_str());
00244       TVolume *starHall = (TVolume *)GetDataSet("HALL");
00245       mBTofGeom->Init(this, starHall);
00246   }
00247 
00248   geomTimer.stop();
00249   LOG_INFO << "CPU time for StBTofGeometry initialization : " << geomTimer.elapsedTime() << " sec" << endm;
00250 
00251   return kStOK;
00252 }
00253 
00254 //----------------------------------------------------------------------------
00255 Int_t StBTofMatchMaker::FinishRun(Int_t runnumber){
00256 
00257   LOG_INFO << "StBTofMatchMaker -- cleaning up geometry (Finish)" << endm;
00258   if (!mInitFromOther && mBTofGeom) delete mBTofGeom;
00259   mBTofGeom=0;
00260 
00261   return kStOK;
00262 }
00263 
00264 
00265 //---------------------------------------------------------------------------
00266 Int_t StBTofMatchMaker::Finish(){
00267 
00268   LOG_DEBUG << "StBTofMatchMaker -----  RUN SUMMARY ----- (Finish)\n"
00269        << "\tProcessed "  << mEventCounter << " events."
00270        << " Accepted  "   << mAcceptedEventCounter << " events."
00271        << " Rejected  "   << mEventCounter - mAcceptedEventCounter << " events\n"
00272        << "\tTOF events " << mTofEventCounter
00273        << "\t Accept & Beam   "   << mAcceptAndBeam   << " events" << endm;
00274   
00275   //if (mHisto) writeHistogramsToFile();
00276   if (mHistoFileName!="") writeHistogramsToFile();
00277 
00278 
00279   return kStOK;
00280 }
00281 
00282 //---------------------------------------------------------------------------
00283 Int_t StBTofMatchMaker::Make(){
00284   LOG_INFO << "StBTofMatchMaker -- welcome" << endm;
00285 
00286   if(mMuDstIn) processMuDst();
00287   else         processStEvent();
00288 
00289   return kStOK;
00290 }
00291 
00292 //---------------------------------------------------------------------------
00293 void StBTofMatchMaker::processStEvent(){
00294 
00295   if(mHisto) mEventCounterHisto->Fill(0);
00296   // event selection ...
00297   mEvent = (StEvent *) GetInputDS("StEvent");
00298   if(!mEvent || !(mEvent->btofCollection()) || !(mEvent->btofCollection()->hitsPresent()) ) {
00299     if (!mEvent) {LOG_INFO << "no StEvent" << endm;}
00300     else
00301       if (!(mEvent->btofCollection())) {LOG_INFO << "no BTof Collection" << endm;}
00302     else
00303       if (!(mEvent->btofCollection()->hitsPresent()) ) LOG_INFO << "no BTOF hits present" << endm;
00304     LOG_INFO << "StBTofMatchMaker -- nothing to do ... bye-bye" << endm;
00305     return;
00306   }
00307   if(mHisto) mEventCounterHisto->Fill(1);
00308 
00309   // timing & memory info -only when requested-
00310   StTimer timer;
00311   if (doPrintCpuInfo) timer.start();
00312   if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
00313 
00314   //.........................................................................
00315   // check for tofCollection and fill local copy with ADC and TDC data
00316   StBTofCollection *theTof = mEvent->btofCollection();
00317 
00318   //.........................................................................
00319   // read data from StBTofHit
00321   //
00322   tofCellHitVector daqCellsHitVec;
00323   idVector validModuleVec;
00324 
00325   // multi-tray system
00326   StSPtrVecBTofHit& tofHits = theTof->tofHits();
00327   LOG_INFO << " Number of BTOF Hits = " << tofHits.size() << endm;
00328   for(size_t i=0;i<tofHits.size();i++) {
00329     //fg StBTofHit* aHit = dynamic_cast<StBTofHit *>(tofHits[i]);
00330     StBTofHit* aHit = tofHits[i];
00331      if(!aHit) continue;
00332      if(aHit->tray()<=0||aHit->tray()>mNTray) continue;  // barrel tray hits
00333 
00334      int trayId = aHit->tray();
00335      int moduleId = aHit->module();
00336      int cellId = aHit->cell();
00337 
00338      if(Debug()) { LOG_INFO <<"A: fired hit in " << " tray="<< trayId << " module="<<moduleId<<" cell="<<cellId<<endm; }
00339 
00340      StructCellHit aDaqCellHit;
00341      aDaqCellHit.tray = trayId;
00342      aDaqCellHit.module = moduleId;
00343      aDaqCellHit.cell = cellId;
00344      aDaqCellHit.tot = aHit->tot();
00345      aDaqCellHit.index2BTofHit = i;
00346      daqCellsHitVec.push_back(aDaqCellHit);
00347 
00348      // additional valid number configuration
00349      int id = trayId*100+moduleId;
00350      //fg bool ifind = kFALSE;
00351      //fg for(size_t im=0;im<validModuleVec.size();im++) {
00352      //fg   if(id==validModuleVec[im]) {
00353      //fg     ifind = kTRUE;
00354      //fg     break;
00355      //fg   }
00356      //fg }
00357      //fg if(!ifind) validModuleVec.push_back(id);
00358      if (find(validModuleVec.begin(), validModuleVec.end(), id) == validModuleVec.end())
00359        validModuleVec.push_back(id);
00360 
00361      if(mHisto) {
00362        mDaqOccupancy[trayId-1]->Fill((moduleId-1)*mNCell+(cellId-1));
00363      }      
00364   }
00365       
00366   // end of Sect.A
00367   if(Debug()) { 
00368     LOG_INFO << "    total # of cells = " << daqCellsHitVec.size() << endm;
00369     for(size_t iv = 0;iv<validModuleVec.size();iv++) {
00370       LOG_DEBUG << " module # " << validModuleVec[iv] << " Valid! " << endm;
00371     }
00372   }
00373   if(mHisto) {
00374     mCellsMultInEvent->Fill(daqCellsHitVec.size());
00375     if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
00376   }
00377 //  if(!daqCellsHitVec.size()) return;
00378   if (doPrintCpuInfo) {
00379     timer.stop();
00380     LOG_INFO << "CPU time after Step A - loading hits : "
00381          << timer.elapsedTime() << " sec" << endm;
00382     timer.start();
00383   }
00384 
00385   //.........................................................................
00387   //
00388   tofCellHitVector allCellsHitVec;
00389   StructCellHit cellHit;
00390 
00391   StSPtrVecTrackNode& nodes = mEvent->trackNodes();
00392   Int_t nAllTracks=0;
00393   Int_t nPrimaryHits = 0;
00394   for (unsigned int iNode=0; iNode<nodes.size(); iNode++){
00395     tofCellHitVector cellHitVec;
00396     //    cellHitVec.clear();
00397     StGlobalTrack *theTrack = dynamic_cast<StGlobalTrack*>(nodes[iNode]->track(global));
00398     if(!theTrack) continue;
00399 
00400     bool isPrimary = kFALSE;
00401     StPrimaryTrack *pTrack = dynamic_cast<StPrimaryTrack*>(theTrack->node()->track(primary));
00402     if(pTrack) isPrimary = kTRUE;
00403 
00404     StThreeVectorF mom = theTrack->geometry()->momentum();
00405     float pt = mom.perp();
00406     float eta = mom.pseudoRapidity();
00407     float phi = mom.phi();
00408     //fg if (phi<0.) phi += 2.*3.14159;
00409     if (phi<0.) phi += 2.*M_PI;
00410 
00411     float dEdx = -999.;
00412     int ndEdxpts = 0;
00413     float nSigmaPion = -999.;
00414     static StTpcDedxPidAlgorithm PidAlgorithm;
00415     static StPionPlus* Pion = StPionPlus::instance();
00416 
00417     const StParticleDefinition* pd = theTrack->pidTraits(PidAlgorithm);
00418     if (pd) {
00419       nSigmaPion = PidAlgorithm.numberOfSigma(Pion);
00420     }
00421     if ( PidAlgorithm.traits() ) {
00422       dEdx = PidAlgorithm.traits()->mean();
00423       ndEdxpts = PidAlgorithm.traits()->numberOfPoints();
00424     }
00425     int nfitpts = theTrack->fitTraits().numberOfFitPoints(kTpcId);
00426 
00427     // make sure we have a track, a miniDST might have removed it...
00428     if (validTrack(theTrack)){
00429       if(mHisto) {
00430         mTrackPtEta->Fill(pt, eta);
00431         mTrackPtPhi->Fill(pt, phi);
00432         mTrackNFitPts->Fill(nfitpts);
00433         if(dEdx>0.) mTrackdEdxvsp->Fill(mom.mag(), dEdx*1.e6);
00434         if(fabs(nSigmaPion)<5.) mNSigmaPivsPt->Fill(pt, nSigmaPion+5.*theTrack->geometry()->charge());
00435       }
00436 
00437       if(mSaveTree) {
00438         trackTree.pt = pt;
00439         trackTree.eta = eta;
00440         trackTree.phi = phi;
00441         trackTree.nfitpts = nfitpts;
00442         trackTree.dEdx = dEdx*1.e6;
00443         trackTree.ndEdxpts = ndEdxpts;
00444         trackTree.charge = theTrack->geometry()->charge();      
00445         trackTree.projTrayId = 0;
00446         trackTree.projCellChan = -1;
00447         trackTree.projY = -999.;
00448         trackTree.projZ = -999.;
00449       }
00450 
00451       nAllTracks++;
00452       StPhysicalHelixD theHelix = trackGeometry(theTrack)->helix();
00453 
00454 //      IntVec projTrayVec;
00455 //      if(!mBTofGeom->projTrayVector(theHelix, projTrayVec)) continue;
00456 
00457       IntVec idVec;
00458       DoubleVec pathVec, thetaVec;
00459       PointVec  crossVec;
00460 
00461 //       idVec.clear();
00462 //       pathVec.clear();
00463 //       crossVec.clear();
00464 
00465       Int_t ncells = 0;
00466       if(mBTofGeom->HelixCrossCellIds(theHelix,idVec,pathVec,crossVec,thetaVec) ) {
00467 //      if(mBTofGeom->HelixCrossCellIds(theHelix, validModuleVec, projTrayVec, idVec, pathVec, crossVec)) {
00468         Int_t cells = idVec.size();
00469         for (Int_t i=0; i<cells; i++) {
00470             Int_t icell,imodule,itray;
00471             Double_t local[3],global[3];
00472             for(Int_t i2=0;i2<3;i2++){
00473                  local[i2]=0;
00474             }
00475             global[0]=crossVec[i].x();
00476             global[1]=crossVec[i].y();
00477             global[2]=crossVec[i].z();
00478             mBTofGeom->DecodeCellId(idVec[i], icell, imodule, itray);
00479 //          LOG_INFO << " decode " << idVec[i] << "  to tray#" << itray << " module#" << imodule << " cell#" << icell << endm;
00480             StBTofGeomSensor* sensor = mBTofGeom->GetGeomSensor(imodule,itray);
00481             if(!sensor) {
00482               LOG_WARN << " No sensitive module in the projection??? -- Something weird!!! " << endm;
00483               continue;
00484             }
00485             sensor->Master2Local(&global[0],&local[0]);
00486             icell = sensor->FindCellIndex(local);
00487             //      StThreeVectorD glo=sensor->GetCenterPosition();
00488             StThreeVectorD glo(global[0], global[1], global[2]);
00489             StThreeVectorD hitPos(local[0], local[1], local[2]);
00490 //          delete sensor;   /// function in StBTofGeometry modified. Don't delete here.
00491             if (fabs(local[2])<mZLocalCut) {   // to be consistent with GEANT geometry
00492                                          // and the alignment calibration
00493 //          if (local[2]<=2.7&&local[2]>=-3.4) {
00494 //            if (fabs(local[2])<4.) {   // loose local z cut at first step
00495               ncells++;
00496               cellHit.tray = itray;
00497               cellHit.module = imodule;
00498               cellHit.cell = icell;
00499               cellHit.trackIdVec.push_back((Int_t)iNode);
00500               cellHit.hitPosition = glo;        // global position
00501               cellHit.zhit = (Float_t)hitPos.z();
00502               cellHit.yhit = (Float_t)hitPos.y();
00503               cellHit.theta = (Double_t)thetaVec[i];
00504               cellHitVec.push_back(cellHit);
00505               allCellsHitVec.push_back(cellHit);
00506 
00507               if(isPrimary) nPrimaryHits++;
00508 
00509               if(mHisto) {
00510                 mDaqOccupancyProj[itray-1]->Fill((imodule-1)*mNCell+(icell-1));
00511                 mHitsPosition->Fill(hitPos.y(), hitPos.z());
00512               }
00513               
00514               if(mSaveTree) {
00515                 trackTree.projTrayId = itray;
00516                 trackTree.projCellChan = (imodule-1)*mNCell+(icell-1);
00517                 trackTree.projY = local[1];
00518                 trackTree.projZ = local[2];
00519               }
00520 
00521               if(Debug()) {
00522                 LOG_DEBUG <<"B: nodeid=" << iNode << "  projected in " << " tray="<< itray << " module="<<imodule<<" cell="<<icell<<endm;
00523                 LOG_DEBUG <<"   hit position " << hitPos << endm;
00524                 // LOG_DEBUG <<"   momemtum= " << pt << " " << eta << " " << phi << endm;
00525               }
00526             }
00527         } // for (Int_t i=0...)
00528       } // endif(helixcross...)
00529       if(ncells>0&&mHisto) mHitsMultPerTrack->Fill(ncells);
00530 
00531       if(mHisto && mSaveTree) mTrackTree->Fill();
00532 
00533     } // if(ValidTrack).. 
00534   } // loop over nodes
00535   if(Debug()) { LOG_INFO << "B:  matched/available/total #tracknodes: " <<allCellsHitVec.size() << "/" <<nAllTracks << "/" << nodes.size() << endm; }
00536   if(mHisto) {
00537     mHitsMultInEvent->Fill(allCellsHitVec.size());
00538     mHitsPrimaryInEvent->Fill(nPrimaryHits);
00539     if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
00540   }
00541   // end of Sect.B
00542   if (doPrintCpuInfo) {
00543     timer.stop();
00544     LOG_INFO << "CPU time after Step B - projection : "
00545          << timer.elapsedTime() << " sec" << endm;
00546     timer.start();
00547   }
00548 
00549   //.........................................................................
00551   //
00552   tofCellHitVector matchHitCellsVec;
00553 
00554   tofCellHitVectorIter daqIter = daqCellsHitVec.begin();
00555   for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
00556     tofCellHitVectorIter proIter = allCellsHitVec.begin();
00557     for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
00558 
00559       int daqIndex = (daqIter->module-1)*6 + (daqIter->cell-1);
00560       int proIndex = (proIter->module-1)*6 + (proIter->cell-1);
00561       int hisIndex = daqIter->tray - 1;
00562 //      int daqAllIndex = (daqIter->tray - 1)*mNTOF + daqIndex;
00563 //      int proAllIndex = (proIter->tray - 1)*mNTOF + proIndex;
00564       if(daqIter->tray==proIter->tray) {
00565         if (mHisto) {
00566           if(hisIndex>=0&&hisIndex<mNTray) {
00567             mHitCorr[hisIndex]->Fill(proIndex,daqIndex);
00568             mHitCorrModule[hisIndex]->Fill(proIter->module-1,daqIter->module-1);
00569           } else {
00570             LOG_WARN << " weird tray # " << daqIter->tray << endm;
00571           }
00572         }
00573       }
00574       if( (daqIter->tray==proIter->tray)&& 
00575           (daqIter->module==proIter->module) &&
00576           ( ( (proIter->cell==6)&&((proIter->cell==daqIter->cell) ||
00577                                    (proIter->cell==daqIter->cell+1)) )
00578             || ( (proIter->cell==1)&&((proIter->cell==daqIter->cell) ||
00579                                       (proIter->cell==daqIter->cell-1)) )
00580             || ( (proIter->cell>=2&&proIter->cell<=6) &&
00581                  ( (proIter->cell==daqIter->cell) ||
00582                    (proIter->cell==daqIter->cell-1) ||
00583                    (proIter->cell==daqIter->cell+1) ) ) ) ) {
00584         cellHit.tray = daqIter->tray;
00585         cellHit.module = daqIter->module;
00586         cellHit.cell = daqIter->cell;
00587         cellHit.hitPosition = proIter->hitPosition;
00588         cellHit.trackIdVec = proIter->trackIdVec;
00589         cellHit.zhit = proIter->zhit;
00590         cellHit.yhit = proIter->yhit;
00591         cellHit.tot = daqIter->tot;
00592         cellHit.index2BTofHit = daqIter->index2BTofHit;
00593         cellHit.theta = proIter->theta;
00594         matchHitCellsVec.push_back(cellHit);
00595       }
00596     }
00597   } //end {sec. C}
00598   if(Debug()) { LOG_INFO << "C: before/after: " << allCellsHitVec.size() << "/" << matchHitCellsVec.size() << endm; }
00599   if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
00600   if (doPrintCpuInfo) {
00601     timer.stop();
00602     LOG_INFO << "CPU time after Step C - matching : "
00603          << timer.elapsedTime() << " sec" << endm;
00604     timer.start();
00605   }
00606 
00607   //.........................................................................
00609   //
00610   Int_t nSingleHitCells(0);
00611   Int_t nMultiHitsCells(0);
00612 
00613   tofCellHitVector singleHitCellsVec;
00614   tofCellHitVector multiHitsCellsVec;
00615 
00616   tofCellHitVector tempVec = matchHitCellsVec;
00617   tofCellHitVector erasedVec = tempVec;
00618   while (tempVec.size() != 0) {
00619     Int_t nTracks = 0;
00620     idVector trackIdVec;
00621 
00622     tofCellHitVectorIter tempIter=tempVec.begin();
00623     tofCellHitVectorIter erasedIter=erasedVec.begin();
00624     while(erasedIter!= erasedVec.end()) {
00625       if(tempIter->tray == erasedIter->tray &&
00626          tempIter->module == erasedIter->module &&
00627          tempIter->cell == erasedIter->cell) {
00628         nTracks++;
00629         trackIdVec.push_back(erasedIter->trackIdVec.back());  // merge
00630         erasedVec.erase(erasedIter);
00631         erasedIter--;
00632       }
00633       erasedIter++;
00634     }
00635 
00636     cellHit.cell = tempIter->cell;
00637     cellHit.module = tempIter->module;
00638     cellHit.tray = tempIter->tray;
00639     cellHit.hitPosition = tempIter->hitPosition;
00640     cellHit.trackIdVec = trackIdVec;
00641     cellHit.zhit = tempIter->zhit;
00642     cellHit.yhit = tempIter->yhit;
00643     cellHit.tot = tempIter->tot;
00644     cellHit.index2BTofHit = tempIter->index2BTofHit;
00645     cellHit.theta = tempIter->theta;
00646 
00647     Float_t ycenter = (tempIter->cell-1-2.5)*mWidthPad;
00648     Float_t dy = tempIter->yhit - ycenter;
00649     Float_t dz = tempIter->zhit;
00650 
00651     if(mHisto) {
00652       mTracksPerCellMatch1->Fill(trackIdVec.size());
00653 //      mDaqOccupancyMatch1->Fill((tempIter->module-1)*mNCell+(tempIter->cell-1));
00654       mDeltaHitMatch1->Fill(dy, dz);
00655     }
00656 
00657     if (nTracks==1){
00658       nSingleHitCells++;      
00659       singleHitCellsVec.push_back(cellHit);
00660     } else if (nTracks>1){
00661       nMultiHitsCells++;
00662       multiHitsCellsVec.push_back(cellHit);
00663       // for multiple hit cells either discard (yes) or
00664       // find the most likely candidate.
00665     } else {
00666       LOG_WARN << "D: no tracks extrapolate to matched cell ... should not happen!" << endm;
00667     }
00668     
00669     if(Debug()) { 
00670       LOG_DEBUG << "D: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
00671       idVectorIter ij=trackIdVec.begin();
00672       while (ij != trackIdVec.end()) { LOG_DEBUG << " " << *ij; ij++; }
00673       LOG_DEBUG << endm;
00674     }
00675 
00676     tempVec = erasedVec;
00677   }
00678   if(Debug()) { LOG_INFO << "D: before/after: " << matchHitCellsVec.size() << "/" << singleHitCellsVec.size() << endm; }
00679   //end of Sect.D
00680 
00681   if(mHisto) {
00682     mCellsPerEventMatch1->Fill(singleHitCellsVec.size()+multiHitsCellsVec.size());
00683     if(singleHitCellsVec.size()) mEventCounterHisto->Fill(9);
00684   } 
00685   if (doPrintCpuInfo) {
00686     timer.stop();
00687     LOG_INFO << "CPU time after Step D - erasing : "
00688          << timer.elapsedTime() << " sec" << endm;
00689     timer.start();
00690   }
00691 
00692   //.........................................................................
00694   //
00695   tofCellHitVector FinalMatchedCellsVec;
00696   //  FinalMatchedCellsVec.clear();
00697   tempVec = singleHitCellsVec;
00698   if(mHisto) {
00699     mCellsPerEventMatch2->Fill(tempVec.size());
00700     for(unsigned int ii=0;ii<tempVec.size();ii++) {
00701       mTracksPerCellMatch2->Fill(tempVec[ii].trackIdVec.size());
00702 //      mDaqOccupancyMatch2->Fill((tempVec[ii].module-1)*mNCell+(tempVec[ii].cell-1));
00703       Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
00704       Float_t dy = tempVec[ii].yhit-ycenter;
00705       Float_t dz = tempVec[ii].zhit;
00706       mDeltaHitMatch2->Fill(dy, dz);
00707     }
00708   }
00709 
00710   erasedVec = tempVec;
00711   while (tempVec.size() != 0) {
00712     StructCellHit cellHit;
00713     Int_t nCells = 0;
00714     idVector vTrackId;
00715     vector<StThreeVectorD> vPosition;
00716     vector<Int_t> vchannel, vtray, vmodule, vcell;
00717     vector<Float_t> vzhit, vyhit;
00718     vector<Double_t> vtot, vtheta;
00719     vector<Int_t> vindex2BTofHit;
00720 
00721     tofCellHitVectorIter tempIter=tempVec.begin();
00722     tofCellHitVectorIter erasedIter=erasedVec.begin();
00723     while(erasedIter!= erasedVec.end()) {
00724       if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
00725         nCells++;
00726         vtray.push_back(erasedIter->tray);
00727         vmodule.push_back(erasedIter->module);
00728         vcell.push_back(erasedIter->cell);
00729         vPosition.push_back(erasedIter->hitPosition);
00730         vTrackId.push_back(erasedIter->trackIdVec.back());
00731         vzhit.push_back(erasedIter->zhit);
00732         vyhit.push_back(erasedIter->yhit);
00733         vtot.push_back(erasedIter->tot);
00734         vindex2BTofHit.push_back(erasedIter->index2BTofHit);
00735         vtheta.push_back(erasedIter->theta);
00736 
00737         erasedVec.erase(erasedIter);
00738         erasedIter--;
00739       }
00740       erasedIter++;
00741     }
00742 
00743     if (nCells==1){
00744       // for singly hit cell, copy data in singleHitCellsVec
00745       cellHit.tray = vtray[0];
00746       cellHit.module = vmodule[0];
00747       cellHit.cell = vcell[0];
00748       cellHit.trackIdVec.push_back(vTrackId[0]);
00749       cellHit.hitPosition = vPosition[0];
00750       cellHit.matchFlag = 1; 
00751       cellHit.zhit = vzhit[0];
00752       cellHit.yhit = vyhit[0];
00753       cellHit.tot = vtot[0];
00754       cellHit.index2BTofHit = vindex2BTofHit[0];
00755       cellHit.theta = vtheta[0];
00756 
00757       FinalMatchedCellsVec.push_back(cellHit);
00758 
00759       // debugging output
00760       if(Debug()) {
00761         LOG_DEBUG << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
00762         idVectorIter ij=vTrackId.begin();
00763         while (ij != vTrackId.end()) { LOG_DEBUG << " " << *ij; ij++; }
00764         LOG_DEBUG << endm;
00765       }
00766       
00767     }
00768     else if (nCells>1){   // for multiple hit cells  find the most likely candidate.
00769       Int_t thiscandidate(-99);
00770       Int_t thisMatchFlag(0);
00771 
00772       // sort on tot
00773       Float_t tot(0.);
00774       vector<Int_t> ttCandidates;
00775       for (Int_t i=0;i<nCells;i++) {
00776         Double_t tt = vtot[i];
00777         if(tt<40.&&tt>tot) {    // open the ToT cut to 40 ns
00778           tot = tt;
00779           ttCandidates.clear();
00780           ttCandidates.push_back(i);
00781         } else if (tt==tot) {
00782           ttCandidates.push_back(i);
00783         }
00784       }
00785       if (ttCandidates.size()==1) {
00786         thiscandidate = ttCandidates[0];
00787         thisMatchFlag = 2;
00788       } else if (ttCandidates.size()>1) {  // sort on hitposition
00789         Float_t ss(99.);
00790         vector<Int_t> ssCandidates;
00791         for(size_t j=0;j<ttCandidates.size();j++) {
00792           Float_t yy = vyhit[ttCandidates[j]];
00793           Float_t ycell = (vcell[ttCandidates[j]]-1-2.5)*mWidthPad;
00794           Float_t ll = fabs(yy-ycell);
00795           if(ll<ss) {
00796             ss = ll; 
00797             ssCandidates.clear();
00798             ssCandidates.push_back(ttCandidates[j]);
00799           }else if  (ll==ss)
00800             ssCandidates.push_back(ttCandidates[j]);
00801         }
00802         if (ssCandidates.size()==1){
00803           thiscandidate = ssCandidates[0];
00804           thisMatchFlag = 3;
00805         }
00806       }
00807 
00808       if (thiscandidate>=0) {
00809         cellHit.tray = vtray[thiscandidate];
00810         cellHit.module = vmodule[thiscandidate];
00811         cellHit.cell = vcell[thiscandidate];
00812         cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
00813         cellHit.hitPosition = vPosition[thiscandidate];
00814         cellHit.matchFlag = thisMatchFlag;
00815         cellHit.zhit = vzhit[thiscandidate];
00816         cellHit.yhit = vyhit[thiscandidate];
00817         cellHit.tot = vtot[thiscandidate];
00818         cellHit.index2BTofHit = vindex2BTofHit[thiscandidate];
00819         cellHit.theta = vtheta[thiscandidate];
00820 
00821         FinalMatchedCellsVec.push_back(cellHit);
00822         
00823         // debugging output
00824         if(Debug()) { LOG_DEBUG << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:" << vTrackId[thiscandidate] << endm; }
00825       }
00826 
00827     } else {
00828       LOG_WARN << "E: no cells belong to this track ... should not happen!" << endm;
00829     }
00830 
00831     tempVec = erasedVec;
00832   }
00833 
00834   if(Debug()) { LOG_INFO << "E: before/after: " << singleHitCellsVec.size() << "/" << FinalMatchedCellsVec.size() << endm; }
00835   // end of Sect.E
00836   if (doPrintCpuInfo) {
00837     timer.stop();
00838     LOG_INFO << "CPU time after Step E - sorting : "
00839          << timer.elapsedTime() << " sec" << endm;
00840     timer.start();
00841   }
00842 
00843   //.........................................................................
00845   //
00846   if(mHisto) {
00847     if(FinalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
00848     mCellsPerEventMatch3->Fill(FinalMatchedCellsVec.size());
00849   }
00850 
00851 //  StSPtrVecBTofHit& tofHits = theTof->tofHits();
00852   Int_t nValidSingleHitCells(0), nValidSinglePrimHitCells(0);
00853 
00854   for (size_t ii=0; ii < FinalMatchedCellsVec.size(); ii++){
00855     Int_t tray = FinalMatchedCellsVec[ii].tray;
00856     Int_t module = FinalMatchedCellsVec[ii].module;
00857     Int_t cell = FinalMatchedCellsVec[ii].cell;
00858 
00859     Float_t ycenter = (cell-1-2.5)*mWidthPad;
00860     Float_t dy = FinalMatchedCellsVec[ii].yhit - ycenter;
00861     Float_t dz = FinalMatchedCellsVec[ii].zhit;
00862     if (FinalMatchedCellsVec[ii].trackIdVec.size()!=1)
00863       LOG_WARN << "F: WHAT!?!  mult.matched cell in single cell list " << tray << " " << module << " " << cell << endm;
00864 
00865     if(mHisto) {
00866       mTracksPerCellMatch3->Fill(FinalMatchedCellsVec[ii].trackIdVec.size());
00867 //      mDaqOccupancyMatch3->Fill((module-1)*mNCell+(cell-1));
00868       mDeltaHitMatch3->Fill(dy, dz);
00869       mDeltaHitFinal[tray-1]->Fill(dy,dz);      
00870     }
00871 
00872     // get track-id from cell hit vector
00873     int trackNode = FinalMatchedCellsVec[ii].trackIdVec[0];
00874     StGlobalTrack *globalTrack = dynamic_cast<StGlobalTrack*>(nodes[trackNode]->track(global));
00875     if(!globalTrack) {
00876       LOG_WARN << "Wrong global track!" << endm;
00877       continue;
00878     }
00879 
00880     // Fill association in TOF Hit Collection
00881     //fg StBTofHit *tofHit = dynamic_cast<StBTofHit *>(tofHits[FinalMatchedCellsVec[ii].index2BTofHit]);
00882     StBTofHit *tofHit = tofHits[FinalMatchedCellsVec[ii].index2BTofHit];
00883     if(tofHit->tray()!=tray || tofHit->module()!=module || tofHit->cell()!=cell) {
00884       LOG_WARN << "Wrong hit in the BTofHitCollection!" << endm;
00885       continue;
00886     }
00887     nValidSingleHitCells++;
00888 
00890     tofHit->setAssociatedTrack(globalTrack);
00891 
00892     // Fill the matched data in StBTofPidTraits
00893     StBTofPidTraits *pidTof = new StBTofPidTraits();
00894     pidTof->setTofHit(tofHit);
00895     pidTof->setMatchFlag(FinalMatchedCellsVec[ii].matchFlag);
00896     pidTof->setYLocal(dy);
00897     pidTof->setZLocal(FinalMatchedCellsVec[ii].zhit);
00898     pidTof->setPosition(FinalMatchedCellsVec[ii].hitPosition);
00899     pidTof->setThetaLocal(FinalMatchedCellsVec[ii].theta);
00900     globalTrack->addPidTraits(pidTof);
00901 
00902     StPrimaryTrack *pTrack = dynamic_cast<StPrimaryTrack*>(nodes[trackNode]->track(primary));
00903     if(pTrack) {
00904       nValidSinglePrimHitCells++;
00905       StBTofPidTraits *ppidTof = new StBTofPidTraits();
00906       ppidTof->setTofHit(tofHit);
00907       ppidTof->setMatchFlag(FinalMatchedCellsVec[ii].matchFlag);
00908       ppidTof->setYLocal(dy);
00909       ppidTof->setZLocal(FinalMatchedCellsVec[ii].zhit);
00910       ppidTof->setPosition(FinalMatchedCellsVec[ii].hitPosition);
00911       ppidTof->setThetaLocal(FinalMatchedCellsVec[ii].theta);
00912       pTrack->addPidTraits(ppidTof);
00913     }
00914 
00915   } // end final matched cells
00916   
00917   if(mHisto) {
00918     mCellsPrimaryPerEventMatch3->Fill(nValidSinglePrimHitCells);
00919   }
00920   
00921   if(Debug()) { LOG_INFO << "F: before/after" << FinalMatchedCellsVec.size() << "/" <<nValidSinglePrimHitCells << endm; }
00922  // end of Sect.F
00923   if (doPrintCpuInfo) {
00924     timer.stop();
00925     LOG_INFO << "CPU time after Step F - final : "
00926          << timer.elapsedTime() << " sec" << endm;
00927     timer.start();
00928   }
00929 
00930   LOG_INFO << "#(daq): " << daqCellsHitVec.size()
00931        << " #(proj): " << allCellsHitVec.size()
00932        << " #(prim proj): " << nPrimaryHits
00933        << " #(matched): " << FinalMatchedCellsVec.size() 
00934        << " #(single valid): " << nValidSingleHitCells
00935        << " #(single prim valid): " << nValidSinglePrimHitCells
00936        << endm;
00937 
00938   //check StEvent collections --
00939   if (theTof->hitsPresent()){
00940     LOG_DEBUG << " BTofCollection: hit container present."<<endm;
00941     if (Debug()){
00942       StSPtrVecBTofHit& tmpCellTofVec = theTof->tofHits();
00943       LOG_INFO << " # of hits in this event:" << tmpCellTofVec.size() << endm;
00944       for (size_t i = 0; i < tmpCellTofVec.size(); i++) {
00945         StBTofHit* p = tmpCellTofVec[i];
00946         LOG_INFO << (*p) << endm;
00947       }
00948     }
00949   }
00950   //-- end check
00951 
00952 
00953   if (doPrintMemoryInfo) {
00954         StMemoryInfo::instance()->snapshot();
00955         StMemoryInfo::instance()->print();
00956   }
00957   if (doPrintCpuInfo) {
00958     timer.stop();
00959     LOG_INFO << "CPU time for StBTofMatchMaker::Make(): "
00960          << timer.elapsedTime() << " sec" << endm;
00961   }
00962 
00963   LOG_DEBUG << "StBTofMatchMaker -- bye-bye" << endm;
00964 
00965   return;
00966 }
00967 
00968 //---------------------------------------------------------------------------
00969 void StBTofMatchMaker::processMuDst(){
00970 
00971   if(mHisto) mEventCounterHisto->Fill(0);
00972   // event selection ...
00973   StMuDstMaker *mMuDstMaker = (StMuDstMaker *)GetMaker("MuDst");
00974   if(!mMuDstMaker) {
00975     LOG_WARN << " No MuDstMaker ...  bye-bye ... " << endm;
00976     return;
00977   }
00978   mMuDst = mMuDstMaker->muDst();
00979   if(!mMuDst) {
00980     LOG_WARN << " No MuDst ... bye-bye" << endm;
00981     return;
00982   }
00983 
00984   // timing & memory info -only when requested-
00985   StTimer timer;
00986   if (doPrintCpuInfo) timer.start();
00987   if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
00988 
00989   //.........................................................................
00990   // read data from StBTofHit
00992   //
00993   tofCellHitVector daqCellsHitVec;
00994   idVector validModuleVec;
00995 
00996   // multi-tray system
00997   Int_t nhits = mMuDst->numberOfBTofHit();
00998   LOG_INFO << " Number of BTOF Hits = " << nhits << endm;
00999   if(mHisto&&nhits>0) mEventCounterHisto->Fill(1);
01000   for(int i=0;i<nhits;i++) {
01001      StMuBTofHit *aHit = (StMuBTofHit*)mMuDst->btofHit(i);
01002      if(!aHit) continue;
01003      if(aHit->tray()<=0||aHit->tray()>mNTray) continue;  // barrel tray hits
01004 
01005      // clean up any association done before
01006      aHit->setIndex2Primary(-1);
01007      aHit->setIndex2Global(-1);
01008      aHit->setAssociatedTrackId(-1);
01009 
01010      int trayId = aHit->tray();
01011      int moduleId = aHit->module();
01012      int cellId = aHit->cell();
01013 
01014      if(Debug()) { LOG_INFO <<"A: fired hit in " << " tray="<< trayId << " module="<<moduleId<<" cell="<<cellId<<endm; }
01015 
01016      StructCellHit aDaqCellHit;
01017      aDaqCellHit.tray = trayId;
01018      aDaqCellHit.module = moduleId;
01019      aDaqCellHit.cell = cellId;
01020      aDaqCellHit.tot = aHit->tot();
01021      aDaqCellHit.index2BTofHit = i;
01022      daqCellsHitVec.push_back(aDaqCellHit);
01023 
01024      // additional valid number configuration
01025      int id = trayId*100+moduleId;
01026      //fg bool ifind = kFALSE;
01027      //fg for(size_t im=0;im<validModuleVec.size();im++) {
01028      //fg   if(id==validModuleVec[im]) {
01029      //fg     ifind = kTRUE;
01030      //fg     break;
01031      //fg   }
01032      //fg }
01033      //fg if(!ifind) validModuleVec.push_back(id);
01034      if (find(validModuleVec.begin(), validModuleVec.end(), id) == validModuleVec.end())
01035        validModuleVec.push_back(id);
01036 
01037      if(mHisto) {
01038        mDaqOccupancy[trayId-1]->Fill((moduleId-1)*mNCell+(cellId-1));
01039      }      
01040   }
01041       
01042   // end of Sect.A
01043   if(Debug()) {
01044     LOG_DEBUG << "    total # of cells = " << daqCellsHitVec.size() << endm;
01045     for(size_t iv = 0;iv<validModuleVec.size();iv++) {
01046       LOG_DEBUG << " module # " << validModuleVec[iv] << " Valid! " << endm;
01047     }
01048   }
01049   if(mHisto) {
01050     mCellsMultInEvent->Fill(daqCellsHitVec.size());
01051     if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
01052   }
01053 //  if(!daqCellsHitVec.size()) return;
01054   if (doPrintCpuInfo) {
01055     timer.stop();
01056     LOG_INFO << "CPU time after Step A - loading hits : "
01057          << timer.elapsedTime() << " sec" << endm;
01058     timer.start();    
01059   }
01060 
01061 
01062   //.........................................................................
01064   //
01065   Int_t index2Primary[50000]; // map
01066   memset(index2Primary, -1, sizeof(index2Primary));
01067   for(int ip=0;ip<(int)mMuDst->array(muPrimary)->GetEntries();ip++) {
01068     StMuTrack *pTrack = (StMuTrack *)mMuDst->array(muPrimary)->UncheckedAt(ip);
01069     if(!pTrack) continue;
01070     int gIndex = pTrack->index2Global();
01071     if(gIndex<0) continue;
01072     index2Primary[gIndex] = ip;
01073   }
01074 
01075   tofCellHitVector allCellsHitVec;
01076   StructCellHit cellHit;
01077 
01078   Int_t nGlobals = mMuDst->numberOfGlobalTracks();
01079   Int_t nAllTracks=0;
01080   Int_t nPrimaryHits = 0;
01081   for (int iNode=0; iNode<nGlobals; iNode++){
01082     tofCellHitVector cellHitVec;
01083     //    cellHitVec.clear();
01084     StMuTrack *theTrack = mMuDst->globalTracks(iNode);
01085     if(!theTrack) continue;
01086 
01087     // clean up any association done before
01088     StMuBTofPidTraits pidTof;
01089     theTrack->setBTofPidTraits(pidTof);
01090     theTrack->setIndex2BTofHit(-1);
01091 
01092     bool isPrimary = kFALSE;
01093     int pIndex = index2Primary[iNode];
01094     if(pIndex>=0) { 
01095       isPrimary = kTRUE;
01096       StMuTrack *pTrack = (StMuTrack *)mMuDst->array(muPrimary)->UncheckedAt(pIndex);
01097       if(pTrack) {
01098         pTrack->setBTofPidTraits(pidTof);
01099         pTrack->setIndex2BTofHit(-1);
01100       }
01101     }
01102 
01103     StThreeVectorF mom = theTrack->momentum();
01104     float pt = mom.perp();
01105     float eta = mom.pseudoRapidity();
01106     float phi = mom.phi();
01107     //fg if (phi<0.) phi += 2.*3.14159;
01108     if (phi<0.) phi += 2.*M_PI;
01109 
01110     float nSigmaPion = theTrack->nSigmaPion();
01111     float dEdx = theTrack->dEdx();
01112     int ndEdxpts = theTrack->nHitsDedx();
01113     int nfitpts = theTrack->nHitsFit(kTpcId);
01114 
01115     // make sure we have a track, a miniDST might have removed it...
01116     if (validTrack(theTrack)){
01117       if(mHisto) {
01118         mTrackPtEta->Fill(pt, eta);
01119         mTrackPtPhi->Fill(pt, phi);
01120         mTrackNFitPts->Fill(nfitpts);
01121         if(dEdx>0.) mTrackdEdxvsp->Fill(mom.mag(), dEdx*1.e6);
01122         if(fabs(nSigmaPion)<5.) mNSigmaPivsPt->Fill(pt, nSigmaPion+5.*theTrack->charge());
01123       }
01124 
01125       if(mSaveTree) {
01126         trackTree.pt = pt;
01127         trackTree.eta = eta;
01128         trackTree.phi = phi;
01129         trackTree.nfitpts = nfitpts;
01130         trackTree.dEdx = dEdx*1.e6;
01131         trackTree.ndEdxpts = ndEdxpts;
01132         trackTree.charge = theTrack->charge();      
01133         trackTree.projTrayId = 0;
01134         trackTree.projCellChan = -1;
01135         trackTree.projY = -999.;
01136         trackTree.projZ = -999.;
01137       }
01138 
01139       nAllTracks++;
01140       StPhysicalHelixD theHelix = mOuterTrackGeometry ? theTrack->outerHelix() : theTrack->helix();
01141 
01142 //      IntVec projTrayVec;
01143 //      if(!mBTofGeom->projTrayVector(theHelix, projTrayVec)) continue;
01144 
01145       IntVec idVec;
01146       DoubleVec pathVec, thetaVec;
01147       PointVec  crossVec;
01148 
01149 //       idVec.clear();
01150 //       pathVec.clear();
01151 //       crossVec.clear();
01152 
01153       Int_t ncells = 0;
01154       if(mBTofGeom->HelixCrossCellIds(theHelix,idVec,pathVec,crossVec,thetaVec) ) {
01155 //      if(mBTofGeom->HelixCrossCellIds(theHelix, validModuleVec, projTrayVec, idVec, pathVec, crossVec)) {
01156         Int_t cells = idVec.size();
01157         for (Int_t i=0; i<cells; i++) {
01158             Int_t icell,imodule,itray;
01159             Double_t local[3],global[3];
01160             for(Int_t i2=0;i2<3;i2++){
01161                  local[i2]=0;
01162             }
01163             global[0]=crossVec[i].x();
01164             global[1]=crossVec[i].y();
01165             global[2]=crossVec[i].z();
01166             mBTofGeom->DecodeCellId(idVec[i], icell, imodule, itray);
01167 //          LOG_INFO << " decode " << idVec[i] << "  to tray#" << itray << " module#" << imodule << " cell#" << icell << endm;
01168             StBTofGeomSensor* sensor = mBTofGeom->GetGeomSensor(imodule,itray);
01169             if(!sensor) {
01170               LOG_WARN << " No sensitive module in the projection??? -- Something weird!!! " << endm;
01171               continue;
01172             }
01173             sensor->Master2Local(&global[0],&local[0]);
01174             icell = sensor->FindCellIndex(local);
01175             //      StThreeVectorD glo=sensor->GetCenterPosition();
01176             StThreeVectorD glo(global[0], global[1], global[2]);
01177             StThreeVectorD hitPos(local[0], local[1], local[2]);
01178 //          delete sensor;   /// function in StBTofGeometry modified. Don't delete here.
01179             if (fabs(local[2])<mZLocalCut) {   // to be consistent with GEANT geometry
01180                                          // and the alignment calibration
01181 //          if (local[2]<=2.7&&local[2]>=-3.4) {
01182 //            if (fabs(local[2])<4.) {   // loose local z cut at first step
01183               ncells++;
01184               cellHit.tray = itray;
01185               cellHit.module = imodule;
01186               cellHit.cell = icell;
01187               cellHit.trackIdVec.push_back(iNode);
01188               cellHit.hitPosition = glo;        // global position
01189               cellHit.zhit = (Float_t)hitPos.z();
01190               cellHit.yhit = (Float_t)hitPos.y();
01191               cellHit.theta = (Double_t)thetaVec[i];
01192               cellHitVec.push_back(cellHit);
01193               allCellsHitVec.push_back(cellHit);
01194 
01195               if(isPrimary) nPrimaryHits++;
01196 
01197               if(mHisto) {
01198                 mDaqOccupancyProj[itray-1]->Fill((imodule-1)*mNCell+(icell-1));
01199                 mHitsPosition->Fill(hitPos.y(), hitPos.z());
01200               }
01201               
01202               if(mSaveTree) {
01203                 trackTree.projTrayId = itray;
01204                 trackTree.projCellChan = (imodule-1)*mNCell+(icell-1);
01205                 trackTree.projY = local[1];
01206                 trackTree.projZ = local[2];
01207               }
01208 
01209               if(Debug()) {
01210                 LOG_DEBUG <<"B: nodeid=" << iNode << "  projected in " << " tray="<< itray << " module="<<imodule<<" cell="<<icell<<endm;
01211                 LOG_DEBUG <<"   hit position " << hitPos << endm;
01212                 // LOG_DEBUG <<"   momemtum= " << pt << " " << eta << " " << phi << endm;
01213               }
01214             }
01215         } // for (Int_t i=0...)
01216       } // endif(helixcross...)
01217       if(ncells>0&&mHisto) mHitsMultPerTrack->Fill(ncells);
01218 
01219       if(mHisto && mSaveTree) mTrackTree->Fill();
01220 
01221     } // if(ValidTrack).. 
01222   } // loop over nodes
01223   if(Debug()) { LOG_INFO << "B:  matched/available/total #tracknodes: " <<allCellsHitVec.size() << "/" <<nAllTracks << "/" << nGlobals << endm; }
01224   if(mHisto) {
01225     mHitsMultInEvent->Fill(allCellsHitVec.size());
01226     mHitsPrimaryInEvent->Fill(nPrimaryHits);
01227     if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
01228   }
01229   // end of Sect.B
01230   if (doPrintCpuInfo) {
01231     timer.stop();
01232     LOG_INFO << "CPU time after Step B - projection : "
01233          << timer.elapsedTime() << " sec" << endm;
01234     timer.start();    
01235   }
01236 
01237   //.........................................................................
01239   //
01240   tofCellHitVector matchHitCellsVec;
01241 
01242   tofCellHitVectorIter daqIter = daqCellsHitVec.begin();
01243   for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
01244     tofCellHitVectorIter proIter = allCellsHitVec.begin();
01245     for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
01246 
01247       int daqIndex = (daqIter->module-1)*6 + (daqIter->cell-1);
01248       int proIndex = (proIter->module-1)*6 + (proIter->cell-1);
01249       int hisIndex = daqIter->tray - 1;
01250 //      int daqAllIndex = (daqIter->tray - 1)*mNTOF + daqIndex;
01251 //      int proAllIndex = (proIter->tray - 1)*mNTOF + proIndex;
01252       if(daqIter->tray==proIter->tray) {
01253         if (mHisto) {
01254           if(hisIndex>=0&&hisIndex<mNTray) {
01255             mHitCorr[hisIndex]->Fill(proIndex,daqIndex);
01256             mHitCorrModule[hisIndex]->Fill(proIter->module-1,daqIter->module-1);
01257           } else {
01258             LOG_WARN << " weird tray # " << daqIter->tray << endm;
01259           }
01260         }
01261       }
01262       if( (daqIter->tray==proIter->tray)&& 
01263           (daqIter->module==proIter->module) &&
01264           ( ( (proIter->cell==6)&&((proIter->cell==daqIter->cell) ||
01265                                    (proIter->cell==daqIter->cell+1)) )
01266             || ( (proIter->cell==1)&&((proIter->cell==daqIter->cell) ||
01267                                       (proIter->cell==daqIter->cell-1)) )
01268             || ( (proIter->cell>=2&&proIter->cell<=6) &&
01269                  ( (proIter->cell==daqIter->cell) ||
01270                    (proIter->cell==daqIter->cell-1) ||
01271                    (proIter->cell==daqIter->cell+1) ) ) ) ) {
01272         cellHit.tray = daqIter->tray;
01273         cellHit.module = daqIter->module;
01274         cellHit.cell = daqIter->cell;
01275         cellHit.hitPosition = proIter->hitPosition;
01276         cellHit.trackIdVec = proIter->trackIdVec;
01277         cellHit.zhit = proIter->zhit;
01278         cellHit.yhit = proIter->yhit;
01279     cellHit.tot = daqIter->tot;
01280     cellHit.index2BTofHit = daqIter->index2BTofHit;
01281     cellHit.theta = proIter->theta;
01282         matchHitCellsVec.push_back(cellHit);
01283       }
01284     }
01285   } //end {sec. C}
01286   if(Debug()) { LOG_INFO << "C: before/after: " << allCellsHitVec.size() << "/" << matchHitCellsVec.size() << endm; }
01287   if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
01288   if (doPrintCpuInfo) {
01289     timer.stop();
01290     LOG_INFO << "CPU time after Step C - matching : "
01291          << timer.elapsedTime() << " sec" << endm;
01292     timer.start();    
01293   }
01294 
01295   //.........................................................................
01297   //
01298   Int_t nSingleHitCells(0);
01299   Int_t nMultiHitsCells(0);
01300 
01301   tofCellHitVector singleHitCellsVec;
01302   tofCellHitVector multiHitsCellsVec;
01303 
01304   tofCellHitVector tempVec = matchHitCellsVec;
01305   tofCellHitVector erasedVec = tempVec;
01306   while (tempVec.size() != 0) {
01307     Int_t nTracks = 0;
01308     idVector trackIdVec;
01309 
01310     tofCellHitVectorIter tempIter=tempVec.begin();
01311     tofCellHitVectorIter erasedIter=erasedVec.begin();
01312     while(erasedIter!= erasedVec.end()) {
01313       if(tempIter->tray == erasedIter->tray &&
01314          tempIter->module == erasedIter->module &&
01315          tempIter->cell == erasedIter->cell) {
01316         nTracks++;
01317         trackIdVec.push_back(erasedIter->trackIdVec.back());  // merge
01318         erasedVec.erase(erasedIter);
01319         erasedIter--;
01320       }
01321       erasedIter++;
01322     }
01323 
01324     cellHit.cell = tempIter->cell;
01325     cellHit.module = tempIter->module;
01326     cellHit.tray = tempIter->tray;
01327     cellHit.hitPosition = tempIter->hitPosition;
01328     cellHit.trackIdVec = trackIdVec;
01329     cellHit.zhit = tempIter->zhit;
01330     cellHit.yhit = tempIter->yhit;
01331     cellHit.tot = tempIter->tot;
01332     cellHit.index2BTofHit = tempIter->index2BTofHit;
01333     cellHit.theta = tempIter->theta;
01334 
01335     Float_t ycenter = (tempIter->cell-1-2.5)*mWidthPad;
01336     Float_t dy = tempIter->yhit - ycenter;
01337     Float_t dz = tempIter->zhit;
01338 
01339     if(mHisto) {
01340       mTracksPerCellMatch1->Fill(trackIdVec.size());
01341 //      mDaqOccupancyMatch1->Fill((tempIter->module-1)*mNCell+(tempIter->cell-1));
01342       mDeltaHitMatch1->Fill(dy, dz);
01343     }
01344 
01345     if (nTracks==1){
01346       nSingleHitCells++;      
01347       singleHitCellsVec.push_back(cellHit);
01348     } else if (nTracks>1){
01349       nMultiHitsCells++;
01350       multiHitsCellsVec.push_back(cellHit);
01351       // for multiple hit cells either discard (yes) or
01352       // find the most likely candidate.
01353     } else {
01354       LOG_WARN << "D: no tracks extrapolate to matched cell ... should not happen!" << endm;
01355     }
01356     
01357     if(Debug()) {
01358       LOG_DEBUG << "D: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
01359       idVectorIter ij=trackIdVec.begin();
01360       while (ij != trackIdVec.end()) { LOG_DEBUG << " " << *ij; ij++; }
01361       LOG_DEBUG << endm;
01362     }
01363 
01364     tempVec = erasedVec;
01365   }
01366   if(Debug()) { LOG_INFO << "D: before/after: " << matchHitCellsVec.size() << "/" << singleHitCellsVec.size() << endm; }
01367   //end of Sect.D
01368 
01369   if(mHisto) {
01370     mCellsPerEventMatch1->Fill(singleHitCellsVec.size()+multiHitsCellsVec.size());
01371     if(singleHitCellsVec.size()) mEventCounterHisto->Fill(9);
01372   } 
01373   if (doPrintCpuInfo) {
01374     timer.stop();
01375     LOG_INFO << "CPU time after Step D - erasing : "
01376          << timer.elapsedTime() << " sec" << endm;
01377     timer.start();    
01378   }
01379 
01380   //.........................................................................
01382   //
01383   tofCellHitVector FinalMatchedCellsVec;
01384   //  FinalMatchedCellsVec.clear();
01385   tempVec = singleHitCellsVec;
01386   if(mHisto) {
01387     mCellsPerEventMatch2->Fill(tempVec.size());
01388     for(unsigned int ii=0;ii<tempVec.size();ii++) {
01389       mTracksPerCellMatch2->Fill(tempVec[ii].trackIdVec.size());
01390 //      mDaqOccupancyMatch2->Fill((tempVec[ii].module-1)*mNCell+(tempVec[ii].cell-1));
01391       Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
01392       Float_t dy = tempVec[ii].yhit-ycenter;
01393       Float_t dz = tempVec[ii].zhit;
01394       mDeltaHitMatch2->Fill(dy, dz);
01395     }
01396   }
01397 
01398   erasedVec = tempVec;
01399   while (tempVec.size() != 0) {
01400     StructCellHit cellHit;
01401     Int_t nCells = 0;
01402     idVector vTrackId;
01403     vector<StThreeVectorD> vPosition;
01404     vector<Int_t> vchannel, vtray, vmodule, vcell;
01405     vector<Float_t> vzhit, vyhit;
01406     vector<Double_t> vtot, vtheta;
01407     vector<Int_t> vindex2BTofHit;
01408 
01409     tofCellHitVectorIter tempIter=tempVec.begin();
01410     tofCellHitVectorIter erasedIter=erasedVec.begin();
01411     while(erasedIter!= erasedVec.end()) {
01412       if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
01413         nCells++;
01414         vtray.push_back(erasedIter->tray);
01415         vmodule.push_back(erasedIter->module);
01416         vcell.push_back(erasedIter->cell);
01417         vPosition.push_back(erasedIter->hitPosition);
01418         vTrackId.push_back(erasedIter->trackIdVec.back());
01419         vzhit.push_back(erasedIter->zhit);
01420         vyhit.push_back(erasedIter->yhit);
01421     vtot.push_back(erasedIter->tot);
01422     vindex2BTofHit.push_back(erasedIter->index2BTofHit);
01423     vtheta.push_back(erasedIter->theta);
01424 
01425         erasedVec.erase(erasedIter);
01426         erasedIter--;
01427       }
01428       erasedIter++;
01429     }
01430 
01431     if (nCells==1){
01432       // for singly hit cell, copy data in singleHitCellsVec
01433       cellHit.tray = vtray[0];
01434       cellHit.module = vmodule[0];
01435       cellHit.cell = vcell[0];
01436       cellHit.trackIdVec.push_back(vTrackId[0]);
01437       cellHit.hitPosition = vPosition[0];
01438       cellHit.matchFlag = 1; 
01439       cellHit.zhit = vzhit[0];
01440       cellHit.yhit = vyhit[0];
01441       cellHit.tot = vtot[0];
01442       cellHit.index2BTofHit = vindex2BTofHit[0];
01443       cellHit.theta = vtheta[0];
01444 
01445       FinalMatchedCellsVec.push_back(cellHit);
01446 
01447       // debugging output
01448       if(Debug()) {
01449         LOG_DEBUG << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
01450         idVectorIter ij=vTrackId.begin();
01451         while (ij != vTrackId.end()) { LOG_DEBUG << " " << *ij; ij++; }
01452         LOG_DEBUG << endm;
01453       }
01454       
01455     }
01456     else if (nCells>1){   // for multiple hit cells  find the most likely candidate.
01457       Int_t thiscandidate(-99);
01458       Int_t thisMatchFlag(0);
01459 
01460       // sort on tot
01461       Float_t tot(0.);
01462       vector<Int_t> ttCandidates;
01463       for (Int_t i=0;i<nCells;i++) {
01464         Double_t tt = vtot[i];
01465         if(tt<40.&&tt>tot) {    // open the ToT to 40 ns
01466           tot = tt;
01467           ttCandidates.clear();
01468           ttCandidates.push_back(i);
01469         } else if (tt==tot) {
01470           ttCandidates.push_back(i);
01471         }
01472       }
01473       if (ttCandidates.size()==1) {
01474         thiscandidate = ttCandidates[0];
01475         thisMatchFlag = 2;
01476       } else if (ttCandidates.size()>1) {  // sort on hitposition
01477         Float_t ss(99.);
01478         vector<Int_t> ssCandidates;
01479         for(size_t j=0;j<ttCandidates.size();j++) {
01480           Float_t yy = vyhit[ttCandidates[j]];
01481           Float_t ycell = (vcell[ttCandidates[j]]-1-2.5)*mWidthPad;
01482           Float_t ll = fabs(yy-ycell);
01483           if(ll<ss) {
01484             ss = ll; 
01485             ssCandidates.clear();
01486             ssCandidates.push_back(ttCandidates[j]);
01487           }else if  (ll==ss)
01488             ssCandidates.push_back(ttCandidates[j]);
01489         }
01490         if (ssCandidates.size()==1){
01491           thiscandidate = ssCandidates[0];
01492           thisMatchFlag = 3;
01493         }
01494       }
01495 
01496       if (thiscandidate>=0) {
01497         cellHit.tray = vtray[thiscandidate];
01498         cellHit.module = vmodule[thiscandidate];
01499         cellHit.cell = vcell[thiscandidate];
01500         cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
01501         cellHit.hitPosition = vPosition[thiscandidate];
01502         cellHit.matchFlag = thisMatchFlag;
01503         cellHit.zhit = vzhit[thiscandidate];
01504         cellHit.yhit = vyhit[thiscandidate];
01505     cellHit.tot = vtot[thiscandidate];
01506     cellHit.index2BTofHit = vindex2BTofHit[thiscandidate];
01507     cellHit.theta = vtheta[thiscandidate];
01508 
01509         FinalMatchedCellsVec.push_back(cellHit);
01510         
01511         // debugging output
01512         if(Debug()) { LOG_DEBUG << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:" << vTrackId[thiscandidate] << endm; }
01513       }
01514 
01515     } else {
01516       LOG_WARN << "E: no cells belong to this track ... should not happen!" << endm;
01517     }
01518 
01519     tempVec = erasedVec;
01520   }
01521 
01522   if(Debug()) { LOG_INFO << "E: before/after: " << singleHitCellsVec.size() << "/" << FinalMatchedCellsVec.size() << endm; }
01523   // end of Sect.E
01524   if (doPrintCpuInfo) {
01525     timer.stop();
01526     LOG_INFO << "CPU time after Step E - sorting : "
01527          << timer.elapsedTime() << " sec" << endm;
01528     timer.start();    
01529   }
01530 
01531   //.........................................................................
01533   //
01534   if(mHisto) {
01535     if(FinalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
01536     mCellsPerEventMatch3->Fill(FinalMatchedCellsVec.size());
01537   }
01538 
01539 //  StSPtrVecBTofHit& tofHits = theTof->tofHits();
01540   Int_t nValidSingleHitCells(0), nValidSinglePrimHitCells(0);
01541 
01542   for (size_t ii=0; ii < FinalMatchedCellsVec.size(); ii++){
01543     Int_t tray = FinalMatchedCellsVec[ii].tray;
01544     Int_t module = FinalMatchedCellsVec[ii].module;
01545     Int_t cell = FinalMatchedCellsVec[ii].cell;
01546 
01547     Float_t ycenter = (cell-1-2.5)*mWidthPad;
01548     Float_t dy = FinalMatchedCellsVec[ii].yhit - ycenter;
01549     Float_t dz = FinalMatchedCellsVec[ii].zhit;
01550     if (FinalMatchedCellsVec[ii].trackIdVec.size()!=1)
01551       LOG_WARN << "F: WHAT!?!  mult.matched cell in single cell list " << tray << " " << module << " " << cell << endm;
01552 
01553     if(mHisto) {
01554       mTracksPerCellMatch3->Fill(FinalMatchedCellsVec[ii].trackIdVec.size());
01555 //      mDaqOccupancyMatch3->Fill((module-1)*mNCell+(cell-1));
01556       mDeltaHitMatch3->Fill(dy, dz);
01557       mDeltaHitFinal[tray-1]->Fill(dy,dz);      
01558     }
01559 
01560     // get track-id from cell hit vector
01561     int trackNode = FinalMatchedCellsVec[ii].trackIdVec[0];
01562     StMuTrack *globalTrack = mMuDst->globalTracks(trackNode);
01563     if(!globalTrack) {
01564       LOG_WARN << "Wrong global track!" << endm;
01565       continue;
01566     }
01567 
01568     // Fill association in TOF Hit Collection
01569     StMuBTofHit *tofHit = mMuDst->btofHit(FinalMatchedCellsVec[ii].index2BTofHit);
01570     if(tofHit->tray()!=tray || tofHit->module()!=module || tofHit->cell()!=cell) {
01571       LOG_WARN << "Wrong hit in the MuBTofHit!" << endm;
01572       continue;
01573     }
01574     nValidSingleHitCells++;
01575 
01577     tofHit->setAssociatedTrackId(globalTrack->id());
01578     tofHit->setIndex2Global(trackNode);
01579     globalTrack->setIndex2BTofHit(FinalMatchedCellsVec[ii].index2BTofHit);
01580 
01581     int ip = index2Primary[trackNode];
01582     StMuTrack *primaryTrack = 0;
01583     if(ip>=0) {
01584       nValidSinglePrimHitCells++;
01585       tofHit->setIndex2Primary(ip);
01586       primaryTrack = (StMuTrack *)mMuDst->array(muPrimary)->UncheckedAt(ip);
01587       if(primaryTrack) {
01588         primaryTrack->setIndex2BTofHit(FinalMatchedCellsVec[ii].index2BTofHit);
01589       }
01590     }
01591 
01592     // Fill the matched data in StBTofPidTraits
01593     StMuBTofPidTraits pidTof = globalTrack->btofPidTraits();
01594     pidTof.setMatchFlag(FinalMatchedCellsVec[ii].matchFlag);
01595     pidTof.setYLocal(dy);
01596     pidTof.setZLocal(FinalMatchedCellsVec[ii].zhit);
01597     pidTof.setPosition(FinalMatchedCellsVec[ii].hitPosition);
01598     pidTof.setThetaLocal(FinalMatchedCellsVec[ii].theta);
01599     globalTrack->setBTofPidTraits(pidTof);
01600 
01601     if(primaryTrack) {
01602       StMuBTofPidTraits ppidTof = primaryTrack->btofPidTraits();
01603       ppidTof.setMatchFlag(FinalMatchedCellsVec[ii].matchFlag);
01604       ppidTof.setYLocal(dy);
01605       ppidTof.setZLocal(FinalMatchedCellsVec[ii].zhit);
01606       ppidTof.setPosition(FinalMatchedCellsVec[ii].hitPosition);
01607       ppidTof.setThetaLocal(FinalMatchedCellsVec[ii].theta);
01608       primaryTrack->setBTofPidTraits(ppidTof);
01609     }
01610 
01611   } // end final matched cells
01612   
01613   if(mHisto) {
01614     mCellsPrimaryPerEventMatch3->Fill(nValidSinglePrimHitCells);
01615   }
01616   
01617   if(Debug()) { LOG_INFO << "F: before/after" << FinalMatchedCellsVec.size() << "/" <<nValidSinglePrimHitCells << endm; }
01618  // end of Sect.F
01619   if (doPrintCpuInfo) {
01620     timer.stop();
01621     LOG_INFO << "CPU time after Step F - final : "
01622          << timer.elapsedTime() << " sec" << endm;
01623     timer.start();    
01624   }
01625 
01626   LOG_INFO << " #(daq hits): " << daqCellsHitVec.size()
01627        << "\t#(proj hits): " << allCellsHitVec.size()
01628        << "\t#(prim proj hits): " << nPrimaryHits
01629        << "\n#(matched hits): " << FinalMatchedCellsVec.size() 
01630        << "\n#(single valid hits): " << nValidSingleHitCells
01631        << "\t#(single prim valid hits): " << nValidSinglePrimHitCells
01632        << endm;
01633 
01634   if (doPrintMemoryInfo) {
01635         StMemoryInfo::instance()->snapshot();
01636         StMemoryInfo::instance()->print();
01637   }
01638   if (doPrintCpuInfo) {
01639     timer.stop();
01640     LOG_INFO << "CPU time for StBTofMatchMaker::Make(): "
01641          << timer.elapsedTime() << " sec" << endm;
01642   }
01643 
01644   LOG_INFO << "StBTofMatchMaker -- bye-bye" << endm;
01645 
01646   return;
01647 }
01648 
01649 //---------------------------------------------------------------------------
01650 // Book histograms and create ordered collection for easy manipulation
01651 void StBTofMatchMaker::bookHistograms(void){
01652 
01653   mEventCounterHisto = new TH1D("eventCounter","eventCounter",20,0,20);
01654 
01655   mCellsMultInEvent = new TH1D("cellsPerEvent","cellsPerEvent",1000,0,1000);
01656   mHitsMultInEvent  = new TH1D("hitsPerEvent","hitsPerEvent",1000,0,1000);
01657   mHitsPrimaryInEvent  = new TH1D("hitsPrimaryPerEvent","hitsPrimaryPerEvent",1000,0,1000);
01658   mHitsMultPerTrack = new TH1D("hitsPerTrack","hitsPerTrack",10,0,10);
01659   mHitsPosition     = new TH2D("hitsPosition","hitsPositions",300,-15.,15.,200,-5.,5.);
01660 
01661   // occupancy
01662   for(int i=0;i<mNTray;i++) {
01663     char hisname[100];
01664     sprintf(hisname,"Occupancy_Tray_%d",i+1);
01665     mDaqOccupancy[i]     = new TH1D(hisname,"",192,0,192);
01666     sprintf(hisname,"OccupancyProj_Tray_%d",i+1);
01667     mDaqOccupancyProj[i] = new TH1D(hisname,"",192,0,192);
01668   }
01669 
01670   // correlation
01671   for(int i=0;i<mNTray;i++) {
01672     char hisname[100];
01673     sprintf(hisname,"Corr_Tray_%d",i+1);
01674     mHitCorr[i] = new TH2D(hisname,"",192,0,192,192,0,192);
01675     sprintf(hisname,"Corr_Tray_%d_module",i+1);
01676     mHitCorrModule[i] = new TH2D(hisname,"",32,0,32,32,0,32);
01677   }
01678 
01679   // project hit position
01680   for(int i=0;i<mNTray;i++) {
01681     char hisname[100];
01682     sprintf(hisname,"LocalYZ_Tray_%d",i+1);
01683     mDeltaHitFinal[i] = new TH2D(hisname,"",300,-15.,15.,200,-5.,5.);
01684   }
01685 
01686   // TPC tracks
01687   if(mSaveTree) {
01688     mTrackTree = new TTree("track","track");
01689     mTrackTree->Branch("pt",&trackTree.pt,"pt/F");
01690     mTrackTree->Branch("eta",&trackTree.eta,"eta/F");
01691     mTrackTree->Branch("phi",&trackTree.phi,"phi/F");
01692     mTrackTree->Branch("nfitpts",&trackTree.nfitpts,"nfitpts/I");
01693     mTrackTree->Branch("dEdx",&trackTree.dEdx,"dEdx/F");
01694     mTrackTree->Branch("ndEdxpts",&trackTree.ndEdxpts,"ndEdxpts/I");
01695     mTrackTree->Branch("charge",&trackTree.charge,"charge/I");
01696     mTrackTree->Branch("projTrayId",&trackTree.projTrayId,"projTrayId/I");
01697     mTrackTree->Branch("projCellChan",&trackTree.projCellChan,"projCellChan/I");
01698     mTrackTree->Branch("projY",&trackTree.projY,"projY/F");
01699     mTrackTree->Branch("projZ",&trackTree.projZ,"projZ/F");
01700   }
01701 
01702   mTrackPtEta = new TH2D("trackPtEta","",100,0.,5.,60,-1.5,1.5);
01703   //mTrackPtPhi = new TH2D("trackPtPhi","",100,0.,5.,120,0.,2*3.14159);
01704   mTrackPtPhi = new TH2D("trackPtPhi","",100,0.,5.,120,0.,2.*M_PI);
01705   mTrackNFitPts = new TH1D("trackNFitPts","",50,0.,50.);
01706   mTrackdEdxvsp = new TH2D("trackdEdxvsp","",500,0.,5.,1000,0.,10.);
01707   mNSigmaPivsPt = new TH2D("nSigmaPivsPt","",500,0.,5.,1000,-10.,10.);
01708 
01709   // association  
01710   mCellsPerEventMatch1 = new TH1D("cellsPerEventMatch1","cellPerEventMatch1",100,0,100);
01711   mHitsPerEventMatch1 = new TH1D("hitsPerEventMatch1","hitsPerEventMatch1",100,0,100);
01712   mCellsPerTrackMatch1 = new TH1D("cellsPerTrackMatch1","cellsPerTrackMatch1",100,0,100);
01713   mTracksPerCellMatch1 = new TH1D("tracksPerCellMatch1","tracksPerCellMatch1",100,0,100);
01714   mDeltaHitMatch1 = new TH2D("deltaHitMatch1","deltaHitMatch1",300,-15,15,200,-5.,5.);
01715 
01716   // kick out multi-hit
01717   mCellsPerEventMatch2 = new TH1D("cellsPerEventMatch2","cellPerEventMatch2",100,0,100);
01718   mHitsPerEventMatch2 = new TH1D("hitsPerEventMatch2","hitsPerEventMatch2",100,0,100);
01719   mCellsPerTrackMatch2 = new TH1D("cellsPerTrackMatch2","cellsPerTrackMatch2",100,0,100);
01720   mTracksPerCellMatch2 = new TH1D("tracksPerCellMatch2","tracksPerCellMatch2",100,0,100);
01721   mDeltaHitMatch2 = new TH2D("deltaHitMatch2","deltaHitMatch2",300,-15,15,200,-5.,5.);
01722 
01723   // sort out multi matched cells
01724   mCellsPerEventMatch3 = new TH1D("cellsPerEventMatch3","cellsPerEventMatch3",100,0,100);
01725   mHitsPerEventMatch3 = new TH1D("hitsPerEventMatch3","hitsPerEventMatch3",100,0,100);
01726   mCellsPerTrackMatch3 = new TH1D("cellsPerTrackMatch3","cellsPerTrackMatch3",100,0,100);
01727   mTracksPerCellMatch3 = new TH1D("tracksPerCellMatch3","tracksPerCellMatch3",100,0,100);
01728   mDeltaHitMatch3 = new TH2D("deltaHitMatch3","deltaHitMatch3",300,-15,15,200,-5.,5.);
01729 
01730   mCellsPrimaryPerEventMatch3 = new TH1D("cellsPrimaryPerEventMatch3","cellsPrimaryPerEventMatch3",100,0,100);
01731 
01732   return;
01733 }
01734 
01735 
01736 //---------------------------------------------------------------------------
01737 // store histograms in a seperate root file
01738 void StBTofMatchMaker::writeHistogramsToFile(){
01739   // Output file
01740   TFile *theHistoFile =  new TFile(mHistoFileName.c_str(), "RECREATE");
01741   LOG_INFO << "StBTofMatchMaker::writeHistogramsToFile()"
01742        << " histogram file " <<  mHistoFileName << endm;
01743 
01744   theHistoFile->cd();
01745 
01746   if(mHisto) {
01747 
01748 
01749     for(int i=0;i<mNTray;i++) {
01750       mDaqOccupancy[i]->Write();
01751       mDaqOccupancyProj[i]->Write();
01752       mHitCorr[i]->Write();
01753       mHitCorrModule[i]->Write();
01754       mDeltaHitFinal[i]->Write();
01755     }
01756 
01757     mEventCounterHisto->Write();
01758     mCellsMultInEvent->Write();
01759     mHitsMultInEvent->Write();
01760     mHitsPrimaryInEvent->Write();
01761     mHitsMultPerTrack->Write();
01762     mHitsPosition->Write();
01763 
01764     mTrackPtEta->Write();
01765     mTrackPtPhi->Write();
01766     mTrackNFitPts->Write();
01767     mTrackdEdxvsp->Write();
01768     mNSigmaPivsPt->Write();
01769 
01770     mCellsPerEventMatch1->Write();
01771     mHitsPerEventMatch1->Write();
01772     mCellsPerTrackMatch1->Write();
01773     mTracksPerCellMatch1->Write();
01774     mDeltaHitMatch1->Write();
01775     
01776     mCellsPerEventMatch2->Write();
01777     mHitsPerEventMatch2->Write();
01778     mCellsPerTrackMatch2->Write();
01779     mTracksPerCellMatch2->Write();
01780     mDeltaHitMatch2->Write();
01781 
01782     mCellsPerEventMatch3->Write();
01783     mHitsPerEventMatch3->Write();
01784     mCellsPerTrackMatch3->Write();
01785     mTracksPerCellMatch3->Write();
01786     mDeltaHitMatch3->Write();
01787 
01788     mCellsPrimaryPerEventMatch3->Write();
01789     
01790     theHistoFile->Write();  
01791     theHistoFile->Close(); 
01792    
01793   }
01794 
01795   if(mSaveTree) {
01796     TFile *theTreeFile =  new TFile((mHistoFileName+".tree.root").c_str(), "RECREATE");
01797     theTreeFile->cd();
01798     mTrackTree->Write();
01799     theTreeFile->Write();
01800     theTreeFile->Close();
01801   }
01802 
01803   return;
01804 }
01805 
01806 //---------------------------------------------------------------------------
01807 // determine whether this is a valid TPC track
01808 bool StBTofMatchMaker::validTrack(StTrack *track){
01809   // 1. no track, no go.
01810   if (!track) return false;
01811 
01812   // 2. track quality flag, should be >0
01813   if (track->flag()<=0 || track->flag()>=1000) return false;
01814 
01815   // 3. minimum #hits per track - obsolete
01816   //  if (track->topologyMap().numberOfHits(kTpcId) < mMinHitsPerTrack) return false;
01817   // 4. minimum #fit points per track
01818   if (track->fitTraits().numberOfFitPoints(kTpcId) < mMinFitPointsPerTrack) return false;
01819   // 5. minimum #fit points over #maximum points
01820   //fg float ratio = (1.0*track->fitTraits().numberOfFitPoints(kTpcId)) / (1.0*track->numberOfPossiblePoints(kTpcId));
01821   float ratio = (float)track->fitTraits().numberOfFitPoints(kTpcId) / (1.0*track->numberOfPossiblePoints(kTpcId));
01822   if (ratio < mMinFitPointsOverMax) return false;
01823 
01824   return true;
01825 }
01826 
01827 //---------------------------------------------------------------------------
01828 // determine whether this is a valid TPC track
01829 bool StBTofMatchMaker::validTrack(StMuTrack *track){
01830   // 1. no track, no go.
01831   if (!track) return false;
01832 
01833   // 2. track quality flag, should be >0
01834   if (track->flag()<=0 || track->flag()>=1000) return false;
01835 
01836   // 3. minimum #hits per track - obsolete
01837   //  if (track->nHits() < mMinHitsPerTrack) return false;
01838   // 4. minimum #fit points per track
01839   if (track->nHitsFit(kTpcId) < mMinFitPointsPerTrack) return false;
01840   // 5. minimum #fit points over #maximum points
01841   //fg float ratio = (1.0*track->fitTraits().numberOfFitPoints(kTpcId)) / (1.0*track->numberOfPossiblePoints(kTpcId));
01842   float ratio = (float)track->nHitsFit(kTpcId) / (1.0*track->nHitsPoss(kTpcId));
01843   if (ratio < mMinFitPointsOverMax) return false;
01844 
01845   return true;
01846 }
01847 
01848 //---------------------------------------------------------------------------
01849 // returns the proper track geometry, based on a global user setting
01850 StTrackGeometry* StBTofMatchMaker::trackGeometry(StTrack* track){
01851   // returns apropriate StTrackGeometry (standard or outerGeometry)
01852   if (!track) return 0;
01853   StTrackGeometry *thisTrackGeometry;
01854   if (mOuterTrackGeometry)
01855     thisTrackGeometry = track->outerGeometry();
01856   else
01857     thisTrackGeometry = track->geometry();
01858   return thisTrackGeometry;
01859 }