StRoot: StTofrMatchMaker/StTofrMatchMaker.cxx Source File

00001 /*******************************************************************
00002  *
00003  * $Id: StTofrMatchMaker.cxx,v 1.30 2009/07/28 16:06:42 geurts Exp $
00004  *
00005  * Author: Xin Dong
00006  *****************************************************************
00007  *
00008  * Description: Tofr Match Maker to do the matching between the 
00009  *              fired celles and TPC tracks ( similar to  Frank's
00010  *              TofpMatchMaker )
00011  *
00012  *****************************************************************
00013  *
00014  * $Log: StTofrMatchMaker.cxx,v $
00015  * Revision 1.30  2009/07/28 16:06:42  geurts
00016  * Bug Ticket #1591: explicit initialization of mStrobeTdcMin, mStrobeTdcMax, and mPedTOFr.
00017  *
00018  * Revision 1.29  2009/07/24 22:42:08  fine
00019  * replace the deprecated API of the STAR messenger
00020  *
00021  * Revision 1.28  2009/07/24 22:33:33  fine
00022  * Make the code C++ compliant
00023  *
00024  * Revision 1.27  2009/06/09 19:45:35  jeromel
00025  * Changes for BT#1428
00026  *
00027  * Revision 1.26  2008/07/23 19:22:03  dongx
00028  * New track quality cuts for Run8
00029  *
00030  * Revision 1.25  2008/06/24 21:58:13  dongx
00031  * fixed a bug of crashing due to potential empty track in trackNodes
00032  *
00033  * Revision 1.24  2008/06/19 16:11:54  dongx
00034  * fixed bug of filling an undefined histogram in case of production
00035  *
00036  * Revision 1.23  2008/05/12 17:16:37  dongx
00037  * letime and tetime in StTofCell stored in nano-seconds
00038  *
00039  * Revision 1.22  2008/05/06 18:41:39  dongx
00040  * - Fixed bug in ouput histogram filename switch
00041  * - Added switch for tpc track tree output
00042  *
00043  * Revision 1.21  2008/04/23 18:20:15  dongx
00044  * vpd letime and tetime stored in double precision
00045  *
00046  * Revision 1.20  2008/04/22 22:31:22  dongx
00047  * leadingEdgeTime and trailingEdgeTime stored as double precision in StTofCell
00048  *
00049  * Revision 1.19  2008/04/22 20:55:13  dongx
00050  * leadingEdgeTime and trailingEdgeTime stored as double precision in StTofCell
00051  *
00052  * Revision 1.18  2008/03/27 00:16:03  dongx
00053  *  update for Run8 finished.
00054  *
00055  * Revision 1.17  2007/11/29 22:43:11  dongx
00056  * changed vpd trayId definition to 121 (East) and 122 (West)
00057  *
00058  * Revision 1.16  2007/11/28 02:17:08  dongx
00059  * trayId for vpd in StTofCell: 901 (E), 902 (W)
00060  * dataIndex for vpd in StTofData:  use 121 and 122 as trayId
00061  *
00062  * Revision 1.15  2007/11/22 00:22:37  dongx
00063  * update for run8 - first version
00064  *
00065  * Revision 1.14  2007/04/17 23:02:20  dongx
00066  * replaced with standard STAR Loggers
00067  *
00068  * Revision 1.13  2007/02/28 23:31:59  dongx
00069  * completion for Run V matching
00070  *   - trailing tdc and leading tdc stored as adc and tdc in StTofCell
00071  *   - multi-hit association cases use hit position priority
00072  *
00073  * Revision 1.12  2005/04/12 17:31:56  dongx
00074  * update for year 5 data - not completed, leave as empty at present
00075  *
00076  * Revision 1.11  2004/08/11 18:57:44  dongx
00077  * deltay quality cut applied
00078  *
00079  * Revision 1.10  2004/08/10 19:25:13  dongx
00080  * updated to be compatible with RunII data
00081  *
00082  * Revision 1.9  2004/07/13 16:12:32  dongx
00083  * continuing update for the DAQ array in Run IV
00084  *
00085  * Revision 1.8  2004/07/12 17:21:19  dongx
00086  * fix the crashing when running events before day 030 (TofDaq not updated to 184 yet)
00087  *
00088  * Revision 1.7  2004/05/03 23:08:50  dongx
00089  * change according to the update of StTofrGeometry, save CPU time by 100 times
00090  *
00091  * Revision 1.5  2004/04/01 20:08:54  dongx
00092  * fix a bug about the hit position stored in TofCell
00093  *
00094  * Revision 1.4  2004/03/16 22:30:49  dongx
00095  * fix the warning message when compiling
00096  *
00097  * Revision 1.3  2004/03/11 22:30:34  dongx
00098  * -move m_Mode control to Init()
00099  * -clear up
00100  *
00101  * Revision 1.2  2004/03/09 17:44:56  dongx
00102  * first release
00103  *
00104  *
00105  *******************************************************************/
00106 #include <iostream>
00107 #include "StEventTypes.h"
00108 #include "Stypes.h"
00109 #include "StThreeVectorD.hh"
00110 #include "StThreeVectorF.hh"
00111 #include "StHelix.hh"
00112 #include "StTrackGeometry.h"
00113 #include "StDcaGeometry.h"
00114 #include "StDedxPidTraits.h"
00115 #include "StTrackPidTraits.h"
00116 #include "StarClassLibrary/StParticleTypes.hh"
00117 #include "StarClassLibrary/StParticleDefinition.hh"
00118 #include "StTpcDedxPidAlgorithm.h"
00119 #include "StEventUtilities/StuRefMult.hh"
00120 #include "PhysicalConstants.h"
00121 #include "StPhysicalHelixD.hh"
00122 #include "StHelix.hh"
00123 #include "StTofUtil/tofPathLength.hh"
00124 #include "StTofUtil/StTofrDaqMap.h"
00125 #include "StTofUtil/StTofINLCorr.h"
00126 #include "StTofUtil/StTofrGeometry.h"
00127 #include "StTofUtil/StTofCellCollection.h"
00128 #include "tables/St_pvpdStrobeDef_Table.h"
00129 #include "tables/St_tofPedestal_Table.h"
00130 #include "tables/St_tofConfig_Table.h"
00131 #include "tables/St_tofTrayConfig_Table.h"
00132 //#include "tables/St_tofr5INLtable_Table.h"
00133 #include "StTofUtil/tofPathLength.hh"
00134 #include "TH1.h"
00135 #include "TH2.h"
00136 #include "TFile.h"
00137 #include "TTree.h"
00138 #include "StMessMgr.h"
00139 #include "StMemoryInfo.hh"
00140 #include "StTimer.hh"
00141 #include "StTofrMatchMaker.h"
00142 //#include "TMemStat.h"
00143 
00144 ClassImp(StTofrMatchMaker)
00145 
00146 // Define the  static constants:
00147 
00148 const Int_t StTofrMatchMaker::mDAQOVERFLOW = 255;
00149 const Int_t StTofrMatchMaker::mNTOFP       = 41;
00150 const Int_t StTofrMatchMaker::mNPVPD       = 6;
00151 const Int_t StTofrMatchMaker::mNTOFR       = 120;
00152 const Int_t StTofrMatchMaker::mNTOFR5      = 192;
00153 
00154 const Int_t StTofrMatchMaker::mNTOF    = 192;    // 192 for tof in Run 8++
00155 const Int_t StTofrMatchMaker::mNModule = 32;  // 32 for tofr5++ 
00156 const Int_t StTofrMatchMaker::mNCell   = 6;
00157 const Int_t StTofrMatchMaker::mNVPD    = 19;    //
00158 
00159 const Int_t StTofrMatchMaker::mEastVpdTrayId = 122;
00160 const Int_t StTofrMatchMaker::mWestVpdTrayId = 121;
00161 
00162 const Int_t StTofrMatchMaker::mNValidTrays_Run3 = 1;
00163 const Int_t StTofrMatchMaker::mNValidTrays_Run4 = 1;
00164 const Int_t StTofrMatchMaker::mNValidTrays_Run5 = 1;
00165 const Int_t StTofrMatchMaker::mNValidTrays_Run6 = 0;
00166 const Int_t StTofrMatchMaker::mNValidTrays_Run7 = 0;
00167 const Int_t StTofrMatchMaker::mNValidTrays_Run8 = 5;
00168     
00169 const Int_t StTofrMatchMaker::mTdigBoard  = 10;
00170 const Int_t StTofrMatchMaker::mTdcOnBoard = 4;
00171 const Int_t StTofrMatchMaker::mTdcChannel = 1024;
00172 
00173 //---------------------------------------------------------------------------
00174 StTofrMatchMaker::StTofrMatchMaker(const Char_t *name): StMaker(name)
00175  , mTofrAdc(mNTOFR,0)
00176  , mTofrTdc(mNTOFR,0)
00177  , mPvpdAdc(mNPVPD,0)
00178  , mPvpdAdcLoRes( mNPVPD,0)
00179  , mPvpdTdc(mNPVPD,0)
00180  , mStrobeTdcMin(mNPVPD,0)
00181  , mStrobeTdcMax(mNPVPD,0)
00182  , mPedTOFr(mNTOFR,0)
00183     //year 5
00184  , mPvpdToT(mNPVPD,0)
00185  , mTofr5Tdc(mNTOFR5,0)
00186  , mTofr5ToT(mNTOFR5,0) // ToT as adc
00187  , mHitCorr(mNValidTrays_Run8,(TH2D*)0)
00188  , mHitCorrModule(mNValidTrays_Run8,(TH2D*)0)
00189 {
00190 
00191   // set default values
00192   mEventCounter = 0;
00193   mAcceptedEventCounter = 0;
00194   mTofEventCounter = 0;
00195   mTofStrobeEventCounter = 0;
00196   mAcceptAndStrobe = 0;
00197   mAcceptAndBeam = 0;
00198 
00199   mTofrGeom = 0;
00200   mDaqMap = 0;
00201   mTofINLCorr = 0;
00202 
00203   mWidthPad = 3.45;
00204 
00205   setValidAdcRange(1,1200);
00206   setValidTdcRange(1,2047);
00207   setOuterTrackGeometry();
00208   setMinHitsPerTrack(15);
00209   setMinFitPointsPerTrack(15);
00210   setMinFitPointsOverMax(0.52);
00211   setMaxDCA(9999.);
00212 
00213   setCreateHistoFlag(kFALSE);
00214   setHistoFileName("tofana.root");
00215   setCreateTreeFlag(kFALSE);
00216   setSaveGeometry(kFALSE);
00217   doPrintMemoryInfo = kFALSE;
00218   doPrintCpuInfo    = kFALSE;
00219 }
00220 
00221 StTofrMatchMaker::~StTofrMatchMaker(){ /* nope */}
00222 
00223 //void StTofrMatchMaker::Clear(Option_t *opt){StMaker::Clear();}
00224 
00225 //---------------------------------------------------------------------------
00226 Int_t StTofrMatchMaker::Init(){
00227   gMessMgr->Info("StTofrMatchMaker -- initializing ...","OS");
00228   LOG_INFO  << "Valid TDC range: " << mMinValidTdc << " " << mMaxValidTdc << endm;
00229   LOG_INFO  << "Valid ADC range: " << mMinValidAdc << " " << mMaxValidAdc << endm;
00230   LOG_INFO  << "Minimum hits per track: " << mMinHitsPerTrack << endm;
00231   LOG_INFO  << "Minimum fitpoints per track: " << mMinFitPointsPerTrack << endm;
00232   LOG_INFO  << "Maximum DCA: " << mMaxDCA << endm;
00233   if (!mOuterTrackGeometry)
00234     gMessMgr->Warning("Warning: using standard trackgeometry()","OS");
00235   
00236   // m_Mode can be set by SetMode() method
00237   if(m_Mode) {
00238 //    setHistoFileName("tofana.root");
00239   } else {
00240     setHistoFileName("");
00241   }
00242 
00243   if (mHisto){
00244     bookHistograms();
00245     LOG_INFO  << "Histograms are booked" << endm;
00246     if (mHistoFileName!="")
00247       LOG_INFO  << "Histograms will be stored in " << mHistoFileName.c_str() << endm;
00248   }
00249 
00250   // reset event counters
00251   mEventCounter = 0;
00252   mAcceptedEventCounter = 0;
00253   mTofEventCounter = 0;
00254   mTofStrobeEventCounter = 0;
00255   mAcceptAndStrobe = 0;
00256   mAcceptAndBeam = 0;
00257   
00258   return kStOK;
00259 }
00260 
00261 
00262 //---------------------------------------------------------------------------
00263 Int_t StTofrMatchMaker::InitRun(Int_t runnumber){
00264 
00265   // determine TOF configuration from run#
00266   mYear2 = (runnumber<4000000);
00267   mYear3 = (runnumber>4000000&&runnumber<5000000);
00268   mYear4 = (runnumber>5000000&&runnumber<6000000);
00269   mYear5 = (runnumber>6000000&&runnumber<7000000);
00270   mYear8 = (runnumber>9000000&&runnumber<10000000);
00271   mYearX = (runnumber>10000000);
00272 
00273   if (mYearX){
00274     Error(":InitRun","Wrong BFC configuration for run %d. Use StBTofHitMaker for Run9+ data.",runnumber);
00275     return 0;
00276   }
00277 
00278   gMessMgr->Info("StTofrMatchMaker -- Initializing TofGeometry (InitRun)","OS");
00280   // TOFr geometry initializtion -- from GEANT geometry directly
00281   //                               need St_geant_Maker be loaded before
00283   mTofrGeom = new StTofrGeometry("tofrGeom","tofrGeom in MatchMaker");
00284   if(!mTofrGeom->IsInitDone()) {
00285     gMessMgr->Info("TofrGemetry initialization..." ,"OS");
00286     TVolume *starHall = (TVolume *)GetDataSet("HALL");
00287     mTofrGeom->Init(starHall);
00288   }
00289   // other makers can get this geometry
00290   if(mGeometrySave) {
00291     if(TDataSet *geom = GetDataSet("tofrGeometry")) delete geom;
00292     AddConst(new TObjectSet("tofrGeometry",mTofrGeom)); 
00293   }
00294 
00295   if(Debug()) {
00296     gMessMgr->Info(" Test the TofrGeometry","OS");
00297     // Test
00298     TVolumeView* topNode = mTofrGeom->GetTopNode();
00299     if(!topNode) {
00300       gMessMgr->Warning("Warning:  no Top node!","OS");
00301       return kStWarn;
00302     }
00303     mTofrGeom->Print();
00304     
00305     topNode->ls(9);
00306     //  TShape *topshape = topNode->GetShape();
00307     LOG_INFO  << " # of trays = " << topNode->GetListSize() << endm;
00308     TList *list = topNode->Nodes();
00309     Int_t ibtoh =0;
00310     TVolumeView *sectorVolume = 0;
00311     for(Int_t i=0;i<list->GetSize();i++) {
00312       sectorVolume = dynamic_cast<TVolumeView*> (list->At(i));
00313       if ( i>=60 ) ibtoh = 1;
00314       LOG_INFO  << " test sector size = " <<sectorVolume->GetListSize() << endm;
00315       StTofrGeomTray *tray = new StTofrGeomTray(ibtoh, sectorVolume, topNode);
00316       tray->Print();
00317       TVolumeView *trayVolume = tray->GetfView();
00318       TList *list1 = trayVolume->Nodes();
00319       LOG_INFO  << "   # of modules in tray " << tray->Index() << " = " << trayVolume->GetListSize() << endm;
00320       if (!list1 ) continue;
00321       TVolumeView *sensorVolume = 0;
00322       for(Int_t j=0;j<list1->GetSize();j++) {
00323         sensorVolume = dynamic_cast<TVolumeView*> (list1->At(j));
00324         StTofrGeomSensor *sensor = new StTofrGeomSensor(sensorVolume, topNode);
00325         sensor->Print();
00326         delete sensor;  
00327       }
00328       delete tray;
00329     }
00330   }
00331 
00333   // TOFr Daq map initialization -- load from StTofUtil
00335   mDaqMap = new StTofrDaqMap();
00336   mTofINLCorr = new StTofINLCorr();
00337   if(mYear2||mYear3||mYear4) {
00338     if(mYear3) mDaqMap->setNValidTrays(mNValidTrays_Run3);
00339     if(mYear4) mDaqMap->setNValidTrays(mNValidTrays_Run4);
00340     mDaqMap->init(this);
00341     //  AddConst(new TObjectSet("tofrDaqMap",mDaqMap));
00342     LOG_INFO  << " Initialize Daq map for run 2,3,4 ... " << endm;
00343     
00344     if(Debug()) {
00345       for(Int_t i=0;i<mNTOFR;i++) {
00346         LOG_INFO  << i << " -- adc=" << mDaqMap->DaqChan2ADCChan(i) << " -- tdc=" << mDaqMap->DaqChan2TDCChan(i) << endm;
00347         IntVec map = mDaqMap->DaqChan2Cell(i);
00348         LOG_INFO  << "    tray=" << map[0] << "  module=" << map[1] << "  cell=" << map[2] << endm;
00349       }
00350     }
00351     
00353       // Load configuration parameters from dbase
00354       //    need "[shell] setenv Calibrations_tof reconV0"
00356       gMessMgr->Info("    -- retrieving run parameters from Calibrations_tof","OS");
00357       TDataSet *mDbDataSet = GetDataBase("Calibrations/tof");
00358       if (!mDbDataSet){
00359         gMessMgr->Error("unable to get TOF run parameters","OS");
00360         //    assert(mDbDataSet);
00361         return kStErr;
00362       }
00363       St_pvpdStrobeDef* pvpdStrobeDef = static_cast<St_pvpdStrobeDef*>(mDbDataSet->Find("pvpdStrobeDef"));
00364       if (!pvpdStrobeDef){
00365         gMessMgr->Error("unable to find pVPD strobe def parameters","OS");
00366         //    assert(pvpdStrobeDef);
00367         return kStErr;
00368       }
00369       pvpdStrobeDef_st *strobeDef = static_cast<pvpdStrobeDef_st*>(pvpdStrobeDef->GetArray());
00370       Int_t numRows = pvpdStrobeDef->GetNRows();
00371       if (mNPVPD != numRows) gMessMgr->Warning("#tubes inconsistency in dbase");
00372       for (Int_t i=0;i<mNPVPD;i++){
00373         Int_t ii = strobeDef[i].id - 1;
00374         mStrobeTdcMin[ii] = strobeDef[i].strobeTdcMin;
00375         mStrobeTdcMax[ii] = strobeDef[i].strobeTdcMax;
00376         if (Debug())
00377           LOG_INFO  << "tube " << strobeDef[i].id << "  min:"<< strobeDef[i].strobeTdcMin
00378                                   <<" max:"<< strobeDef[i].strobeTdcMax<< endm;
00379       }
00380       
00381       //TOFr pedestals
00382       St_tofPedestal* tofPed = static_cast<St_tofPedestal*>(mDbDataSet->Find("tofPedestal"));
00383       if (!tofPed){
00384         gMessMgr->Error("unable to find TOF pedestal table","OS");
00385         return kStErr;
00386         //    assert(tofPed);
00387       }
00388       tofPedestal_st *pedestal = static_cast<tofPedestal_st*>(tofPed->GetArray());
00389       Int_t entries = (Int_t)(pedestal[0].entries);
00390       if (mNTOFR>entries) gMessMgr->Warning("#TOFr channels inconsistency in dbase");
00391       for (Int_t i=0;i<entries;i++) {
00392         Int_t ii = pedestal[0].daqChannel[i];
00393         Int_t adcChan = pedestal[0].adcChan[i];
00394         if(adcChan>=60) 
00395           mPedTOFr[ii] = (Double_t)(pedestal[0].adcPedestal[i]);
00396         if (Debug())
00397           LOG_INFO  << "daqChannel" << ii << " ped:" << pedestal[0].adcPedestal[i] <<endm;
00398       }
00399       
00400       /*   only test
00401       //TOFr configuration
00402       St_tofTrayConfig* toftrayconfig = static_cast<St_tofTrayConfig*>(mDbDataSet->Find("tofTrayConfig"));
00403       if(!toftrayconfig) {
00404       gMessMgr->Error("unable to find TOF tray configuration parameters","OS");
00405       //    assert(toftrayconfig);
00406       return kStErr;    
00407       }
00408       tofTrayConfig_st *toftrayconf = static_cast<tofTrayConfig_st*>(toftrayconfig->GetArray());
00409       entries = toftrayconf[0].entries;
00410       for(Int_t i=0;i<entries;i++) {
00411       Int_t iTray = toftrayconf[0].iTray[i];
00412       Int_t nModules = toftrayconf[0].nModules[i];
00413       if(Debug())
00414       LOG_INFO  <<  "   " << iTray << " tray has " << nModules << " modules " << endm;
00415       }
00416       */
00417       
00418   } else if(mYear5) {
00419     mDaqMap->setNValidTrays(mNValidTrays_Run5);
00420     mDaqMap->initFromDbaseY5(this);
00421     LOG_INFO  << " Initialize Daq map for run 5 ... " << endm;
00422     
00423     if(Debug()) {
00424       //    if(0) {
00425       for(Int_t i=0;i<mNTOFR5;i++) {
00426         IntVec map = mDaqMap->Tofr5TDCChan2Cell(i);
00427         LOG_INFO  << "InitRun():i="<<i<<" tray=" << map[0] << "  module=" << map[1] << "  cell=" << map[2] << endm;
00428       }
00429     }
00430     
00431     /* move INL calibration to StTofCalibMaker
00432     // read INL table
00433     LOG_INFO  << "read in INL table" << endm;
00434     TDataSet *mDbTOFDataSet = this->GetDataBase("Calibrations/tof");
00435     if(!mDbTOFDataSet) {
00436       gMessMgr->Error("unable to access Calibrations TOF parameters","OS");
00437       //    assert(mDbTOFDataSet);
00438       return kStErr;
00439     }
00440     
00441     St_tofr5INLtable* tofr5INLtable = static_cast<St_tofr5INLtable*>(mDbTOFDataSet->Find("tofr5INLtable"));
00442     if(!tofr5INLtable) {
00443       gMessMgr->Error("unable to get tofr5 INL table parameters","OS");
00444       //    assert(tofr5INLtable);
00445       return kStErr;
00446     }
00447     tofr5INLtable_st* inltable = static_cast<tofr5INLtable_st*>(tofr5INLtable->GetArray());
00448     Int_t numRows = tofr5INLtable->GetNRows();
00449     
00450     LOG_INFO  << " Number of rows read in: " << numRows << endm;
00451     
00452     Char_t *boardName;
00453     Short_t boardId;
00454     Short_t tdcId;
00455     Float_t INLcorr[1024];
00456     
00457     for (Int_t i=0;i<numRows;i++) {
00458       boardName = "";
00459       boardId = -1;
00460       tdcId = -1;
00461       for(int j=0;j<mTdcChannel;j++) {
00462         INLcorr[j] = 0.0;
00463       }
00464       
00465       boardName = (Char_t *)(inltable[i].boardID);
00466       boardId = inltable[i].boardNumber;
00467       tdcId = inltable[i].TDCID;
00468       if(Debug())
00469         LOG_INFO  << " name = " << boardName << " bId = " << boardId << " tdcId = " << tdcId << endm;
00470       for(int j=0;j<mTdcChannel;j++) {
00471         INLcorr[j] = inltable[i].INLcorrection[j];
00472         if(Debug()&&j%100==0) LOG_INFO  << " j=" << j << " inlcorr=" << INLcorr[j] << endm;
00473         mINLtable[boardId][tdcId][j] = INLcorr[j];
00474       }
00475             
00476     }
00477     */
00478 
00479   } else if(mYear8) {
00480     mDaqMap->setNValidTrays(mNValidTrays_Run8);
00481     mDaqMap->initFromDbaseGeneral(this);
00482     LOG_INFO  << " Initialize Daq map for run 8 ... " << endm;
00483 
00484     mTofINLCorr->setNValidTrays(mNValidTrays_Run8);
00485     mTofINLCorr->initFromDbase(this);
00486     LOG_INFO  << " Initialize INL table for run 8 ... " << endm;
00487 
00488   }
00489   
00490   return kStOK;
00491 }
00492 
00493 //----------------------------------------------------------------------------
00494 Int_t StTofrMatchMaker::FinishRun(Int_t runnumber){
00495 
00496   LOG_INFO  << "StTofrMatchMaker -- cleaning up geometry (FinishRun)" << endm;
00497   if (mTofrGeom) delete mTofrGeom;
00498   mTofrGeom=0;
00499 
00500   if(mDaqMap) delete mDaqMap;
00501   mDaqMap = 0;
00502 
00503   if(mTofINLCorr) delete mTofINLCorr;
00504   mTofINLCorr = 0;
00505 
00506   return kStOK;
00507 }
00508 
00509 
00510 //---------------------------------------------------------------------------
00511 Int_t StTofrMatchMaker::Finish(){
00512 
00513   LOG_INFO  << "StTofrMatchMaker -----  RUN SUMMARY ----- (Finish)\n"
00514        << "\tProcessed "  << mEventCounter << " events."
00515        << " Accepted  "   << mAcceptedEventCounter << " events."
00516        << " Rejected  "   << mEventCounter - mAcceptedEventCounter << " events\n"
00517        << "\tTOF events " << mTofEventCounter
00518        << ". Beam "       << mTofEventCounter - mTofStrobeEventCounter
00519        << "  Strobe "     << mTofStrobeEventCounter
00520        << "\n\t Accept & Strobe " << mAcceptAndStrobe << " events\n"
00521        << "\t Accept & Beam   "   << mAcceptAndBeam   << " events" << endm;
00522   
00523   //if (mHisto) writeHistogramsToFile();
00524   if (mHistoFileName!="") writeHistogramsToFile();
00525   return kStOK;
00526 }
00527 
00528 
00529 //---------------------------------------------------------------------------
00530 Int_t StTofrMatchMaker::Make(){
00531   gMessMgr->Info("StTofrMatchMaker -- welcome","OS");
00532 
00533   Int_t iret = kStOK;
00534   if(mYear2||mYear3||mYear4) {
00535     iret = processEventYear2to4();
00536   } else if(mYear5) {
00537     iret = processEventYear5();
00538   } else if(mYear8) {
00539     iret = processEventYear8();
00540   }
00541   return iret;
00542 }
00543 
00544 //---------------------------------------------------------------------------
00545 Int_t StTofrMatchMaker::processEventYear2to4(){
00546   if(mHisto) mEventCounterHisto->Fill(0);
00547   // event selection ...
00548   mEvent = (StEvent *) GetInputDS("StEvent");
00549   if (!validEvent(mEvent)){
00550     gMessMgr->Info("StTofrMatchMaker -- nothing to do ... bye-bye","OS");
00551     return kStOK;
00552   }
00553 
00554   if (mYear2) {
00555     gMessMgr->Info("StTofrMatchMaker -- no TOFr in Year2!","OS");
00556     return kStOK;
00557   }
00558 
00559   // timing & memory info -only when requested-
00560   StTimer timer;
00561   if (doPrintCpuInfo) timer.start();
00562   if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
00563 
00564   //.........................................................................
00565   // check for tofCollection and fill local copy with ADC and TDC data
00566   StTofCollection *theTof = mEvent->tofCollection();
00567   getTofData(theTof);
00568 
00569   //.........................................................................
00570   // update Tofr cell counters
00571   Float_t sumAdcTofr(0.); Int_t nAdcTofr(0), nTdcTofr(0), nAdcTdcTofr(0);
00572   for (Int_t i=0;i<mNTOFR;i++){
00573     sumAdcTofr += mTofrAdc[i];
00574     bool tdcValid = validTdc(mTofrTdc[i]);
00575     bool adcValid = validAdc(mTofrAdc[i]);
00576     if (tdcValid) nTdcTofr++;
00577     if (adcValid) nAdcTofr++;
00578     if (adcValid && tdcValid) nAdcTdcTofr++;
00579   }
00580   if (Debug())
00581     LOG_INFO  << " Tofr #Adc:" << nAdcTofr << "   #Tdc:" << nTdcTofr << endm;
00582 
00583   //..................................
00584   // update pVPD tubes counters
00585   Float_t sumAdcPvpd=0; Int_t nAdcPvpd=0, nTdcPvpd=0;
00586   for (Int_t i=0;i<mNPVPD;i++){
00587     sumAdcPvpd += mPvpdAdc[i];
00588     bool tdcValid = validTdc(mPvpdTdc[i]);
00589     bool adcValid = validAdc(mPvpdAdc[i]);
00590     if (tdcValid) nTdcPvpd++;
00591     if (adcValid) nAdcPvpd++;
00592   }
00593   if (Debug())
00594     LOG_INFO  << " pVPD #Adc:" << nAdcPvpd << "   #Tdc:" << nTdcPvpd << endm;
00595 
00596 
00597     
00598   // number of primary tracks
00599   //  (note: different meaning in event.root and richtof.root)
00600   Int_t refmult(0);
00601   bool richTofMuDST = (mEvent->summary()->numberOfExoticTracks() == -999);
00602   if (richTofMuDST)
00603     refmult = mEvent->summary()->numberOfGoodTracks();
00604   else
00605     refmult = uncorrectedNumberOfPrimaries(*mEvent);
00606 
00607   if (Debug()){
00608     LOG_INFO  << " #Tracks           :"      << mEvent->summary()->numberOfTracks()
00609          << "\n #goodPrimaryTracks:"    << mEvent->summary()->numberOfGoodPrimaryTracks()
00610          << "\n #uncorr.prim.tracks  :" << refmult << endm;
00611     if (!richTofMuDST)
00612       LOG_INFO  << " #goodTracks (global):"  << mEvent->summary()->numberOfGoodTracks() << endm;
00613   }
00614 
00615   //.........................................................................
00616   // A. build vector of candidate cells with valid ADC signals 
00617   //  idVector validCellIdVec;
00618   tofCellHitVector daqCellsHitVec;
00619   //  daqCellsHitVec.clear();
00620   idVector validModuleVec;
00621 
00622   for (Int_t i=0;i<mNTOFR;i++){
00623     Float_t rawAdc = mTofrAdc[i];
00624     Float_t rawTdc = mTofrTdc[i];
00625     if (mHisto) mDaqOccupancy->Fill(i);
00626 
00627     if (validAdc(rawAdc) && validTdc(rawTdc) && rawAdc>mPedTOFr[i] ){
00628       IntVec map = mDaqMap->DaqChan2Cell(i);
00629       Int_t IDtray = map[0];
00630       Int_t IDmodule = map[1];
00631       Int_t IDcell = map[2];
00632       
00633       StructCellHit aDaqCellHit;
00634       aDaqCellHit.channel = i;
00635       aDaqCellHit.tray = IDtray;
00636       aDaqCellHit.module = IDmodule;
00637       aDaqCellHit.cell = IDcell;
00638       daqCellsHitVec.push_back(aDaqCellHit);
00639       // additional valid number configuration
00640       int id = IDtray*100+IDmodule;
00641       bool ifind = kFALSE;
00642       for(size_t iv=0;iv<validModuleVec.size();iv++) {
00643         if(id==validModuleVec[iv]) {
00644           ifind = kTRUE;
00645           break;
00646         }
00647       }
00648       if(!ifind) validModuleVec.push_back(id);
00649 
00650       if (mHisto) {
00651         mDaqOccupancyValid->Fill(i);
00652         Int_t adcchan = mDaqMap->DaqChan2ADCChan(i);
00653         Int_t tdcchan = mDaqMap->DaqChan2TDCChan(i);
00654         mADCTDCCorelation->Fill(tdcchan, adcchan);
00655       }
00656       if(Debug()) {
00657         LOG_INFO  <<"A: daqId=" << i << "  rawAdc= " << rawAdc << " rawTdc="<< rawTdc <<endm;
00658       }
00659     }
00660   }
00661   // end of Sect.A
00662   if(Debug()) {
00663     LOG_INFO  << "    total # of cells = " << daqCellsHitVec.size() << endm;
00664     for(size_t iv = 0;iv<validModuleVec.size();iv++) {
00665       LOG_INFO  << " module # " << validModuleVec[iv] << " Valid! " << endm;
00666     }
00667   }
00668   if(mHisto) {
00669     mCellsMultInEvent->Fill(daqCellsHitVec.size());
00670     if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
00671   }
00672   if(!daqCellsHitVec.size()) return kStOK;
00673 
00674   //.........................................................................
00675   // B. loop over global tracks and determine all cell-track matches
00676   //
00677   tofCellHitVector allCellsHitVec;
00678   //  allCellsHitVec.clear();
00679   StructCellHit cellHit;
00680 
00681   StSPtrVecTrackNode& nodes = mEvent->trackNodes();
00682   Int_t nAllTracks=0;
00683   for (unsigned int iNode=0; iNode<nodes.size(); iNode++){
00684     tofCellHitVector cellHitVec;
00685     //    cellHitVec.clear();
00686     StTrack *theTrack = nodes[iNode]->track(global);
00687 
00688     // make sure we have a track, a miniDST might have removed it...
00689     if (validTrack(theTrack)){
00690       nAllTracks++;
00691       StPhysicalHelixD theHelix = trackGeometry(theTrack)->helix();
00692 
00693       IntVec projTrayVec;
00694       if(!mTofrGeom->projTrayVector(theHelix, projTrayVec)) continue;
00695 
00696       IntVec idVec;
00697       DoubleVec pathVec;
00698       PointVec  crossVec;
00699 
00700 //       idVec.clear();
00701 //       pathVec.clear();
00702 //       crossVec.clear();
00703 
00704       Int_t ncells = 0;
00705       //      if(mTofrGeom->HelixCrossCellIds(theHelix,idVec,pathVec,crossVec) ) {
00706       if(mTofrGeom->HelixCrossCellIds(theHelix, validModuleVec, projTrayVec, idVec, pathVec, crossVec)) {
00707         Int_t cells = idVec.size();
00708         for (Int_t i=0; i<cells; i++) {
00709             Int_t icell,imodule,itray;
00710             Double_t local[3],global[3];
00711             for(Int_t i2=0;i2<3;i2++){
00712                  local[i2]=0;
00713             }
00714             global[0]=crossVec[i].x();
00715             global[1]=crossVec[i].y();
00716             global[2]=crossVec[i].z();
00717             mTofrGeom->DecodeCellId(idVec[i], icell, imodule, itray);
00718             //      LOG_INFO << " decode " << idVec[i] << "  to tray#" << itray << " module#" << imodule << " cell#" << icell << endm;
00719             StTofrGeomSensor* sensor = 
00720                   mTofrGeom->GetGeomSensor(imodule,itray);
00721             if(!sensor) {
00722               gMessMgr->Warning("","OS") << " No sensitive module in the projection??? -- Something weird!!! " << endm;
00723               continue;
00724             }
00725             sensor->Master2Local(&global[0],&local[0]);
00726             icell = sensor->FindCellIndex(local);
00727             //      StThreeVectorD glo=sensor->GetCenterPosition();
00728             StThreeVectorD glo(global[0], global[1], global[2]);
00729             StThreeVectorD hitPos(local[0], local[1], local[2]);
00730             delete sensor;
00731             if (local[2]<=3.4&&local[2]>=-2.7) {
00732               Int_t Iarray = mDaqMap->Cell2DaqChan(itray, imodule, icell);
00733               if(Iarray>=mDAQOVERFLOW) continue;
00734               ncells++;
00735               cellHit.channel = Iarray;
00736               cellHit.tray = itray;
00737               cellHit.module = imodule;
00738               cellHit.cell = icell;
00739               cellHit.trackIdVec.push_back(iNode);
00740               cellHit.hitPosition = glo;        // global position
00741               cellHit.zhit = (Float_t)hitPos.z();
00742               cellHit.yhit = (Float_t)hitPos.y();
00743               cellHitVec.push_back(cellHit);
00744               allCellsHitVec.push_back(cellHit);
00745               if(mHisto) {
00746                 mDaqOccupancyProj->Fill(Iarray);
00747                 mHitsPosition->Fill(hitPos.y(), hitPos.z());
00748               }
00749               
00750               if(Debug()) {
00751                 LOG_INFO  <<"B: nodeid=" << iNode << "  projected in " << " tray="<< itray << " module="<<imodule<<" cell="<<icell<<endm;
00752                 LOG_INFO  <<"   hit position " << hitPos << endm;
00753               }
00754             }
00755         } // for (Int_t i=0...)
00756       } // endif(helixcross...)
00757       if(ncells>0&&mHisto) mHitsMultPerTrack->Fill(ncells);
00758 
00759     } // if(ValidTrack).. 
00760   } // loop over nodes
00761   if(Debug())
00762     LOG_INFO  << "B:  matched/available/total #tracknodes: " <<allCellsHitVec.size() << "/" <<nAllTracks << "/" << nodes.size() << endm;
00763   if(mHisto) {
00764     mHitsMultInEvent->Fill(allCellsHitVec.size());
00765     if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
00766   }
00767   // end of Sect.B
00768   
00769   //.........................................................................
00770   // C. Match find Neighbours -- identify crosstalk
00771   //
00772   tofCellHitVector matchHitCellsVec;
00773   //  matchHitCellsVec.clear();
00774 
00775   tofCellHitVectorIter daqIter = daqCellsHitVec.begin();
00776   for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
00777     tofCellHitVectorIter proIter = allCellsHitVec.begin();
00778     for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
00779       if( (daqIter->tray==proIter->tray)&& 
00780           (daqIter->module==proIter->module) &&
00781           ( ( (proIter->cell==6)&&((proIter->cell==daqIter->cell) ||
00782                                    (proIter->cell==daqIter->cell+1)) )
00783             || ( (proIter->cell==1)&&((proIter->cell==daqIter->cell) ||
00784                                       (proIter->cell==daqIter->cell-1)) )
00785             || ( (proIter->cell>=2&&proIter->cell<=6) &&
00786                  ( (proIter->cell==daqIter->cell) ||
00787                    (proIter->cell==daqIter->cell-1) ||
00788                    (proIter->cell==daqIter->cell+1) ) ) ) ) {
00789         cellHit.channel = daqIter->channel;
00790         cellHit.tray = daqIter->tray;
00791         cellHit.module = daqIter->module;
00792         cellHit.cell = daqIter->cell;
00793         cellHit.hitPosition = proIter->hitPosition;
00794         cellHit.trackIdVec = proIter->trackIdVec;
00795         cellHit.zhit = proIter->zhit;
00796         cellHit.yhit = proIter->yhit;
00797         matchHitCellsVec.push_back(cellHit);
00798       }
00799     }
00800   } //end {sec. C}
00801   if(Debug()) {
00802     LOG_INFO  << "C: before/after: " << allCellsHitVec.size() << "/" << matchHitCellsVec.size() << endm;
00803   }
00804   if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
00805 
00806   //.........................................................................
00807   // D. sort hit vectors  and deal with (discard) cells matched by multiple tracks
00808   //
00809   Int_t nSingleHitCells(0);
00810   Int_t nMultiHitsCells(0);
00811 
00812   tofCellHitVector singleHitCellsVec;
00813   tofCellHitVector multiHitsCellsVec;
00814 //   singleHitCellsVec.clear();
00815 //   multiHitsCellsVec.clear();
00816 
00817   tofCellHitVector tempVec = matchHitCellsVec;
00818   tofCellHitVector erasedVec = tempVec;
00819   while (tempVec.size() != 0) {
00820     Int_t nTracks = 0;
00821     idVector trackIdVec;
00822 
00823     tofCellHitVectorIter tempIter=tempVec.begin();
00824     tofCellHitVectorIter erasedIter=erasedVec.begin();
00825     while(erasedIter!= erasedVec.end()) {
00826       if(tempIter->tray == erasedIter->tray &&
00827          tempIter->module == erasedIter->module &&
00828          tempIter->cell == erasedIter->cell) {
00829         nTracks++;
00830         trackIdVec.push_back(erasedIter->trackIdVec.back());  // merge
00831         erasedVec.erase(erasedIter);
00832         erasedIter--;
00833       }
00834       erasedIter++;
00835     }
00836 
00837     cellHit.channel = tempIter->channel;
00838     cellHit.cell = tempIter->cell;
00839     cellHit.module = tempIter->module;
00840     cellHit.tray = tempIter->tray;
00841     cellHit.hitPosition = tempIter->hitPosition;
00842     cellHit.trackIdVec = trackIdVec;
00843     cellHit.zhit = tempIter->zhit;
00844     cellHit.yhit = tempIter->yhit;
00845 
00846     Float_t ycenter = (tempIter->cell-1-2.5)*mWidthPad;
00847     Float_t dy = tempIter->yhit - ycenter;
00848     Float_t dz = tempIter->zhit;
00849 
00850     if(mHisto) {
00851       mTracksPerCellMatch1->Fill(trackIdVec.size());
00852       mDaqOccupancyMatch1->Fill(tempIter->channel);
00853       mDeltaHitMatch1->Fill(dy, dz);
00854     }
00855 
00856     if (nTracks==1){
00857       nSingleHitCells++;      
00858       singleHitCellsVec.push_back(cellHit);
00859     } else if (nTracks>1){
00860       nMultiHitsCells++;
00861       multiHitsCellsVec.push_back(cellHit);
00862       // for multiple hit cells either discard (yes) or
00863       // find the most likely candidate.
00864     } else {
00865       LOG_INFO  << "D: no tracks extrapolate to matched cell ... should not happen!" << endm;
00866     }
00867     
00868     if (Debug()) {
00869       LOG_INFO  << "D: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
00870       idVectorIter ij=trackIdVec.begin();
00871       while (ij != trackIdVec.end()) { LOG_INFO  << " " << *ij; ij++; }
00872       LOG_INFO  <<endm;
00873     }
00874     
00875     tempVec = erasedVec;
00876   }
00877   if(Debug())
00878     LOG_INFO  << "D: before/after: " << matchHitCellsVec.size() << "/" << singleHitCellsVec.size() << endm;
00879   //end of Sect.C
00880   if(mHisto) {
00881     mCellsPerEventMatch1->Fill(singleHitCellsVec.size()+multiHitsCellsVec.size());
00882     if(singleHitCellsVec.size()) mEventCounterHisto->Fill(9);
00883   } 
00884 
00885   //.........................................................................
00886   // E. sort and deal singleHitCellsVector for multiple cells associated to single tracks
00887   //
00888   tofCellHitVector FinalMatchedCellsVec;
00889   //  FinalMatchedCellsVec.clear();
00890   tempVec = singleHitCellsVec;
00891   if(mHisto) {
00892     mCellsPerEventMatch2->Fill(tempVec.size());
00893     for(unsigned int ii=0;ii<tempVec.size();ii++) {
00894       mTracksPerCellMatch2->Fill(tempVec[ii].trackIdVec.size());
00895       mDaqOccupancyMatch2->Fill(tempVec[ii].channel);
00896       Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
00897       Float_t dy = tempVec[ii].yhit-ycenter;
00898       Float_t dz = tempVec[ii].zhit;
00899       mDeltaHitMatch2->Fill(dy, dz);
00900     }
00901   }
00902 
00903   erasedVec = tempVec;
00904   while (tempVec.size() != 0) {
00905     StructCellHit cellHit;
00906     Int_t nCells = 0;
00907     idVector vTrackId;
00908     vector<StThreeVectorD> vPosition;
00909     vector<Int_t> vchannel, vtray, vmodule, vcell;
00910     vector<Float_t> vzhit, vyhit;
00911 
00912     tofCellHitVectorIter tempIter=tempVec.begin();
00913     tofCellHitVectorIter erasedIter=erasedVec.begin();
00914     while(erasedIter!= erasedVec.end()) {
00915       if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
00916         nCells++;
00917         vchannel.push_back(erasedIter->channel);
00918         vtray.push_back(erasedIter->tray);
00919         vmodule.push_back(erasedIter->module);
00920         vcell.push_back(erasedIter->cell);
00921         vPosition.push_back(erasedIter->hitPosition);
00922         vTrackId.push_back(erasedIter->trackIdVec.back());
00923         vzhit.push_back(erasedIter->zhit);
00924         vyhit.push_back(erasedIter->yhit);
00925 
00926         erasedVec.erase(erasedIter);
00927         erasedIter--;
00928       }
00929       erasedIter++;
00930     }
00931 
00932     
00933     if (nCells==1){
00934       // for singly hit cell, copy data in singleHitCellsVec
00935       cellHit.channel = vchannel[0];
00936       cellHit.tray = vtray[0];
00937       cellHit.module = vmodule[0];
00938       cellHit.cell = vcell[0];
00939       cellHit.trackIdVec.push_back(vTrackId[0]);
00940       cellHit.hitPosition = vPosition[0];
00941       cellHit.matchFlag = 0; 
00942       cellHit.zhit = vzhit[0];
00943       cellHit.yhit = vyhit[0];
00944 
00945       FinalMatchedCellsVec.push_back(cellHit);
00946 
00947       // debugging output
00948       if (Debug()) {
00949         LOG_INFO  << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
00950         idVectorIter ij=vTrackId.begin();
00951         while (ij != vTrackId.end()) { LOG_INFO  << " " << *ij; ij++; }
00952         LOG_INFO  <<endm;
00953       }
00954     }
00955     else if (nCells>1){   // for multiple hit cells  find the most likely candidate.
00956       Int_t thiscandidate(-99);
00957       Int_t thisMatchFlag(0);
00958 
00959       // 1. sort on adc weight
00960       Float_t weight(0);
00961       vector<Int_t> weightCandidates;
00962       thisMatchFlag = 1;
00963       if (Debug()) LOG_INFO  << "E: find ... weight ";
00964       for (Int_t i=0;i<nCells;i++){
00965         if(vchannel[i]>=0&&vchannel[i]<mNTOFR) {
00966           Float_t adcwt = mTofrAdc[vchannel[i]]-mPedTOFr[vchannel[i]];
00967           if(adcwt>weight) {
00968             weight = adcwt;
00969             weightCandidates.clear();
00970             weightCandidates.push_back(i);
00971           } else if (adcwt == weight) {
00972             weightCandidates.push_back(i);
00973           }
00974         }
00975       }
00976       if (weightCandidates.size()==1){
00977         thiscandidate = weightCandidates[0];
00978         Int_t daqId = vchannel[thiscandidate];
00979         if (Debug()) LOG_INFO  << "candidate =" << daqId << endm;
00980       }
00981 
00982       // 2. if still undecided check on hitposition
00983       if (weightCandidates.size()>1){
00984         Float_t ss(99.);
00985         vector<Int_t> ssCandidates;
00986         thisMatchFlag = 2;
00987         if (Debug()) LOG_INFO  << " ss ";
00988         for (unsigned int i=0;i<weightCandidates.size();i++){
00989           Int_t ii=weightCandidates[i];   
00990           //      Int_t daqId = vchannel[ii];
00991           StThreeVectorD ahit = vPosition[ii];
00992           Float_t yy = ahit.y();
00993           Float_t ycell = (vcell[ii]-1-2.5)*mWidthPad;
00994           Float_t ll = fabs(yy-ycell);
00995           if(ll<ss) {
00996             ss = ll;
00997             ssCandidates.clear();
00998             ssCandidates.push_back(ii);
00999           }else if  (ll==ss)
01000             ssCandidates.push_back(ii);   
01001         }
01002         if (ssCandidates.size()==1){
01003           thiscandidate = ssCandidates[0];
01004           Int_t daqId = vchannel[thiscandidate];
01005           if (Debug()) LOG_INFO  << "candidate =" << daqId << endm;
01006         }
01007 
01008       }
01009 
01010       if (thiscandidate>=0) {
01011         cellHit.channel = vchannel[thiscandidate];
01012         cellHit.tray = vtray[thiscandidate];
01013         cellHit.module = vmodule[thiscandidate];
01014         cellHit.cell = vcell[thiscandidate];
01015         cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
01016         cellHit.hitPosition = vPosition[thiscandidate];
01017         cellHit.matchFlag = thisMatchFlag;
01018         cellHit.zhit = vzhit[thiscandidate];
01019         cellHit.yhit = vyhit[thiscandidate];
01020 
01021         FinalMatchedCellsVec.push_back(cellHit);
01022         
01023         // debugging output
01024         if (Debug()) {
01025           LOG_INFO  << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:" << vTrackId[thiscandidate] << endm;
01026         }
01027       }
01028 
01029     } else {
01030       LOG_INFO  << "E: no cells belong to this track ... should not happen!" << endm;
01031     }
01032 
01033     tempVec = erasedVec;
01034   }
01035 
01036   LOG_INFO  << "E: before/after: " << singleHitCellsVec.size() << "/" << FinalMatchedCellsVec.size() << endm;
01037   // end of Sect.E
01038 
01039   //.........................................................................
01040   // F. perform further selection and
01041   //    fill valid track histograms, ntuples and CellCollection
01042   //
01043   tempVec.clear();
01044   tempVec = FinalMatchedCellsVec;
01045   if(mHisto) {
01046     if(FinalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
01047     mCellsPerEventMatch3->Fill(tempVec.size());
01048     for(unsigned int ii=0;ii<tempVec.size();ii++) {
01049       mTracksPerCellMatch3->Fill(tempVec[ii].trackIdVec.size());
01050       mDaqOccupancyMatch3->Fill(tempVec[ii].channel);
01051       Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
01052       Float_t dy = tempVec[ii].yhit - ycenter;
01053       Float_t dz = tempVec[ii].zhit;
01054       mDeltaHitMatch3->Fill(dy, dz);
01055     }
01056   }
01057 
01058   StTofCellCollection *mCellCollection =  new StTofCellCollection;
01059   Int_t nValidSingleHitCells(0), nValidSinglePrimHitCells(0);
01060 
01061   for (size_t ii=0; ii < FinalMatchedCellsVec.size(); ii++){
01062     Int_t daqId = FinalMatchedCellsVec[ii].channel;
01063     Int_t jj = daqId;
01064     Int_t tray = FinalMatchedCellsVec[ii].tray;
01065     Int_t module = FinalMatchedCellsVec[ii].module;
01066     Int_t cell = FinalMatchedCellsVec[ii].cell;
01067 
01068     Float_t ycenter = (cell-1-2.5)*mWidthPad;
01069     Float_t dy = FinalMatchedCellsVec[ii].yhit - ycenter;
01070     if (FinalMatchedCellsVec[ii].trackIdVec.size()!=1)
01071       LOG_INFO  << "F: WHAT!?!  mult.matched cell in single cell list " << daqId << endm;
01072 
01073     // 1. fill valid single track AND valid tdc histograms
01074     if (validTdc(mTofrTdc[jj])) nValidSingleHitCells++;
01075 
01076     // get track-id from cell hit vector
01077     unsigned int trackNode = FinalMatchedCellsVec[ii].trackIdVec[0];
01078     StTrack *theTrack = nodes[trackNode]->track(primary);
01079     StTrack *globalTrack = nodes[trackNode]->track(global);
01080 
01081     // 2. continue only if the (primary) track exists
01082     if (validTofTrack(theTrack) && fabs(dy)<1.9 ){
01083       nValidSinglePrimHitCells++;
01084 
01085       //--- store number of hits per track
01086       Int_t nHitsPerTrack = theTrack->topologyMap().numberOfHits(kTpcId);
01087           
01088       // select the apropriate track geometry
01089       StTrackGeometry *theTrackGeometry = trackGeometry(theTrack);
01090 
01091       //--- get momentum from track
01092       const StThreeVectorF momentum = theTrackGeometry->momentum();
01093             
01094       //--- calculate flight path
01095       Double_t pathLength = tofPathLength(&mEvent->primaryVertex()->position(),
01096                                         &FinalMatchedCellsVec[ii].hitPosition,
01097                                         theTrackGeometry->helix().curvature());
01098 
01099       //--- calculate local hit position on cell based first, last and middle plane
01100       //    (middle plane is the average of first and last plane, which is mathematically
01101       //     the same as CellHitVec.hitPosition ... )
01102 //       StThreeVectorD *pInnerLayer, *pOuterLayer;
01103 //       pInnerLayer =  FinalMatchedCellsVec[ii].layerHitPositions.begin();
01104 //       pOuterLayer =  FinalMatchedCellsVec[ii].layerHitPositions.end() - 1;
01105       
01106       //--- dig out from the dedx and rich pid traits
01107       Float_t dedx(0.), cherang(0);
01108       Int_t dedx_np(0), cherang_nph(0);
01109       StSPtrVecTrackPidTraits& traits = theTrack->pidTraits();
01110       for (unsigned int it=0;it<traits.size();it++){
01111         if (traits[it]->detector() == kTpcId){
01112           StDedxPidTraits *pid = dynamic_cast<StDedxPidTraits*>(traits[it]);
01113           if (pid && pid->method() ==kTruncatedMeanId){
01114             dedx    = pid->mean();
01115             dedx_np =  pid->numberOfPoints();
01116           }
01117         } else if  (traits[it]->detector() == kRichId){
01118           StRichPidTraits *pid = dynamic_cast<StRichPidTraits*>(traits[it]);
01119           if (pid){ 
01120             StRichSpectra* pidinfo = pid->getRichSpectra();
01121             if (pidinfo && pidinfo->getCherenkovPhotons()>2){
01122               cherang     = pidinfo->getCherenkovAngle();
01123               cherang_nph = pidinfo->getCherenkovPhotons();
01124             }
01125           }
01126         }
01127       }
01128       
01129       //--- calculate local hit position on cell based on average hitposition
01130       //      Float_t localHitPos = mTofGeom->cellHitPosition(&allMatchedCellsVec[ii].hitPosition);
01131       
01132       // Fill TOF Cell Collection
01133       StTofCell *tofCell = new StTofCell(tray, module, cell, daqId, (Int_t)(mTofrAdc[jj]-mPedTOFr[jj]),(Int_t)mTofrTdc[jj],theTrack,FinalMatchedCellsVec[ii].zhit,FinalMatchedCellsVec[ii].matchFlag,FinalMatchedCellsVec[ii].hitPosition);
01134       mCellCollection->push_back(tofCell);
01135       
01136       // dump debug data
01137       if (Debug()){
01138         LOG_INFO  << "F: itray=" << tray << " imodule=" << module << " icell=" << cell << "\tnodeid:";
01139         idVectorIter ij=FinalMatchedCellsVec[ii].trackIdVec.begin();
01140         while (ij != FinalMatchedCellsVec[ii].trackIdVec.end()) { LOG_INFO  << " " << *ij; ij++; }
01141         LOG_INFO  << "\tR=" << 1/(theTrackGeometry->helix().curvature())
01142              << "\tpT=" << momentum.perp() << "\tp=" << momentum.mag()
01143              << "\thits="<< nHitsPerTrack << "\ts="<< pathLength
01144              << "\t#fitp=" <<theTrack->fitTraits().numberOfFitPoints(kTpcId)
01145           //         << "\t#trkp=" <<theTrack->detectorInfo()->numberOfPoints(kTpcId)
01146              << " \tdedx=" << dedx
01147              << " \tdca="<< globalTrack->geometry()->helix().distance(mEvent->primaryVertex()->position())<<" and "<<theTrackGeometry->helix().distance(mEvent->primaryVertex()->position());
01148         if (cherang!=0) LOG_INFO  << " \trich="<< cherang << " (" << cherang_nph << ")";
01149         LOG_INFO  << endm;
01150       }
01151 
01152     } // track exists 
01153   }
01154   
01155   storeMatchData(mCellCollection,theTof);
01156   delete mCellCollection;
01157   
01158   //check StEvent collections --
01159   if (theTof->dataPresent())
01160     LOG_INFO  << " TofCollection: raw data container present" << endm;
01161   if (theTof->cellsPresent()){
01162     LOG_INFO  << " TofCollection: cell container present."<<endm;
01163     if (Debug()){
01164       StSPtrVecTofCell& tmpCellTofVec = theTof->tofCells();
01165       for (size_t i = 0; i < tmpCellTofVec.size(); i++) {
01166         StTofCell* p = tmpCellTofVec[i];
01167         LOG_INFO  << p->trayIndex() << " " << p->moduleIndex() << " " << p->cellIndex() << " " << p->adc() << " " << p->tdc() << " " << p->associatedTrack() << " " << p->matchFlag() << " " << p->position() << endm;
01168       }
01169     }
01170   }
01171   //-- end check
01172   
01173   LOG_INFO  << "F: before/after" << FinalMatchedCellsVec.size() << "/" <<nValidSinglePrimHitCells << endm;
01174  // end of Sect.F
01175 
01176   LOG_INFO  << "#(cell tracks): " << allCellsHitVec.size()
01177        << " #(hit cells): " << FinalMatchedCellsVec.size()
01178        << " #cells (valid tdc): " << nTdcTofr
01179        << "\n#(single hits): " << nSingleHitCells 
01180        << " #(single valid hits): " << nValidSingleHitCells
01181        << " #(single prim valid hits): " << nValidSinglePrimHitCells
01182        << endm;
01183 
01184 
01185 
01186   if (doPrintMemoryInfo) {
01187         StMemoryInfo::instance()->snapshot();
01188         StMemoryInfo::instance()->print();
01189   }
01190   if (doPrintCpuInfo) {
01191     timer.stop();
01192     LOG_INFO  << "CPU time for StTofrMatchMaker::Make(): "
01193          << timer.elapsedTime() << " sec\n" << endm;
01194   }
01195 
01196   LOG_INFO  << "StTofrMatchMaker -- bye-bye" << endm;
01197   return kStOK;
01198 }
01199 
01200 //---------------------------------------------------------------------------
01201 Int_t StTofrMatchMaker::processEventYear5(){
01202   // leave as empty now
01203 
01204   if(mHisto) mEventCounterHisto->Fill(0);
01205   // event selection ...
01206   mEvent = (StEvent *) GetInputDS("StEvent");
01207   if (!validEvent(mEvent)){
01208     gMessMgr->Info("StTofrMatchMaker -- nothing to do ... bye-bye","OS");
01209     return kStOK;
01210   }
01211 
01212   // timing & memory info -only when requested-
01213   StTimer timer;
01214   if (doPrintCpuInfo) timer.start();
01215   if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
01216 
01217   //.........................................................................
01218   // check for tofCollection and fill local copy with ADC and TDC data
01219   StTofCollection *theTof = mEvent->tofCollection();
01220   //  getTofData(theTof);
01221   // year5
01222   mSortTofRawData = new StSortTofRawData(theTof);
01223 
01224   IntVec validchannel = mSortTofRawData->GetValidChannel();
01225   // Any Hits in TOF+PVPD ?
01226   if(validchannel.size()<=0) {
01227     LOG_INFO  << " No hits in TOF or pVPD! " << endm;
01228     // return kStOK;
01229   } else {
01230     LOG_INFO  << " Number of TOF+pVPD fired hits: " << validchannel.size() << endm;
01231   }
01232   if(Debug()) {
01233     for(size_t iv=0;iv<validchannel.size();iv++) {
01234       //if(validchannel[iv]<0||validchannel[iv]>=mNTOFR5) continue; // skip pvpd
01235       LOG_INFO  << " channel = " << validchannel[iv]<<endm;
01236       IntVec leTdc = mSortTofRawData->GetLeadingTdc(validchannel[iv]);
01237       IntVec teTdc = mSortTofRawData->GetTrailingTdc(validchannel[iv]);
01238       for(size_t iv1=0;iv1<leTdc.size();iv1++) {
01239         LOG_INFO  << " leading Tdc = " << leTdc[iv1]<<endm;
01240       }
01241       for(size_t iv2=0;iv2<teTdc.size();iv2++) {
01242         LOG_INFO  << " trailing Tdc = " << teTdc[iv2] << endm;
01243       }
01244     }
01245   }   // end mdebug 
01246     
01247   // number of primary tracks
01248   //  (note: different meaning in event.root and richtof.root)
01249   Int_t refmult(0);
01250   refmult = uncorrectedNumberOfPrimaries(*mEvent);
01251   
01252   if (Debug()){
01253     LOG_INFO  << " #Tracks           :"      << mEvent->summary()->numberOfTracks()
01254          << "\n #goodPrimaryTracks:"    << mEvent->summary()->numberOfGoodPrimaryTracks()
01255          << "\n #uncorr.prim.tracks  :" << refmult << endm;
01256   }
01257 
01258   //.........................................................................
01259   // A. build vector of candidate cells with valid ADC signals 
01260   //  idVector validCellIdVec;
01261   tofCellHitVector daqCellsHitVec;
01262   //  daqCellsHitVec.clear();
01263   idVector validModuleVec;
01264 
01265   for(size_t ich=0;ich<validchannel.size();ich++) {
01266     int ichan = validchannel[ich];
01267     if(ichan<0||ichan>=mNTOFR5) continue; // tray
01268     IntVec map = mDaqMap->Tofr5TDCChan2Cell(ichan);
01269     if(map.size()!=3) continue;
01270 
01271     Int_t IDtray = map[0];
01272     Int_t IDmodule = map[1];
01273     Int_t IDcell = map[2];
01274     
01275     StructCellHit aDaqCellHit;
01276     aDaqCellHit.channel = ichan;
01277     aDaqCellHit.tray = IDtray;
01278     aDaqCellHit.module = IDmodule;
01279     aDaqCellHit.cell = IDcell;
01280     daqCellsHitVec.push_back(aDaqCellHit);
01281 
01282     // additional valid number configuration
01283     int id = IDtray*100+IDmodule;
01284     bool ifind = kFALSE;
01285     for(size_t iv=0;iv<validModuleVec.size();iv++) {
01286       if(id==validModuleVec[iv]) {
01287         ifind = kTRUE;
01288         break;
01289       }
01290     }
01291     if(!ifind) validModuleVec.push_back(id);
01292   }
01293 
01294   // end of Sect.A
01295   if(Debug()) {
01296     LOG_INFO  << "    total # of cells = " << daqCellsHitVec.size() << endm;
01297     for(size_t iv = 0;iv<validModuleVec.size();iv++) {
01298       LOG_INFO  << " module # " << validModuleVec[iv] << " Valid! " << endm;
01299     }
01300   }
01301   if(mHisto) {
01302     mCellsMultInEvent->Fill(daqCellsHitVec.size());
01303     if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
01304   }
01305   if(!daqCellsHitVec.size()) return kStOK;
01306 
01307   //.........................................................................
01308   // B. loop over global tracks and determine all cell-track matches
01309   //
01310   tofCellHitVector allCellsHitVec;
01311   //  allCellsHitVec.clear();
01312   StructCellHit cellHit;
01313 
01314   StSPtrVecTrackNode& nodes = mEvent->trackNodes();
01315   Int_t nAllTracks=0;
01316   for (unsigned int iNode=0; iNode<nodes.size(); iNode++){
01317     tofCellHitVector cellHitVec;
01318     //    cellHitVec.clear();
01319     StTrack *theTrack = nodes[iNode]->track(global);
01320 
01321     // make sure we have a track, a miniDST might have removed it...
01322     if (validTrack(theTrack)){
01323       nAllTracks++;
01324       StPhysicalHelixD theHelix = trackGeometry(theTrack)->helix();
01325 
01326       IntVec projTrayVec;
01327       if(!mTofrGeom->projTrayVector(theHelix, projTrayVec)) continue;
01328 
01329       IntVec idVec;
01330       DoubleVec pathVec;
01331       PointVec  crossVec;
01332 
01333 //       idVec.clear();
01334 //       pathVec.clear();
01335 //       crossVec.clear();
01336 
01337       Int_t ncells = 0;
01338       //      if(mTofrGeom->HelixCrossCellIds(theHelix,idVec,pathVec,crossVec) ) {
01339       if(mTofrGeom->HelixCrossCellIds(theHelix, validModuleVec, projTrayVec, idVec, pathVec, crossVec)) {
01340         Int_t cells = idVec.size();
01341         for (Int_t i=0; i<cells; i++) {
01342             Int_t icell,imodule,itray;
01343             Double_t local[3],global[3];
01344             for(Int_t i2=0;i2<3;i2++){
01345                  local[i2]=0;
01346             }
01347             global[0]=crossVec[i].x();
01348             global[1]=crossVec[i].y();
01349             global[2]=crossVec[i].z();
01350             mTofrGeom->DecodeCellId(idVec[i], icell, imodule, itray);
01351             LOG_INFO << " decode " << idVec[i] << "  to tray#" << itray << " module#" << imodule << " cell#" << icell << endm;
01352             StTofrGeomSensor* sensor = 
01353                   mTofrGeom->GetGeomSensor(imodule,itray);
01354             if(!sensor) {
01355               gMessMgr->Warning("","OS") << " No sensitive module in the projection??? -- Something weird!!! " << endm;
01356               continue;
01357             }
01358             sensor->Master2Local(&global[0],&local[0]);
01359             icell = sensor->FindCellIndex(local);
01360             //      StThreeVectorD glo=sensor->GetCenterPosition();
01361             StThreeVectorD glo(global[0], global[1], global[2]);
01362             StThreeVectorD hitPos(local[0], local[1], local[2]);
01363             delete sensor;
01364             //      if (local[2]<=3.4&&local[2]>=-2.7) {
01365               Int_t Iarray = mDaqMap->Tofr5Cell2TDCChan(itray, imodule, icell);
01366               if(Iarray>=mDAQOVERFLOW||Iarray<0) continue;
01367               ncells++;
01368               cellHit.channel = Iarray;
01369               cellHit.tray = itray;
01370               cellHit.module = imodule;
01371               cellHit.cell = icell;
01372               cellHit.trackIdVec.push_back(iNode);
01373               cellHit.hitPosition = glo;        // global position
01374               cellHit.zhit = (Float_t)hitPos.z();
01375               cellHit.yhit = (Float_t)hitPos.y();
01376               cellHitVec.push_back(cellHit);
01377               allCellsHitVec.push_back(cellHit);
01378               if(mHisto) {
01379                 mDaqOccupancyProj->Fill(Iarray);
01380                 mHitsPosition->Fill(hitPos.y(), hitPos.z());
01381               }
01382               
01383               if(Debug()) {
01384                 LOG_INFO  <<"B: nodeid=" << iNode << "  projected in " << " tray="<< itray << " module="<<imodule<<" cell="<<icell<<endm;
01385                 LOG_INFO  <<"   hit position " << hitPos << endm;
01386               }
01387               //            }
01388         } // for (Int_t i=0...)
01389       } // endif(helixcross...)
01390       if(ncells>0&&mHisto) mHitsMultPerTrack->Fill(ncells);
01391 
01392     } // if(ValidTrack).. 
01393   } // loop over nodes
01394   if(Debug())
01395     LOG_INFO  << "B:  matched/available/total #tracknodes: " <<allCellsHitVec.size() << "/" <<nAllTracks << "/" << nodes.size() << endm;
01396   if(mHisto) {
01397     mHitsMultInEvent->Fill(allCellsHitVec.size());
01398     if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
01399   }
01400   // end of Sect.B
01401   
01402   //.........................................................................
01403   // C. Match find Neighbours -- identify crosstalk
01404   //
01405   tofCellHitVector matchHitCellsVec;
01406   //  matchHitCellsVec.clear();
01407 
01408   tofCellHitVectorIter daqIter = daqCellsHitVec.begin();
01409   for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
01410     tofCellHitVectorIter proIter = allCellsHitVec.begin();
01411     for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
01412       if( (daqIter->tray==proIter->tray)&& 
01413           (daqIter->module==proIter->module) &&
01414           ( ( (proIter->cell==6)&&((proIter->cell==daqIter->cell) ||
01415                                    (proIter->cell==daqIter->cell+1)) )
01416             || ( (proIter->cell==1)&&((proIter->cell==daqIter->cell) ||
01417                                       (proIter->cell==daqIter->cell-1)) )
01418             || ( (proIter->cell>=2&&proIter->cell<=6) &&
01419                  ( (proIter->cell==daqIter->cell) ||
01420                    (proIter->cell==daqIter->cell-1) ||
01421                    (proIter->cell==daqIter->cell+1) ) ) ) ) {
01422         cellHit.channel = daqIter->channel;
01423         cellHit.tray = daqIter->tray;
01424         cellHit.module = daqIter->module;
01425         cellHit.cell = daqIter->cell;
01426         cellHit.hitPosition = proIter->hitPosition;
01427         cellHit.trackIdVec = proIter->trackIdVec;
01428         cellHit.zhit = proIter->zhit;
01429         cellHit.yhit = proIter->yhit;
01430         matchHitCellsVec.push_back(cellHit);
01431       }
01432     }
01433   } //end {sec. C}
01434   if(Debug()) {
01435     LOG_INFO  << "C: before/after: " << allCellsHitVec.size() << "/" << matchHitCellsVec.size() << endm;
01436   }
01437   if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
01438 
01439   //.........................................................................
01440   // D. sort hit vectors  and deal with (discard) cells matched by multiple tracks
01441   //
01442   Int_t nSingleHitCells(0);
01443   Int_t nMultiHitsCells(0);
01444 
01445   tofCellHitVector singleHitCellsVec;
01446   tofCellHitVector multiHitsCellsVec;
01447 //   singleHitCellsVec.clear();
01448 //   multiHitsCellsVec.clear();
01449 
01450   tofCellHitVector tempVec = matchHitCellsVec;
01451   tofCellHitVector erasedVec = tempVec;
01452   while (tempVec.size() != 0) {
01453     Int_t nTracks = 0;
01454     idVector trackIdVec;
01455 
01456     tofCellHitVectorIter tempIter=tempVec.begin();
01457     tofCellHitVectorIter erasedIter=erasedVec.begin();
01458     while(erasedIter!= erasedVec.end()) {
01459       if(tempIter->tray == erasedIter->tray &&
01460          tempIter->module == erasedIter->module &&
01461          tempIter->cell == erasedIter->cell) {
01462         nTracks++;
01463         trackIdVec.push_back(erasedIter->trackIdVec.back());  // merge
01464         erasedVec.erase(erasedIter);
01465         erasedIter--;
01466       }
01467       erasedIter++;
01468     }
01469 
01470     cellHit.channel = tempIter->channel;
01471     cellHit.cell = tempIter->cell;
01472     cellHit.module = tempIter->module;
01473     cellHit.tray = tempIter->tray;
01474     cellHit.hitPosition = tempIter->hitPosition;
01475     cellHit.trackIdVec = trackIdVec;
01476     cellHit.zhit = tempIter->zhit;
01477     cellHit.yhit = tempIter->yhit;
01478 
01479     Float_t ycenter = (tempIter->cell-1-2.5)*mWidthPad;
01480     Float_t dy = tempIter->yhit - ycenter;
01481     Float_t dz = tempIter->zhit;
01482 
01483     if(mHisto) {
01484       mTracksPerCellMatch1->Fill(trackIdVec.size());
01485       mDaqOccupancyMatch1->Fill(tempIter->channel);
01486       mDeltaHitMatch1->Fill(dy, dz);
01487     }
01488 
01489     if (nTracks==1){
01490       nSingleHitCells++;      
01491       singleHitCellsVec.push_back(cellHit);
01492     } else if (nTracks>1){
01493       nMultiHitsCells++;
01494       multiHitsCellsVec.push_back(cellHit);
01495       // for multiple hit cells either discard (yes) or
01496       // find the most likely candidate.
01497     } else {
01498       LOG_INFO  << "D: no tracks extrapolate to matched cell ... should not happen!" << endm;
01499     }
01500     
01501     if (Debug()) {
01502       LOG_INFO  << "D: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
01503       idVectorIter ij=trackIdVec.begin();
01504       while (ij != trackIdVec.end()) { LOG_INFO  << " " << *ij; ij++; }
01505       LOG_INFO  <<endm;
01506     }
01507     
01508     tempVec = erasedVec;
01509   }
01510   if(Debug())
01511     LOG_INFO  << "D: before/after: " << matchHitCellsVec.size() << "/" << singleHitCellsVec.size() << endm;
01512   //end of Sect.C
01513   if(mHisto) {
01514     mCellsPerEventMatch1->Fill(singleHitCellsVec.size()+multiHitsCellsVec.size());
01515     if(singleHitCellsVec.size()) mEventCounterHisto->Fill(9);
01516   } 
01517 
01518   //.........................................................................
01519   // E. sort and deal singleHitCellsVector for multiple cells associated to single tracks
01520   //
01521   tofCellHitVector FinalMatchedCellsVec;
01522   //  FinalMatchedCellsVec.clear();
01523   tempVec = singleHitCellsVec;
01524   if(mHisto) {
01525     mCellsPerEventMatch2->Fill(tempVec.size());
01526     for(unsigned int ii=0;ii<tempVec.size();ii++) {
01527       mTracksPerCellMatch2->Fill(tempVec[ii].trackIdVec.size());
01528       mDaqOccupancyMatch2->Fill(tempVec[ii].channel);
01529       Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
01530       Float_t dy = tempVec[ii].yhit-ycenter;
01531       Float_t dz = tempVec[ii].zhit;
01532       mDeltaHitMatch2->Fill(dy, dz);
01533     }
01534   }
01535 
01536   erasedVec = tempVec;
01537   while (tempVec.size() != 0) {
01538     StructCellHit cellHit;
01539     Int_t nCells = 0;
01540     idVector vTrackId;
01541     vector<StThreeVectorD> vPosition;
01542     vector<Int_t> vchannel, vtray, vmodule, vcell;
01543     vector<Float_t> vzhit, vyhit;
01544 
01545     tofCellHitVectorIter tempIter=tempVec.begin();
01546     tofCellHitVectorIter erasedIter=erasedVec.begin();
01547     while(erasedIter!= erasedVec.end()) {
01548       if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
01549         nCells++;
01550         vchannel.push_back(erasedIter->channel);
01551         vtray.push_back(erasedIter->tray);
01552         vmodule.push_back(erasedIter->module);
01553         vcell.push_back(erasedIter->cell);
01554         vPosition.push_back(erasedIter->hitPosition);
01555         vTrackId.push_back(erasedIter->trackIdVec.back());
01556         vzhit.push_back(erasedIter->zhit);
01557         vyhit.push_back(erasedIter->yhit);
01558 
01559         erasedVec.erase(erasedIter);
01560         erasedIter--;
01561       }
01562       erasedIter++;
01563     }
01564 
01565     if (nCells==1){
01566       // for singly hit cell, copy data in singleHitCellsVec
01567       cellHit.channel = vchannel[0];
01568       cellHit.tray = vtray[0];
01569       cellHit.module = vmodule[0];
01570       cellHit.cell = vcell[0];
01571       cellHit.trackIdVec.push_back(vTrackId[0]);
01572       cellHit.hitPosition = vPosition[0];
01573       cellHit.matchFlag = 0; 
01574       cellHit.zhit = vzhit[0];
01575       cellHit.yhit = vyhit[0];
01576 
01577       FinalMatchedCellsVec.push_back(cellHit);
01578 
01579       // debugging output
01580       if (Debug()) {
01581         LOG_INFO  << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
01582         idVectorIter ij=vTrackId.begin();
01583         while (ij != vTrackId.end()) { LOG_INFO  << " " << *ij; ij++; }
01584         LOG_INFO  <<endm;
01585       }
01586     }
01587     else if (nCells>1){   // for multiple hit cells  find the most likely candidate.
01588       Int_t thiscandidate(-99);
01589       Int_t thisMatchFlag(0);
01590 
01591       // sort on hitposition
01592       Float_t ss(99.);
01593       vector<Int_t> ssCandidates;
01594       thisMatchFlag = 2;
01595       if (Debug()) LOG_INFO  << " ss " << endm;
01596       for (Int_t i=0;i<nCells;i++){
01597         Float_t yy = vyhit[i];
01598         Float_t ycell = (vcell[i]-1-2.5)*mWidthPad;
01599         Float_t ll = fabs(yy-ycell);
01600         if(ll<ss) {
01601           ss = ll;
01602           ssCandidates.clear();
01603           ssCandidates.push_back(i);
01604         }else if  (ll==ss)
01605           ssCandidates.push_back(i);      
01606       }
01607       if (ssCandidates.size()==1){
01608         thiscandidate = ssCandidates[0];
01609         Int_t daqId = vchannel[thiscandidate];
01610         if (Debug()) LOG_INFO  << "candidate =" << daqId << endm;
01611       }
01612       
01613 
01614       if (thiscandidate>=0) {
01615         cellHit.channel = vchannel[thiscandidate];
01616         cellHit.tray = vtray[thiscandidate];
01617         cellHit.module = vmodule[thiscandidate];
01618         cellHit.cell = vcell[thiscandidate];
01619         cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
01620         cellHit.hitPosition = vPosition[thiscandidate];
01621         cellHit.matchFlag = thisMatchFlag;
01622         cellHit.zhit = vzhit[thiscandidate];
01623         cellHit.yhit = vyhit[thiscandidate];
01624 
01625         FinalMatchedCellsVec.push_back(cellHit);
01626         
01627         // debugging output
01628         if (Debug()) {
01629         LOG_INFO  << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:" << vTrackId[thiscandidate] << endm;
01630         }
01631       }
01632 
01633     } else {
01634       LOG_INFO  << "E: no cells belong to this track ... should not happen!" << endm;
01635     }
01636 
01637     tempVec = erasedVec;
01638   }
01639 
01640   LOG_INFO  << "E: before/after: " << singleHitCellsVec.size() << "/" << FinalMatchedCellsVec.size() << endm;
01641   // end of Sect.E
01642 
01643   //.........................................................................
01644   // F. perform further selection and
01645   //    fill valid track histograms, ntuples and CellCollection
01646   //
01647   tempVec.clear();
01648   tempVec = FinalMatchedCellsVec;
01649   if(mHisto) {
01650     if(FinalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
01651     mCellsPerEventMatch3->Fill(tempVec.size());
01652     for(unsigned int ii=0;ii<tempVec.size();ii++) {
01653       mTracksPerCellMatch3->Fill(tempVec[ii].trackIdVec.size());
01654       mDaqOccupancyMatch3->Fill(tempVec[ii].channel);
01655       Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
01656       Float_t dy = tempVec[ii].yhit - ycenter;
01657       Float_t dz = tempVec[ii].zhit;
01658       mDeltaHitMatch3->Fill(dy, dz);
01659     }
01660   }
01661 
01662   StTofCellCollection *mCellCollection =  new StTofCellCollection;
01663   Int_t nValidSingleHitCells(0), nValidSinglePrimHitCells(0);
01664 
01665   for (size_t ii=0; ii < FinalMatchedCellsVec.size(); ii++){
01666     Int_t daqId = FinalMatchedCellsVec[ii].channel;
01667     Int_t jj = daqId;
01668     Int_t tray = FinalMatchedCellsVec[ii].tray;
01669     Int_t module = FinalMatchedCellsVec[ii].module;
01670     Int_t cell = FinalMatchedCellsVec[ii].cell;
01671 
01672     Float_t ycenter = (cell-1-2.5)*mWidthPad;
01673     Float_t dy = FinalMatchedCellsVec[ii].yhit - ycenter;
01674     if (FinalMatchedCellsVec[ii].trackIdVec.size()!=1)
01675       LOG_INFO  << "F: WHAT!?!  mult.matched cell in single cell list " << daqId << endm;
01676 
01677 
01678     /* move INL correction to calibration maker
01679     // Read in Leading and Trailing edge TDC, apply on INL correction
01680     int tmptdc = (mSortTofRawData->GetLeadingTdc(jj))[0];
01681     // do inl correction
01682     int bin = int(tmptdc)&0x03ff;
01683     float tmptdc_f = tmptdc + GetINLcorr(4, jj, bin);
01684     float letime  = tmptdc_f*25./1024; // ns
01685 
01686     tmptdc=(mSortTofRawData->GetTrailingTdc(jj))[0];    
01687     // do inl correction
01688     bin = int(tmptdc)&0x0ff;
01689     tmptdc_f = tmptdc + GetINLcorr(5, jj, bin);
01690     float tetime  = tmptdc_f*100./1024; // ns
01691     float tot = tetime - letime; // ns
01692     */
01693 
01694     int letdc = (mSortTofRawData->GetLeadingTdc(jj))[0];
01695     int tetdc=(mSortTofRawData->GetTrailingTdc(jj))[0];
01696 
01697     // get track-id from cell hit vector
01698     unsigned int trackNode = FinalMatchedCellsVec[ii].trackIdVec[0];
01699     StTrack *theTrack = nodes[trackNode]->track(primary);
01700     StTrack *globalTrack = nodes[trackNode]->track(global);
01701 
01702     // 2. continue only if the (primary) track exists
01703     if (validTofTrack(theTrack) && fabs(dy)<1.9 ){
01704       nValidSinglePrimHitCells++;
01705 
01706       //--- store number of hits per track
01707       Int_t nHitsPerTrack = theTrack->topologyMap().numberOfHits(kTpcId);
01708           
01709       // select the apropriate track geometry
01710       StTrackGeometry *theTrackGeometry = trackGeometry(theTrack);
01711 
01712       //--- get momentum from track
01713       const StThreeVectorF momentum = theTrackGeometry->momentum();
01714             
01715       //--- calculate flight path
01716       Double_t pathLength = tofPathLength(&mEvent->primaryVertex()->position(),
01717                                         &FinalMatchedCellsVec[ii].hitPosition,
01718                                         theTrackGeometry->helix().curvature());
01719 
01720       //--- calculate local hit position on cell based first, last and middle plane
01721       //    (middle plane is the average of first and last plane, which is mathematically
01722       //     the same as CellHitVec.hitPosition ... )
01723 //       StThreeVectorD *pInnerLayer, *pOuterLayer;
01724 //       pInnerLayer =  FinalMatchedCellsVec[ii].layerHitPositions.begin();
01725 //       pOuterLayer =  FinalMatchedCellsVec[ii].layerHitPositions.end() - 1;
01726       
01727       //--- dig out from the dedx and rich pid traits
01728       Float_t dedx(0.), cherang(0);
01729       Int_t dedx_np(0), cherang_nph(0);
01730       StSPtrVecTrackPidTraits& traits = theTrack->pidTraits();
01731       for (unsigned int it=0;it<traits.size();it++){
01732         if (traits[it]->detector() == kTpcId){
01733           StDedxPidTraits *pid = dynamic_cast<StDedxPidTraits*>(traits[it]);
01734           if (pid && pid->method() ==kTruncatedMeanId){
01735             dedx    = pid->mean();
01736             dedx_np =  pid->numberOfPoints();
01737           }
01738         } else if  (traits[it]->detector() == kRichId){
01739           StRichPidTraits *pid = dynamic_cast<StRichPidTraits*>(traits[it]);
01740           if (pid){ 
01741             StRichSpectra* pidinfo = pid->getRichSpectra();
01742             if (pidinfo && pidinfo->getCherenkovPhotons()>2){
01743               cherang     = pidinfo->getCherenkovAngle();
01744               cherang_nph = pidinfo->getCherenkovPhotons();
01745             }
01746           }
01747         }
01748       }
01749       
01750       //--- calculate local hit position on cell based on average hitposition
01751       //      Float_t localHitPos = mTofGeom->cellHitPosition(&allMatchedCellsVec[ii].hitPosition);
01752       
01753       // Fill TOF Cell Collection
01754       // year5, raw letdc and tetdc stored as "tdc" and "adc" in StTofCell
01755       // 
01756       //      StTofCell *tofCell = new StTofCell(tray, module, cell, daqId, (Int_t)(tot*1000.),(Int_t)(letime*1000.),theTrack,FinalMatchedCellsVec[ii].zhit,FinalMatchedCellsVec[ii].matchFlag,FinalMatchedCellsVec[ii].hitPosition);
01757       StTofCell *tofCell = new StTofCell(tray, module, cell, daqId, tetdc,letdc,theTrack,FinalMatchedCellsVec[ii].zhit,FinalMatchedCellsVec[ii].matchFlag,FinalMatchedCellsVec[ii].hitPosition);
01758       mCellCollection->push_back(tofCell);
01759       
01760       // dump debug data
01761       if (Debug()){
01762         LOG_INFO  << "F: itray=" << tray << " imodule=" << module << " icell=" << cell << "\tnodeid:";
01763         idVectorIter ij=FinalMatchedCellsVec[ii].trackIdVec.begin();
01764         while (ij != FinalMatchedCellsVec[ii].trackIdVec.end()) { LOG_INFO  << " " << *ij; ij++; }
01765         LOG_INFO  << "\tR=" << 1/(theTrackGeometry->helix().curvature())
01766              << "\tpT=" << momentum.perp() << "\tp=" << momentum.mag()
01767              << "\thits="<< nHitsPerTrack << "\ts="<< pathLength
01768              << "\t#fitp=" <<theTrack->fitTraits().numberOfFitPoints(kTpcId)
01769           //         << "\t#trkp=" <<theTrack->detectorInfo()->numberOfPoints(kTpcId)
01770              << " \tdedx=" << dedx
01771              << " \tdca="<< globalTrack->geometry()->helix().distance(mEvent->primaryVertex()->position())<<" and "<<theTrackGeometry->helix().distance(mEvent->primaryVertex()->position());
01772         if (cherang!=0) LOG_INFO  << " \trich="<< cherang << " (" << cherang_nph << ")";
01773         LOG_INFO  << endm;
01774       }
01775 
01776     } // track exists 
01777   } // end final matched cells
01778   
01779   storeMatchData(mCellCollection,theTof);
01780   delete mCellCollection;
01781 
01782   delete mSortTofRawData;
01783   mSortTofRawData = 0;
01784   
01785   //check StEvent collections --
01786   if (theTof->dataPresent())
01787     LOG_INFO  << " TofCollection: raw data container present" << endm;
01788   if (theTof->cellsPresent()){
01789     LOG_INFO  << " TofCollection: cell container present."<<endm;
01790     if (Debug()){
01791       StSPtrVecTofCell& tmpCellTofVec = theTof->tofCells();
01792       LOG_INFO  << " # of matched cells " << tmpCellTofVec.size() << endm;
01793       for (size_t i = 0; i < tmpCellTofVec.size(); i++) {
01794         StTofCell* p = tmpCellTofVec[i];
01795         LOG_INFO  << p->trayIndex() << " " << p->moduleIndex() << " " << p->cellIndex() << " " << p->trailingEdgeTime() << " " << p->leadingEdgeTime() << " " << p->associatedTrack() << " " << p->matchFlag() << " " << p->position() << endm;
01796       }
01797     }
01798   }
01799   //-- end check
01800   
01801   LOG_INFO  << "F: before/after" << FinalMatchedCellsVec.size() << "/" <<nValidSinglePrimHitCells << endm;
01802  // end of Sect.F
01803 
01804   LOG_INFO  << "#(cell tracks): " << allCellsHitVec.size()
01805        << " #(hit cells): " << FinalMatchedCellsVec.size()
01806        << "\n#(single hits): " << nSingleHitCells 
01807        << " #(single valid hits): " << nValidSingleHitCells
01808        << " #(single prim valid hits): " << nValidSinglePrimHitCells
01809        << endm;
01810 
01811 
01812 
01813   if (doPrintMemoryInfo) {
01814         StMemoryInfo::instance()->snapshot();
01815         StMemoryInfo::instance()->print();
01816   }
01817   if (doPrintCpuInfo) {
01818     timer.stop();
01819     LOG_INFO  << "CPU time for StTofrMatchMaker::Make(): "
01820          << timer.elapsedTime() << " sec\n" << endm;
01821   }
01822 
01823   LOG_INFO  << "StTofrMatchMaker -- bye-bye" << endm;
01824 
01825 
01826 
01827   return kStOK;
01828 }
01829 
01830 //---------------------------------------------------------------------------
01831 Int_t StTofrMatchMaker::processEventYear8(){
01832   // leave as empty now
01833 
01834   if(Debug()) LOG_INFO << " processing event in run 8 " << endm;
01835   if(mHisto) mEventCounterHisto->Fill(0);
01836   // event selection ...
01837   mEvent = (StEvent *) GetInputDS("StEvent");
01838 //  if (!validEvent(mEvent)){
01839   if(!mEvent || !(mEvent->tofCollection()) || !(mEvent->tofCollection()->rawdataPresent()) ) {
01840     gMessMgr->Info("StTofrMatchMaker -- nothing to do ... bye-bye","OS");
01841     return kStOK;
01842   }
01843   if(mHisto) mEventCounterHisto->Fill(1);
01844 
01845   // timing & memory info -only when requested-
01846   StTimer timer;
01847   if (doPrintCpuInfo) timer.start();
01848   if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
01849 
01850   //.........................................................................
01851   // check for tofCollection and fill local copy with ADC and TDC data
01852   StTofCollection *theTof = mEvent->tofCollection();
01853 
01854   //.........................................................................
01855   // push raw data into StTofData for QA:
01856   // Be careful, the vpd trayId in TofData is 121 and 122 (**)
01857   //  and
01858   //.........................................................................
01859   // A. build vector of candidate cells with valid ADC signals 
01860   //  idVector validCellIdVec;
01861   tofCellHitVector daqCellsHitVec;
01862   //  daqCellsHitVec.clear();
01863   idVector validModuleVec;
01864 
01865   mSortTofRawData = new StSortTofRawData(theTof, mDaqMap);
01866 
01867   // multi-tray system
01868   IntVec validtray = mDaqMap->ValidTrays();
01869   for(size_t i=0;i<validtray.size();i++) {
01870     int trayId = validtray[i];
01871     IntVec validchannel = mSortTofRawData->GetValidChannel(trayId);
01872     if(Debug()) LOG_INFO  << " Number of fired hits on tray " << trayId << " = " << validchannel.size() << endm;
01873 
01874     for(size_t iv=0;iv<validchannel.size();iv++) {
01875       IntVec leTdc = mSortTofRawData->GetLeadingTdc(trayId, validchannel[iv], kTRUE);
01876       IntVec teTdc = mSortTofRawData->GetTrailingTdc(trayId, validchannel[iv], kTRUE);
01877 
01878       if(!leTdc.size() || !teTdc.size()) continue;
01879 
01880       int chan = validchannel[iv];
01881       IntVec map = mDaqMap->TDIGChan2Cell(chan);
01882       int moduleId = map[0];
01883       int cellId = map[1];
01884       
01885       StructCellHit aDaqCellHit;
01886       aDaqCellHit.channel = chan;
01887       aDaqCellHit.tray = trayId;
01888       aDaqCellHit.module = moduleId;
01889       aDaqCellHit.cell = cellId;
01890       daqCellsHitVec.push_back(aDaqCellHit);
01891       
01892       // additional valid number configuration
01893       int id = trayId*100+moduleId;
01894       bool ifind = kFALSE;
01895       for(size_t im=0;im<validModuleVec.size();im++) {
01896         if(id==validModuleVec[im]) {
01897           ifind = kTRUE;
01898           break;
01899         }
01900       }
01901       if(!ifind) validModuleVec.push_back(id);
01902 
01903       if(mHisto) {
01904   //      mDaqOccupancyValid->Fill(chan);
01905         mDaqOccupancyValid->Fill((moduleId-1)*mNCell+(cellId-1));
01906         mDaqOccupancyValidAll->Fill((trayId-76)*mNTOF+(moduleId-1)*mNCell+(cellId-1));
01907       }      
01908       //
01909       // store data from trays and vpds into StTofData
01910       //
01911       int dataIndex = (trayId-1)*mNTOF + (moduleId-1)*mNCell + (cellId-1);
01912       StTofData *aData = new StTofData(dataIndex,0,0,0,0,leTdc[0],teTdc[0]);
01913       theTof->addData(aData);
01914 
01915       if(Debug()) {
01916         for(size_t iv1=0;iv1<leTdc.size();iv1++) {
01917           LOG_INFO  << " leading Tdc = " << leTdc[iv1]<<endm;
01918         }
01919         for(size_t iv2=0;iv2<teTdc.size();iv2++) {
01920           LOG_INFO  << " trailing Tdc = " << teTdc[iv2] << endm;
01921         }
01922       } // end debug
01923     } // end channel
01924 
01925   } // end tray
01926 
01927   // vpd -> StTofData
01928   for(int ivpd=0;ivpd<2;ivpd++) { // west and east sides
01929     int trayId = (ivpd==0) ? mWestVpdTrayId : mEastVpdTrayId;
01930     IntVec validtube = mSortTofRawData->GetValidChannel(trayId);
01931     if(Debug()) LOG_INFO  << " Number of fired hits on tray(vpd) " << trayId << " = " << validtube.size() << endm;
01932 
01933     if(!validtube.size()) continue;
01934     for(int i=0;i<mNVPD;i++) {
01935       int tubeId = i+1;
01936       int lechan = (ivpd==0) ? mDaqMap->WestPMT2TDIGLeChan(tubeId) : mDaqMap->EastPMT2TDIGLeChan(tubeId);
01937       IntVec leTdc = mSortTofRawData->GetLeadingTdc(trayId, lechan, kTRUE);
01938       IntVec teTdc = mSortTofRawData->GetTrailingTdc(trayId, lechan, kTRUE);  // channel number should be le, sorted in StSortTofRawData
01939 
01940       if(leTdc.size()&&mHisto) mDaqOccupancyVpd->Fill(ivpd*mNVPD+i);
01941 
01942       if(leTdc.size() && teTdc.size()) {
01943         int dataIndex = (ivpd+120)*mNTOF + (tubeId-1);
01944         StTofData *aData = new StTofData(dataIndex,0,0,0,0,leTdc[0],teTdc[0]);
01945         theTof->addData(aData);
01946 
01947         if(mHisto) mDaqOccupancyValidVpd->Fill(ivpd*mNVPD+i);
01948       }
01949     }
01950   }
01951   //
01952 
01953   // end of Sect.A
01954   if(Debug()) {
01955     LOG_INFO  << "    total # of cells = " << daqCellsHitVec.size() << endm;
01956     for(size_t iv = 0;iv<validModuleVec.size();iv++) {
01957       LOG_INFO  << " module # " << validModuleVec[iv] << " Valid! " << endm;
01958     }
01959   }
01960   if(mHisto) {
01961     mCellsMultInEvent->Fill(daqCellsHitVec.size());
01962     if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
01963   }
01964 //  if(!daqCellsHitVec.size()) return kStOK;
01965 
01966   //.........................................................................
01967   // B. loop over global tracks and determine all cell-track matches
01968   //
01969   tofCellHitVector allCellsHitVec;
01970   //  allCellsHitVec.clear();
01971   StructCellHit cellHit;
01972 
01973   StSPtrVecTrackNode& nodes = mEvent->trackNodes();
01974   Int_t nAllTracks=0;
01975   for (unsigned int iNode=0; iNode<nodes.size(); iNode++){
01976     tofCellHitVector cellHitVec;
01977     //    cellHitVec.clear();
01978     StGlobalTrack *theTrack = dynamic_cast<StGlobalTrack*>(nodes[iNode]->track(global));
01979     if(!theTrack) continue;
01980 
01981     StThreeVectorF mom = theTrack->geometry()->momentum();
01982     float pt = mom.perp();
01983     float eta = mom.pseudoRapidity();
01984     float phi = mom.phi();
01985     if (phi<0.) phi += 2.*3.14159;
01986 
01987     float dEdx = -999.;
01988     int ndEdxpts = 0;
01989     float nSigmaPion = -999.;
01990     static StTpcDedxPidAlgorithm PidAlgorithm;
01991     static StPionPlus* Pion = StPionPlus::instance();
01992 
01993     const StParticleDefinition* pd = theTrack->pidTraits(PidAlgorithm);
01994     if (pd) {
01995       nSigmaPion = PidAlgorithm.numberOfSigma(Pion);
01996     }
01997     if ( PidAlgorithm.traits() ) {
01998       dEdx = PidAlgorithm.traits()->mean();
01999       ndEdxpts = PidAlgorithm.traits()->numberOfPoints();
02000     }
02001     int nfitpts = theTrack->fitTraits().numberOfFitPoints(kTpcId);
02002 //    cout << " dEdx = " << dEdx << endl;
02003 //    cout << " nSigmaPi = " << nSigmaPion << endl;
02004 //    if (pt<0.2) continue;
02005 
02006     // make sure we have a track, a miniDST might have removed it...
02007     if (validTrackRun8(theTrack)){
02008       if(mHisto) {
02009         mTrackPtEta->Fill(pt, eta);
02010         mTrackPtPhi->Fill(pt, phi);
02011         mTrackNFitPts->Fill(nfitpts);
02012         if(dEdx>0.) mTrackdEdxvsp->Fill(mom.mag(), dEdx*1.e6);
02013         if(fabs(nSigmaPion)<5.) mNSigmaPivsPt->Fill(pt, nSigmaPion+5.*theTrack->geometry()->charge());
02014       }
02015 
02016       if(mSaveTree) {
02017         trackTree.pt = pt;
02018         trackTree.eta = eta;
02019         trackTree.phi = phi;
02020         trackTree.nfitpts = nfitpts;
02021         trackTree.dEdx = dEdx*1.e6;
02022         trackTree.ndEdxpts = ndEdxpts;
02023         trackTree.charge = theTrack->geometry()->charge();      
02024         trackTree.projTrayId = 0;
02025         trackTree.projCellChan = -1;
02026         trackTree.projY = -999.;
02027         trackTree.projZ = -999.;
02028       }
02029 
02030       nAllTracks++;
02031       StPhysicalHelixD theHelix = trackGeometry(theTrack)->helix();
02032 
02033 //      IntVec projTrayVec;
02034 //      if(!mTofrGeom->projTrayVector(theHelix, projTrayVec)) continue;
02035 
02036       IntVec idVec;
02037       DoubleVec pathVec;
02038       PointVec  crossVec;
02039 
02040 //       idVec.clear();
02041 //       pathVec.clear();
02042 //       crossVec.clear();
02043 
02044       Int_t ncells = 0;
02045       if(mTofrGeom->HelixCrossCellIds(theHelix,idVec,pathVec,crossVec) ) {
02046 //      if(mTofrGeom->HelixCrossCellIds(theHelix, validModuleVec, projTrayVec, idVec, pathVec, crossVec)) {
02047         Int_t cells = idVec.size();
02048         for (Int_t i=0; i<cells; i++) {
02049             Int_t icell,imodule,itray;
02050             Double_t local[3],global[3];
02051             for(Int_t i2=0;i2<3;i2++){
02052                  local[i2]=0;
02053             }
02054             global[0]=crossVec[i].x();
02055             global[1]=crossVec[i].y();
02056             global[2]=crossVec[i].z();
02057             mTofrGeom->DecodeCellId(idVec[i], icell, imodule, itray);
02058             LOG_INFO << " decode " << idVec[i] << "  to tray#" << itray << " module#" << imodule << " cell#" << icell << endm;
02059             StTofrGeomSensor* sensor = 
02060                   mTofrGeom->GetGeomSensor(imodule,itray);
02061             if(!sensor) {
02062               gMessMgr->Warning("","OS") << " No sensitive module in the projection??? -- Something weird!!! " << endm;
02063               continue;
02064             }
02065             sensor->Master2Local(&global[0],&local[0]);
02066             icell = sensor->FindCellIndex(local);
02067             //      StThreeVectorD glo=sensor->GetCenterPosition();
02068             StThreeVectorD glo(global[0], global[1], global[2]);
02069             StThreeVectorD hitPos(local[0], local[1], local[2]);
02070             delete sensor;
02071             if (local[2]<=2.7&&local[2]>=-3.4) {
02072               Int_t Iarray = mDaqMap->Cell2TDIGChan(imodule, icell);
02073               if(Iarray>=mDAQOVERFLOW||Iarray<0) continue;
02074               ncells++;
02075               cellHit.channel = Iarray;
02076               cellHit.tray = itray;
02077               cellHit.module = imodule;
02078               cellHit.cell = icell;
02079               cellHit.trackIdVec.push_back(iNode);
02080               cellHit.hitPosition = glo;        // global position
02081               cellHit.zhit = (Float_t)hitPos.z();
02082               cellHit.yhit = (Float_t)hitPos.y();
02083               cellHitVec.push_back(cellHit);
02084               allCellsHitVec.push_back(cellHit);
02085               if(mHisto) {
02086 //              mDaqOccupancyProj->Fill(Iarray);
02087                 if(itray>=76&&itray<=80) {
02088                   mDaqOccupancyProj->Fill((imodule-1)*mNCell+(icell-1));
02089                   mDaqOccupancyProjAll->Fill((itray-76)*mNTOF+(imodule-1)*mNCell+(icell-1));
02090                   mHitsPosition->Fill(hitPos.y(), hitPos.z());
02091                 }
02092               }
02093               
02094               if(mSaveTree) {
02095                 trackTree.projTrayId = itray;
02096                 trackTree.projCellChan = (imodule-1)*mNCell+(icell-1);
02097                 trackTree.projY = local[1];
02098                 trackTree.projZ = local[2];
02099               }
02100 
02101               if(Debug()) {
02102                 LOG_INFO  <<"B: nodeid=" << iNode << "  projected in " << " tray="<< itray << " module="<<imodule<<" cell="<<icell<<endm;
02103                 LOG_INFO  <<"   hit position " << hitPos << endm;
02104               }
02105             }
02106         } // for (Int_t i=0...)
02107       } // endif(helixcross...)
02108       if(ncells>0&&mHisto) mHitsMultPerTrack->Fill(ncells);
02109 
02110       if(mHisto && mSaveTree) mTrackTree->Fill();
02111 
02112     } // if(ValidTrack).. 
02113   } // loop over nodes
02114   if(Debug())
02115     LOG_INFO  << "B:  matched/available/total #tracknodes: " <<allCellsHitVec.size() << "/" <<nAllTracks << "/" << nodes.size() << endm;
02116   if(mHisto) {
02117     mHitsMultInEvent->Fill(allCellsHitVec.size());
02118     if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
02119   }
02120   // end of Sect.B
02121 
02122   //.........................................................................
02123   // C. Match find Neighbours -- identify crosstalk
02124   //
02125   tofCellHitVector matchHitCellsVec;
02126   //  matchHitCellsVec.clear();
02127 
02128   tofCellHitVectorIter daqIter = daqCellsHitVec.begin();
02129   for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
02130     tofCellHitVectorIter proIter = allCellsHitVec.begin();
02131     for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
02132 
02133       int daqIndex = (daqIter->module-1)*6 + (daqIter->cell-1);
02134       int proIndex = (proIter->module-1)*6 + (proIter->cell-1);
02135       int hisIndex = daqIter->tray - 76;
02136       int daqAllIndex = (daqIter->tray - 76)*192 + daqIndex;
02137       int proAllIndex = (proIter->tray - 76)*192 + proIndex;
02138       if (mHisto) mHitCorrAll->Fill(proAllIndex,daqAllIndex);
02139       if(daqIter->tray==proIter->tray) {
02140         if (mHisto) {
02141           if(hisIndex>=0&&hisIndex<5) {
02142             mHitCorr[hisIndex]->Fill(proIndex,daqIndex);
02143             mHitCorrModule[hisIndex]->Fill(proIter->module-1,daqIter->module-1);
02144           } else {
02145             cout << " weird tray # " << daqIter->tray << endl;
02146           }
02147         }
02148       }
02149       if( (daqIter->tray==proIter->tray)&& 
02150           (daqIter->module==proIter->module) &&
02151           ( ( (proIter->cell==6)&&((proIter->cell==daqIter->cell) ||
02152                                    (proIter->cell==daqIter->cell+1)) )
02153             || ( (proIter->cell==1)&&((proIter->cell==daqIter->cell) ||
02154                                       (proIter->cell==daqIter->cell-1)) )
02155             || ( (proIter->cell>=2&&proIter->cell<=6) &&
02156                  ( (proIter->cell==daqIter->cell) ||
02157                    (proIter->cell==daqIter->cell-1) ||
02158                    (proIter->cell==daqIter->cell+1) ) ) ) ) {
02159         cellHit.channel = daqIter->channel;
02160         cellHit.tray = daqIter->tray;
02161         cellHit.module = daqIter->module;
02162         cellHit.cell = daqIter->cell;
02163         cellHit.hitPosition = proIter->hitPosition;
02164         cellHit.trackIdVec = proIter->trackIdVec;
02165         cellHit.zhit = proIter->zhit;
02166         cellHit.yhit = proIter->yhit;
02167         matchHitCellsVec.push_back(cellHit);
02168       }
02169     }
02170   } //end {sec. C}
02171   if(Debug()) {
02172     LOG_INFO  << "C: before/after: " << allCellsHitVec.size() << "/" << matchHitCellsVec.size() << endm;
02173   }
02174   if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
02175 
02176   //.........................................................................
02177   // D. sort hit vectors  and deal with (discard) cells matched by multiple tracks
02178   //
02179   Int_t nSingleHitCells(0);
02180   Int_t nMultiHitsCells(0);
02181 
02182   tofCellHitVector singleHitCellsVec;
02183   tofCellHitVector multiHitsCellsVec;
02184 //   singleHitCellsVec.clear();
02185 //   multiHitsCellsVec.clear();
02186 
02187   tofCellHitVector tempVec = matchHitCellsVec;
02188   tofCellHitVector erasedVec = tempVec;
02189   while (tempVec.size() != 0) {
02190     Int_t nTracks = 0;
02191     idVector trackIdVec;
02192 
02193     tofCellHitVectorIter tempIter=tempVec.begin();
02194     tofCellHitVectorIter erasedIter=erasedVec.begin();
02195     while(erasedIter!= erasedVec.end()) {
02196       if(tempIter->tray == erasedIter->tray &&
02197          tempIter->module == erasedIter->module &&
02198          tempIter->cell == erasedIter->cell) {
02199         nTracks++;
02200         trackIdVec.push_back(erasedIter->trackIdVec.back());  // merge
02201         erasedVec.erase(erasedIter);
02202         erasedIter--;
02203       }
02204       erasedIter++;
02205     }
02206 
02207     cellHit.channel = tempIter->channel;
02208     cellHit.cell = tempIter->cell;
02209     cellHit.module = tempIter->module;
02210     cellHit.tray = tempIter->tray;
02211     cellHit.hitPosition = tempIter->hitPosition;
02212     cellHit.trackIdVec = trackIdVec;
02213     cellHit.zhit = tempIter->zhit;
02214     cellHit.yhit = tempIter->yhit;
02215 
02216     Float_t ycenter = (tempIter->cell-1-2.5)*mWidthPad;
02217     Float_t dy = tempIter->yhit - ycenter;
02218     Float_t dz = tempIter->zhit;
02219 
02220     if(mHisto) {
02221       mTracksPerCellMatch1->Fill(trackIdVec.size());
02222 //      mDaqOccupancyMatch1->Fill(tempIter->channel);
02223       mDaqOccupancyMatch1->Fill((tempIter->module-1)*mNCell+(tempIter->cell-1));
02224       mDeltaHitMatch1->Fill(dy, dz);
02225     }
02226 
02227     if (nTracks==1){
02228       nSingleHitCells++;      
02229       singleHitCellsVec.push_back(cellHit);
02230     } else if (nTracks>1){
02231       nMultiHitsCells++;
02232       multiHitsCellsVec.push_back(cellHit);
02233       // for multiple hit cells either discard (yes) or
02234       // find the most likely candidate.
02235     } else {
02236       LOG_INFO  << "D: no tracks extrapolate to matched cell ... should not happen!" << endm;
02237     }
02238     
02239     if (Debug()) {
02240       LOG_INFO  << "D: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
02241       idVectorIter ij=trackIdVec.begin();
02242       while (ij != trackIdVec.end()) { LOG_INFO  << " " << *ij; ij++; }
02243       LOG_INFO  <<endm;
02244     }
02245     
02246     tempVec = erasedVec;
02247   }
02248   if(Debug())
02249     LOG_INFO  << "D: before/after: " << matchHitCellsVec.size() << "/" << singleHitCellsVec.size() << endm;
02250   //end of Sect.C
02251   if(mHisto) {
02252     mCellsPerEventMatch1->Fill(singleHitCellsVec.size()+multiHitsCellsVec.size());
02253     if(singleHitCellsVec.size()) mEventCounterHisto->Fill(9);
02254   } 
02255 
02256   //.........................................................................
02257   // E. sort and deal singleHitCellsVector for multiple cells associated to single tracks
02258   //
02259   tofCellHitVector FinalMatchedCellsVec;
02260   //  FinalMatchedCellsVec.clear();
02261   tempVec = singleHitCellsVec;
02262   if(mHisto) {
02263     mCellsPerEventMatch2->Fill(tempVec.size());
02264     for(unsigned int ii=0;ii<tempVec.size();ii++) {
02265       mTracksPerCellMatch2->Fill(tempVec[ii].trackIdVec.size());
02266  //     mDaqOccupancyMatch2->Fill(tempVec[ii].channel);
02267       mDaqOccupancyMatch2->Fill((tempVec[ii].module-1)*mNCell+(tempVec[ii].cell-1));
02268       Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
02269       Float_t dy = tempVec[ii].yhit-ycenter;
02270       Float_t dz = tempVec[ii].zhit;
02271       mDeltaHitMatch2->Fill(dy, dz);
02272     }
02273   }
02274 
02275   erasedVec = tempVec;
02276   while (tempVec.size() != 0) {
02277     StructCellHit cellHit;
02278     Int_t nCells = 0;
02279     idVector vTrackId;
02280     vector<StThreeVectorD> vPosition;
02281     vector<Int_t> vchannel, vtray, vmodule, vcell;
02282     vector<Float_t> vzhit, vyhit;
02283 
02284     tofCellHitVectorIter tempIter=tempVec.begin();
02285     tofCellHitVectorIter erasedIter=erasedVec.begin();
02286     while(erasedIter!= erasedVec.end()) {
02287       if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
02288         nCells++;
02289         vchannel.push_back(erasedIter->channel);
02290         vtray.push_back(erasedIter->tray);
02291         vmodule.push_back(erasedIter->module);
02292         vcell.push_back(erasedIter->cell);
02293         vPosition.push_back(erasedIter->hitPosition);
02294         vTrackId.push_back(erasedIter->trackIdVec.back());
02295         vzhit.push_back(erasedIter->zhit);
02296         vyhit.push_back(erasedIter->yhit);
02297 
02298         erasedVec.erase(erasedIter);
02299         erasedIter--;
02300       }
02301       erasedIter++;
02302     }
02303 
02304     if (nCells==1){
02305       // for singly hit cell, copy data in singleHitCellsVec
02306       cellHit.channel = vchannel[0];
02307       cellHit.tray = vtray[0];
02308       cellHit.module = vmodule[0];
02309       cellHit.cell = vcell[0];
02310       cellHit.trackIdVec.push_back(vTrackId[0]);
02311       cellHit.hitPosition = vPosition[0];
02312       cellHit.matchFlag = 0; 
02313       cellHit.zhit = vzhit[0];
02314       cellHit.yhit = vyhit[0];
02315 
02316       FinalMatchedCellsVec.push_back(cellHit);
02317 
02318       // debugging output
02319       if (Debug()) {
02320         LOG_INFO  << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
02321         idVectorIter ij=vTrackId.begin();
02322         while (ij != vTrackId.end()) { LOG_INFO  << " " << *ij; ij++; }
02323         LOG_INFO  <<endm;
02324       }
02325     }
02326     else if (nCells>1){   // for multiple hit cells  find the most likely candidate.
02327       Int_t thiscandidate(-99);
02328       Int_t thisMatchFlag(0);
02329 
02330       // sort on hitposition
02331       Float_t ss(99.);
02332       vector<Int_t> ssCandidates;
02333       thisMatchFlag = 2;
02334       if (Debug()) LOG_INFO  << " ss " << endm;
02335       for (Int_t i=0;i<nCells;i++){
02336         Float_t yy = vyhit[i];
02337         Float_t ycell = (vcell[i]-1-2.5)*mWidthPad;
02338         Float_t ll = fabs(yy-ycell);
02339         if(ll<ss) {
02340           ss = ll;
02341           ssCandidates.clear();
02342           ssCandidates.push_back(i);
02343         }else if  (ll==ss)
02344           ssCandidates.push_back(i);      
02345       }
02346       if (ssCandidates.size()==1){
02347         thiscandidate = ssCandidates[0];
02348         Int_t daqId = vchannel[thiscandidate];
02349         if (Debug()) LOG_INFO  << "candidate =" << daqId << endm;
02350       }
02351       
02352 
02353       if (thiscandidate>=0) {
02354         cellHit.channel = vchannel[thiscandidate];
02355         cellHit.tray = vtray[thiscandidate];
02356         cellHit.module = vmodule[thiscandidate];
02357         cellHit.cell = vcell[thiscandidate];
02358         cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
02359         cellHit.hitPosition = vPosition[thiscandidate];
02360         cellHit.matchFlag = thisMatchFlag;
02361         cellHit.zhit = vzhit[thiscandidate];
02362         cellHit.yhit = vyhit[thiscandidate];
02363 
02364         FinalMatchedCellsVec.push_back(cellHit);
02365         
02366         // debugging output
02367         if (Debug()) {
02368         LOG_INFO  << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:" << vTrackId[thiscandidate] << endm;
02369         }
02370       }
02371 
02372     } else {
02373       LOG_INFO  << "E: no cells belong to this track ... should not happen!" << endm;
02374     }
02375 
02376     tempVec = erasedVec;
02377   }
02378 
02379   LOG_INFO  << "E: before/after: " << singleHitCellsVec.size() << "/" << FinalMatchedCellsVec.size() << endm;
02380   // end of Sect.E
02381 
02382   //.........................................................................
02383   // F. perform further selection and
02384   //    fill valid track histograms, ntuples and CellCollection
02385   //
02386   tempVec.clear();
02387   tempVec = FinalMatchedCellsVec;
02388   if(mHisto) {
02389     if(FinalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
02390     mCellsPerEventMatch3->Fill(tempVec.size());
02391     for(unsigned int ii=0;ii<tempVec.size();ii++) {
02392       mTracksPerCellMatch3->Fill(tempVec[ii].trackIdVec.size());
02393 //      mDaqOccupancyMatch3->Fill(tempVec[ii].channel);
02394       mDaqOccupancyMatch3->Fill((tempVec[ii].module-1)*mNCell+(tempVec[ii].cell-1));
02395       Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
02396       Float_t dy = tempVec[ii].yhit - ycenter;
02397       Float_t dz = tempVec[ii].zhit;
02398       mDeltaHitMatch3->Fill(dy, dz);
02399     }
02400   }
02401 
02402   StTofCellCollection *mCellCollection =  new StTofCellCollection;
02403   Int_t nValidSingleHitCells(0), nValidSinglePrimHitCells(0);
02404 
02405   for (size_t ii=0; ii < FinalMatchedCellsVec.size(); ii++){
02406     Int_t daqId = FinalMatchedCellsVec[ii].channel;
02407     Int_t jj = daqId;
02408     Int_t tray = FinalMatchedCellsVec[ii].tray;
02409     Int_t module = FinalMatchedCellsVec[ii].module;
02410     Int_t cell = FinalMatchedCellsVec[ii].cell;
02411 
02412 //    Float_t ycenter = (cell-1-2.5)*mWidthPad;
02413 //    Float_t dy = FinalMatchedCellsVec[ii].yhit - ycenter;
02414     if (FinalMatchedCellsVec[ii].trackIdVec.size()!=1)
02415       LOG_INFO  << "F: WHAT!?!  mult.matched cell in single cell list " << daqId << endm;
02416 
02417 
02418     // Read in Leading and Trailing edge TDC, apply on INL correction
02419     //
02420     int tmptdc = (mSortTofRawData->GetLeadingTdc(tray,jj,kTRUE))[0];
02421     int bin = (int)tmptdc&0x3ff;
02422     double tmptdc_f = tmptdc + mTofINLCorr->getTrayINLCorr(tray, jj, bin);
02423     double letime = tmptdc_f*VHRBIN2PS / 1000.;
02424 
02425     tmptdc=(mSortTofRawData->GetTrailingTdc(tray,jj,kTRUE))[0];
02426     bin = (int)tmptdc&0x3ff;
02427     tmptdc_f = tmptdc + mTofINLCorr->getTrayINLCorr(tray, jj, bin);
02428     double tetime = tmptdc_f*VHRBIN2PS / 1000.;
02429 
02430     // get track-id from cell hit vector
02431     unsigned int trackNode = FinalMatchedCellsVec[ii].trackIdVec[0];
02432 //    StTrack *theTrack = nodes[trackNode]->track(primary);
02433     StTrack *globalTrack = nodes[trackNode]->track(global);
02434 
02435     // 2. continue only if the (primary) track exists
02436 //    if (validTofTrack(theTrack) && fabs(dy)<1.9 ){
02437       nValidSinglePrimHitCells++;
02438 
02439       //--- store number of hits per track
02440  //     Int_t nHitsPerTrack = theTrack->topologyMap().numberOfHits(kTpcId);
02441           
02442       // select the apropriate track geometry
02443  //     StTrackGeometry *theTrackGeometry = trackGeometry(globalTrack);
02444 
02445       //--- calculate local hit position on cell based on average hitposition
02446       //      Float_t localHitPos = mTofGeom->cellHitPosition(&allMatchedCellsVec[ii].hitPosition);
02447       
02448       // Fill TOF Cell Collection
02449       StTofCell *tofCell = new StTofCell(tray, module, cell, daqId, globalTrack, FinalMatchedCellsVec[ii].zhit, FinalMatchedCellsVec[ii].matchFlag, FinalMatchedCellsVec[ii].hitPosition);
02450       tofCell->setLeadingEdgeTime(letime);
02451       tofCell->setTrailingEdgeTime(tetime);
02452       mCellCollection->push_back(tofCell);
02453       
02454       // dump debug data
02455       if (Debug()){
02456         LOG_INFO  << "F: itray=" << tray << " imodule=" << module << " icell=" << cell << "\tnodeid:";
02457         idVectorIter ij=FinalMatchedCellsVec[ii].trackIdVec.begin();
02458         while (ij != FinalMatchedCellsVec[ii].trackIdVec.end()) { LOG_INFO  << " " << *ij; ij++; }
02459         LOG_INFO  << endm;
02460       }
02461 
02462 //    } // track exists 
02463   } // end final matched cells
02464 
02465 
02466   // put INL corrected vpd information in StTofCell as well
02467   StSPtrVecTofData &tofData = theTof->tofData();
02468   for(size_t ii = 0; ii<tofData.size(); ii++) {
02469     StTofData *aData = dynamic_cast<StTofData *>(tofData[ii]);
02470     if(!aData) continue;
02471 
02472     int dataIndex = aData->dataIndex();
02473     int trayId = dataIndex / mNTOF;
02474     if(trayId<120) continue;  // only vpd selected
02475     int ewId = trayId - 120 + 1;  // 1: west,  2: east
02476     
02477     int tubeId = dataIndex % mNTOF + 1;
02478     int lechan = (ewId==1) ? mDaqMap->WestPMT2TDIGLeChan(tubeId) : mDaqMap->EastPMT2TDIGLeChan(tubeId);
02479     int techan = (ewId==1) ? mDaqMap->WestPMT2TDIGTeChan(tubeId) : mDaqMap->EastPMT2TDIGTeChan(tubeId);
02480 
02481     int tmptdc = aData->leadingTdc();
02482     int bin = (int)tmptdc&0x3ff;
02483     double tmptdc_f = tmptdc + mTofINLCorr->getVpdINLCorr(ewId, lechan, bin);
02484     double letime = tmptdc_f*VHRBIN2PS / 1000.;
02485 
02486     tmptdc = aData->trailingTdc();
02487     bin = (int)tmptdc&0x3ff;
02488     tmptdc_f = tmptdc + mTofINLCorr->getVpdINLCorr(ewId, techan, bin);
02489     double tetime = tmptdc_f*VHRBIN2PS / 1000.;
02490 
02491     StThreeVectorF zero(0.,0.,0.);
02492     StTofCell *tofCell = new StTofCell(120+ewId, 0, tubeId, lechan, 0, 0, 0, zero);
02493     tofCell->setLeadingEdgeTime(letime);
02494     tofCell->setTrailingEdgeTime(tetime);
02495     mCellCollection->push_back(tofCell);
02496 
02497     if (Debug()){
02498       LOG_INFO  << "F: itray=" << trayId << " imodule=" << 0 << " itube=" << tubeId << " letime=" << letime << " tetime=" << tetime << endm;
02499     }
02500   }
02501 
02502   
02503   storeMatchData(mCellCollection,theTof);
02504   delete mCellCollection;
02505 
02506   delete mSortTofRawData;
02507   mSortTofRawData = 0;
02508   
02509   //check StEvent collections --
02510   if (theTof->dataPresent())
02511     LOG_INFO  << " TofCollection: raw data container present" << endm;
02512   if (theTof->cellsPresent()){
02513     LOG_INFO  << " TofCollection: cell container present."<<endm;
02514     if (Debug()){
02515       StSPtrVecTofCell& tmpCellTofVec = theTof->tofCells();
02516       LOG_INFO  << " # of matched cells " << tmpCellTofVec.size() << endm;
02517       for (size_t i = 0; i < tmpCellTofVec.size(); i++) {
02518         StTofCell* p = tmpCellTofVec[i];
02519         LOG_INFO  << p->trayIndex() << " " << p->moduleIndex() << " " << p->cellIndex() << " " << p->trailingEdgeTime() << " " << p->leadingEdgeTime() << " " << p->associatedTrack() << " " << p->matchFlag() << " " << p->position() << endm;
02520       }
02521     }
02522   }
02523   //-- end check
02524   
02525   LOG_INFO  << "F: before/after" << FinalMatchedCellsVec.size() << "/" <<nValidSinglePrimHitCells << endm;
02526  // end of Sect.F
02527 
02528   LOG_INFO  << "#(cell tracks): " << allCellsHitVec.size()
02529        << " #(hit cells): " << FinalMatchedCellsVec.size()
02530        << "\n#(single hits): " << nSingleHitCells 
02531        << " #(single valid hits): " << nValidSingleHitCells
02532        << " #(single prim valid hits): " << nValidSinglePrimHitCells
02533        << endm;
02534 
02535 
02536 
02537   if (doPrintMemoryInfo) {
02538         StMemoryInfo::instance()->snapshot();
02539         StMemoryInfo::instance()->print();
02540   }
02541   if (doPrintCpuInfo) {
02542     timer.stop();
02543     LOG_INFO  << "CPU time for StTofrMatchMaker::Make(): "
02544          << timer.elapsedTime() << " sec\n" << endm;
02545   }
02546 
02547   LOG_INFO  << "StTofrMatchMaker -- bye-bye" << endm;
02548 
02549 
02550 
02551   return kStOK;
02552 }
02553 
02554 //---------------------------------------------------------------------------
02555 // store local slat collection in StEvent's tofCollection
02556 Int_t StTofrMatchMaker::storeMatchData(StTofCellCollection *cellCollection,
02557                                          StTofCollection* tofCollection){
02558   if(!tofCollection){
02559     LOG_INFO  << "Error: no TofCollection -- returning" << endm;
02560     return kStErr;
02561   }
02562 
02563   for (size_t j=0;j<cellCollection->size();j++){
02564     tofCollection->addCell(cellCollection->getCell(j)); 
02565     if (Debug())
02566       LOG_INFO  << "storing " << j << "  " << "  tray:"
02567            << cellCollection->getCell(j)->trayIndex() << "  module:"
02568            << cellCollection->getCell(j)->moduleIndex() << "  cell:"
02569            << cellCollection->getCell(j)->cellIndex() << endm;
02570   }
02571   return kStOK;
02572 }
02573 
02574 
02575 //---------------------------------------------------------------------------
02576 // create a local copy of the raw tofp data tofData in StEvent's tofCollection
02577 Int_t StTofrMatchMaker::getTofData(StTofCollection* tofCollection){
02578   if (!tofCollection) return kStERR;
02579   LOG_INFO  << " Read in tof data ... " << endm;
02580   StSPtrVecTofData &tofData = tofCollection->tofData();
02581 
02582   // perform consistency check
02583   bool dataOK(true);
02584   LOG_INFO  << "TOF raw data consistency test ...";
02585 
02586   Int_t tofsize = tofData.size();
02587   Int_t nTOFr = 0;
02588   if(mYear4) {
02589     if(tofsize<184) {
02590       gMessMgr->Warning("The size of tofData is NOT 184!","OS");
02591       nTOFr = 72; // DAQ not updated
02592     } else {
02593       nTOFr = mNTOFR;
02594     }
02595   } else if(mYear3) {
02596     nTOFr = 72;
02597   } else if(mYear2) {
02598     nTOFr = 0;
02599   }
02600 
02601   for (Int_t i=0;i<nTOFr;i++){
02602     Int_t iAdc = mDaqMap->DaqChan2ADCChan(i);
02603     mTofrAdc[i] = tofData[iAdc]->adc();
02604     Int_t iTdc = mDaqMap->DaqChan2TDCChan(i);
02605     mTofrTdc[i] = tofData[iTdc]->tdc();
02606   }
02607 
02608   for (Int_t i=0;i<6;i++){
02609     mPvpdAdc[i] = tofData[42+i]->adc(); 
02610     mPvpdTdc[i] = tofData[42+i]->tdc(); 
02611     if (mYear3||mYear4)
02612       mPvpdAdcLoRes[i] = tofData[54+i]->adc();
02613   }
02614   
02615   if (!dataOK) return kStWarn;
02616 
02617   return kStOK;
02618 }
02619 
02620 
02621 //---------------------------------------------------------------------------
02622 // Book histograms and create ordered collection for easy manipulation
02623 void StTofrMatchMaker::bookHistograms(void){
02624 
02625   mEventCounterHisto = new TH1D("eventCounter","eventCounter",20,0,20);
02626 
02627   mADCTDCCorelation = new TH2D("ADCTDCCorelation","ADCTDCCorelation",200,0,200,200,0,200);
02628   mCellsMultInEvent = new TH1D("cellsPerEvent","cellsPerEvent",100,0,100);
02629   mHitsMultInEvent  = new TH1D("hitsPerEvent","hitsPerEvent",100,0,100);
02630   mHitsMultPerTrack = new TH1D("hitsPerTrack","hitsPerTrack",10,0,10);
02631   mDaqOccupancy     = new TH1D("daqOccupancy","daqOccupancy",192,0,192);
02632   mDaqOccupancyValid= new TH1D("daqOccupancyValid","daqOccupancyValid",192,0,192);
02633   mDaqOccupancyProj = new TH1D("daqOccupancyProj","daqOccupancyProj",192,0,192);
02634   mHitsPosition     = new TH2D("hitsPosition","hitsPositions",300,-15.,15.,200,-5.,5.);
02635   
02636   mDaqOccupancyValidAll = new TH1D("daqOccupancyValidAll","",960,0,960);
02637   mDaqOccupancyProjAll = new TH1D("daqOccupancyProjAll","",960,0,960);
02638 
02639   mDaqOccupancyVpd = new TH1D("daqOccupancyVpd","",38,0,38);
02640   mDaqOccupancyValidVpd = new TH1D("daqOccupancyValidVpd","",38,0,38);
02641 
02642   // correlation
02643   for(int i=0;i<mNValidTrays_Run8;i++) {
02644     char hisname[100];
02645     sprintf(hisname,"Tray_%d",i+76);
02646     mHitCorr[i] = new TH2D(hisname,"",192,0,192,192,0,192);
02647     sprintf(hisname,"Tray_%d_module",i+76);
02648     mHitCorrModule[i] = new TH2D(hisname,"",32,0,32,32,0,32);
02649   }
02650   mHitCorrAll = new TH2D("Tray_All","",960,0,960,960,0,960);
02651 
02652   // TPC tracks
02653   if(mSaveTree) {
02654   mTrackTree = new TTree("track","track");
02655   mTrackTree->Branch("pt",&trackTree.pt,"pt/F");
02656   mTrackTree->Branch("eta",&trackTree.eta,"eta/F");
02657   mTrackTree->Branch("phi",&trackTree.phi,"phi/F");
02658   mTrackTree->Branch("nfitpts",&trackTree.nfitpts,"nfitpts/I");
02659   mTrackTree->Branch("dEdx",&trackTree.dEdx,"dEdx/F");
02660   mTrackTree->Branch("ndEdxpts",&trackTree.ndEdxpts,"ndEdxpts/I");
02661   mTrackTree->Branch("charge",&trackTree.charge,"charge/I");
02662   mTrackTree->Branch("projTrayId",&trackTree.projTrayId,"projTrayId/I");
02663   mTrackTree->Branch("projCellChan",&trackTree.projCellChan,"projCellChan/I");
02664   mTrackTree->Branch("projY",&trackTree.projY,"projY/F");
02665   mTrackTree->Branch("projZ",&trackTree.projZ,"projZ/F");
02666   }
02667 
02668   mTrackPtEta = new TH2D("trackPtEta","",100,0.,5.,60,-1.5,1.5);
02669   mTrackPtPhi = new TH2D("trackPtPhi","",100,0.,5.,120,0.,2*3.14159);
02670   mTrackNFitPts = new TH1D("trackNFitPts","",50,0.,50.);
02671   mTrackdEdxvsp = new TH2D("trackdEdxvsp","",500,0.,5.,1000,0.,10.);
02672   mNSigmaPivsPt = new TH2D("nSigmaPivsPt","",500,0.,5.,1000,-10.,10.);
02673 
02674   // primary association  
02675   mCellsPerEventMatch1 = new TH1D("cellsPerEventMatch1","cellPerEventMatch1",100,0,100);
02676   mHitsPerEventMatch1 = new TH1D("hitsPerEventMatch1","hitsPerEventMatch1",100,0,100);
02677   mCellsPerTrackMatch1 = new TH1D("cellsPerTrackMatch1","cellsPerTrackMatch1",100,0,100);
02678   mTracksPerCellMatch1 = new TH1D("tracksPerCellMatch1","tracksPerCellMatch1",100,0,100);
02679   mDaqOccupancyMatch1 = new TH1D("daqOccupancyMatch1","daqOccupancyMatch1",192,0,192);
02680   mDeltaHitMatch1 = new TH2D("deltaHitMatch1","deltaHitMatch1",300,-15,15,200,-5.,5.);
02681 
02682   // kick out multi-hit
02683   mCellsPerEventMatch2 = new TH1D("cellsPerEventMatch2","cellPerEventMatch2",100,0,100);
02684   mHitsPerEventMatch2 = new TH1D("hitsPerEventMatch2","hitsPerEventMatch2",100,0,100);
02685   mCellsPerTrackMatch2 = new TH1D("cellsPerTrackMatch2","cellsPerTrackMatch2",100,0,100);
02686   mTracksPerCellMatch2 = new TH1D("tracksPerCellMatch2","tracksPerCellMatch2",100,0,100);
02687   mDaqOccupancyMatch2 = new TH1D("daqOccupancyMatch2","daqOccupancyMatch2",192,0,192);
02688   mDeltaHitMatch2 = new TH2D("deltaHitMatch2","deltaHitMatch2",300,-15,15,200,-5.,5.);
02689 
02690   // sort out multi matched cells
02691   mCellsPerEventMatch3 = new TH1D("cellsPerEventMatch3","cellPerEventMatch3",100,0,100);
02692   mHitsPerEventMatch3 = new TH1D("hitsPerEventMatch3","hitsPerEventMatch3",100,0,100);
02693   mCellsPerTrackMatch3 = new TH1D("cellsPerTrackMatch3","cellsPerTrackMatch3",100,0,100);
02694   mTracksPerCellMatch3 = new TH1D("tracksPerCellMatch3","tracksPerCellMatch3",100,0,100);
02695   mDaqOccupancyMatch3 = new TH1D("daqOccupancyMatch3","daqOccupancyMatch3",192,0,192);
02696   mDeltaHitMatch3 = new TH2D("deltaHitMatch3","deltaHitMatch3",300,-15,15,200,-5.,5.);
02697 
02698   return;
02699 }
02700 
02701 
02702 //---------------------------------------------------------------------------
02703 // store histograms in a seperate root file
02704 void StTofrMatchMaker::writeHistogramsToFile(){
02705   // Output file
02706   TFile *theHistoFile =  new TFile(mHistoFileName.c_str(), "RECREATE");
02707   LOG_INFO  << "StTofrMatchMaker::writeHistogramsToFile()"
02708        << " histogram file " <<  mHistoFileName << endm;
02709 
02710   theHistoFile->cd();
02711 
02712   if(mHisto) {
02713 
02714     mEventCounterHisto->Write();
02715     mADCTDCCorelation->Write();
02716     mCellsMultInEvent->Write();
02717     mHitsMultInEvent->Write();
02718     mHitsMultPerTrack->Write();
02719     mDaqOccupancy->Write();
02720     mDaqOccupancyValid->Write();
02721     mDaqOccupancyProj->Write();
02722     mHitsPosition->Write();
02723 
02724     mDaqOccupancyValidAll->Write();
02725     mDaqOccupancyProjAll->Write();
02726     mDaqOccupancyVpd->Write();
02727     mDaqOccupancyValidVpd->Write();
02728     
02729     for(int i=0;i<mNValidTrays_Run8;i++) {
02730       mHitCorr[i]->Write();
02731       mHitCorrModule[i]->Write();
02732     }
02733     mHitCorrAll->Write();
02734 
02735     mTrackPtEta->Write();
02736     mTrackPtPhi->Write();
02737     mTrackNFitPts->Write();
02738     mTrackdEdxvsp->Write();
02739     mNSigmaPivsPt->Write();
02740 
02741     mCellsPerEventMatch1->Write();
02742     mHitsPerEventMatch1->Write();
02743     mCellsPerTrackMatch1->Write();
02744     mTracksPerCellMatch1->Write();
02745     mDaqOccupancyMatch1->Write();
02746     mDeltaHitMatch1->Write();
02747     
02748     mCellsPerEventMatch2->Write();
02749     mHitsPerEventMatch2->Write();
02750     mCellsPerTrackMatch2->Write();
02751     mTracksPerCellMatch2->Write();
02752     mDaqOccupancyMatch2->Write();
02753     mDeltaHitMatch2->Write();
02754 
02755     mCellsPerEventMatch3->Write();
02756     mHitsPerEventMatch3->Write();
02757     mCellsPerTrackMatch3->Write();
02758     mTracksPerCellMatch3->Write();
02759     mDaqOccupancyMatch3->Write();
02760     mDeltaHitMatch3->Write();
02761     
02762     theHistoFile->Write();  
02763     theHistoFile->Close(); 
02764    
02765   }
02766 
02767   if(mSaveTree) {
02768     TFile *theTreeFile =  new TFile((mHistoFileName+".tree.root").c_str(), "RECREATE");
02769     theTreeFile->cd();
02770     mTrackTree->Write();
02771     theTreeFile->Write();
02772     theTreeFile->Close();
02773   }
02774 
02775   return;
02776 }
02777 
02778 
02779 //---------------------------------------------------------------------------
02780 // determine pVPD event type (strobe or beam)
02781 bool StTofrMatchMaker::strobeEvent(StSPtrVecTofData& tofData){
02782   // determine strobe event from pVPD TDC data
02783 
02784   Int_t nStrobedPvpdTdcs=0;
02785   for(Int_t i=0;i<mNPVPD;i++)
02786     if((tofData[42+i]->tdc()>mStrobeTdcMin[i]) &&
02787        (tofData[42+i]->tdc()<mStrobeTdcMax[i]))
02788         nStrobedPvpdTdcs++;
02789   
02790   if (nStrobedPvpdTdcs==mNPVPD) return true;
02791 
02792   return false;
02793 }
02794 
02795 
02796 //---------------------------------------------------------------------------
02797 // determine whether this is a valid TOF beam event
02798 bool StTofrMatchMaker::validEvent(StEvent *event){
02799   mEventCounter++;
02800   // 1. must have non-zero pointer
02801   if (!event) return false;
02802   if(mHisto) mEventCounterHisto->Fill(1);
02803 
02804   // 2. must have a valid primary vertex 
02805 //  if (!event->primaryVertex()) return false;
02806   mAcceptedEventCounter++;
02807   if(mHisto) mEventCounterHisto->Fill(2);
02808 
02809   // 3a. must have TOF collection
02810   if (!event->tofCollection()){
02811     LOG_INFO  << "TOF is not present" << endm;
02812     return false;
02813   }
02814   if(mHisto) mEventCounterHisto->Fill(3);
02815 
02816   // 3b. must have TOF raw data available
02817   if (!(event->tofCollection()->dataPresent()||event->tofCollection()->rawdataPresent())){
02818     LOG_INFO  << "TOF is present but no Raw Data" << endm;
02819     if  (!(event->tofCollection()->cellsPresent())){
02820       LOG_INFO  << "              and no Cell Data" << endm;
02821     }
02822     return false;
02823   }
02824   mTofEventCounter++;
02825   if(mHisto) mEventCounterHisto->Fill(4);
02826 
02827   
02828   
02829   // 4. must be a TOF beam event, i.e. a non-strobe event
02830   if(event->tofCollection()->dataPresent()) {
02831     StSPtrVecTofData  &tofData = event->tofCollection()->tofData();
02832     LOG_INFO  << " tofData size = " << tofData.size() << endm;
02833     if (strobeEvent(tofData)){
02834       mTofStrobeEventCounter++;
02835       if (event->primaryVertex()) mAcceptAndStrobe++; // keep track of #valid strobed evts
02836       gMessMgr->Info("strobe event","OTS");
02837       return false;
02838     }
02839   }
02840   if(mHisto) mEventCounterHisto->Fill(5);
02841   
02842   mAcceptAndBeam++;
02843 
02844   // and we have a winner!
02845   gMessMgr->Info("TOF present ... and valid beam event","OTS");
02846 
02847   return true;
02848 }
02849 
02850 
02851 //---------------------------------------------------------------------------
02852 // determine whether this is a valid TPC track
02853 bool StTofrMatchMaker::validTrack(StTrack *track){
02854   // 1. no track, no go.
02855   if (!track) return false;
02856 
02857   // 2. track quality flag, should be >0
02858   if (track->flag()<=0) return false;
02859 
02860   // 3. minimum #hits per track
02861   if (track->topologyMap().numberOfHits(kTpcId) < mMinHitsPerTrack) return false;
02862   // 4. minimum #fit points per track
02863   if (track->fitTraits().numberOfFitPoints(kTpcId) < mMinFitPointsPerTrack) return false;
02864   // 5. minimum #fit points over #maximum points
02865   float ratio = (1.0*track->fitTraits().numberOfFitPoints(kTpcId)) / (1.0*track->numberOfPossiblePoints(kTpcId));
02866   if (ratio < mMinFitPointsOverMax) return false;
02867 
02868   return true;
02869 }
02870 
02871 //---------------------------------------------------------------------------
02872 // determine whether this is a valid TPC track
02873 bool StTofrMatchMaker::validTrackRun8(StGlobalTrack *track){
02874   // 1. no track, no go.
02875   if (!track) return false;
02876 
02877   // 2. track quality flag, should be >0
02878   if (track->flag()<=0) return false;
02879 
02880   // 3. minimum #hits per track
02881   if (track->topologyMap().numberOfHits(kTpcId) < mMinHitsPerTrack) return false;
02882   // 4. minimum #fit points per track
02883   if (track->fitTraits().numberOfFitPoints(kTpcId) < mMinFitPointsPerTrack) return false;
02884   // 5. chi2 per track
02885   float chi2 = track->fitTraits().chi2(0);
02886   if ( chi2 > 10. ) return false;
02887   
02888   // 6. error
02889   StDcaGeometry *dcaGeom = track->dcaGeometry();
02890   if(!dcaGeom) return false;
02891   float pt = dcaGeom->pt();
02892   float sigma_pT = sqrt(dcaGeom->errMatrix()[9])*pt*pt;
02893   if ( sigma_pT/pt > 0.3 ) return false;
02894 
02895   return true;
02896 }
02897 
02898 
02899 //---------------------------------------------------------------------------
02900 // determine whether this is a valid TOF track
02901 bool StTofrMatchMaker::validTofTrack(StTrack *track){
02902   // select valid tracks for time-of-flight calculations
02903 
02904   // 1. track must exist
02905   if (!track) return false;
02906 
02907   // 2. track quality flag, should be >0
02908   if (track->flag()<=0) return false;
02909 
02910   // 3. track must be a primary track
02911   if (!dynamic_cast<StPrimaryTrack*>(track)) return false;
02912 
02913   // 4. DCA cut (obsolete?)
02914   Double_t DCA= track->impactParameter();
02915   Int_t charge = track->geometry()->charge();
02916   if (DCA > mMaxDCA) {LOG_INFO  << "dca>max:" << DCA<< endm; return false;}
02917   if (charge==0) { LOG_INFO  << " neutral charge" << endm; return false; }
02918 
02919   return true;
02920 }  
02921 
02922 
02923 //---------------------------------------------------------------------------
02924 // returns the proper track geometry, based on a global user setting
02925 StTrackGeometry* StTofrMatchMaker::trackGeometry(StTrack* track){
02926   // returns apropriate StTrackGeometry (standard or outerGeometry)
02927   if (!track) return 0;
02928   StTrackGeometry *thisTrackGeometry;
02929   if (mOuterTrackGeometry)
02930     thisTrackGeometry = track->outerGeometry();
02931   else
02932     thisTrackGeometry = track->geometry();
02933   return thisTrackGeometry;
02934 }