StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StTofrMatchMaker.cxx
1 /*******************************************************************
2  *
3  * $Id: StTofrMatchMaker.cxx,v 1.34 2018/02/26 23:26:52 smirnovd Exp $
4  *
5  * Author: Xin Dong
6  *****************************************************************
7  *
8  * Description: Tofr Match Maker to do the matching between the
9  * fired celles and TPC tracks ( similar to Frank's
10  * TofpMatchMaker )
11  *
12  *******************************************************************/
13 #include <iostream>
14 #include "StEventTypes.h"
15 #include "Stypes.h"
16 #include "StThreeVectorD.hh"
17 #include "StThreeVectorF.hh"
18 #include "StHelix.hh"
19 #include "StTrackGeometry.h"
20 #include "StDcaGeometry.h"
21 #include "StDedxPidTraits.h"
22 #include "StTrackPidTraits.h"
23 #include "StarClassLibrary/StParticleTypes.hh"
24 #include "StarClassLibrary/StParticleDefinition.hh"
25 #include "StTpcDedxPidAlgorithm.h"
26 #include "StEventUtilities/StuRefMult.hh"
27 #include "PhysicalConstants.h"
28 #include "StPhysicalHelixD.hh"
29 #include "StHelix.hh"
30 #include "StTofUtil/tofPathLength.hh"
31 #include "StTofUtil/StTofrDaqMap.h"
32 #include "StTofUtil/StTofINLCorr.h"
33 #include "StTofUtil/StTofrGeometry.h"
34 #include "StTofUtil/StTofCellCollection.h"
35 #include "tables/St_pvpdStrobeDef_Table.h"
36 #include "tables/St_tofPedestal_Table.h"
37 #include "tables/St_tofConfig_Table.h"
38 #include "tables/St_tofTrayConfig_Table.h"
39 //#include "tables/St_tofr5INLtable_Table.h"
40 #include "StTofUtil/tofPathLength.hh"
41 #include "TH1.h"
42 #include "TH2.h"
43 #include "TFile.h"
44 #include "TTree.h"
45 #include "StMessMgr.h"
46 #include "StMemoryInfo.hh"
47 #include "StTimer.hh"
48 #include "StTofrMatchMaker.h"
49 //#include "TMemStat.h"
50 
51 
52 // Define the static constants:
53 
54 const Int_t StTofrMatchMaker::mDAQOVERFLOW = 255;
55 const Int_t StTofrMatchMaker::mNTOFP = 41;
56 const Int_t StTofrMatchMaker::mNPVPD = 6;
57 const Int_t StTofrMatchMaker::mNTOFR = 120;
58 const Int_t StTofrMatchMaker::mNTOFR5 = 192;
59 
60 const Int_t StTofrMatchMaker::mNTOF = 192; // 192 for tof in Run 8++
61 const Int_t StTofrMatchMaker::mNModule = 32; // 32 for tofr5++
62 const Int_t StTofrMatchMaker::mNCell = 6;
63 const Int_t StTofrMatchMaker::mNVPD = 19; //
64 
65 const Int_t StTofrMatchMaker::mEastVpdTrayId = 122;
66 const Int_t StTofrMatchMaker::mWestVpdTrayId = 121;
67 
68 const Int_t StTofrMatchMaker::mNValidTrays_Run3 = 1;
69 const Int_t StTofrMatchMaker::mNValidTrays_Run4 = 1;
70 const Int_t StTofrMatchMaker::mNValidTrays_Run5 = 1;
71 const Int_t StTofrMatchMaker::mNValidTrays_Run6 = 0;
72 const Int_t StTofrMatchMaker::mNValidTrays_Run7 = 0;
73 const Int_t StTofrMatchMaker::mNValidTrays_Run8 = 5;
74 
75 const Int_t StTofrMatchMaker::mTdigBoard = 10;
76 const Int_t StTofrMatchMaker::mTdcOnBoard = 4;
77 const Int_t StTofrMatchMaker::mTdcChannel = 1024;
78 
79 //---------------------------------------------------------------------------
80 StTofrMatchMaker::StTofrMatchMaker(const Char_t *name): StMaker(name)
81  , mTofrAdc(mNTOFR,0)
82  , mTofrTdc(mNTOFR,0)
83  , mPvpdAdc(mNPVPD,0)
84  , mPvpdAdcLoRes( mNPVPD,0)
85  , mPvpdTdc(mNPVPD,0)
86  , mStrobeTdcMin(mNPVPD,0)
87  , mStrobeTdcMax(mNPVPD,0)
88  , mPedTOFr(mNTOFR,0)
89  //year 5
90  , mPvpdToT(mNPVPD,0)
91  , mTofr5Tdc(mNTOFR5,0)
92  , mTofr5ToT(mNTOFR5,0) // ToT as adc
93  , mHitCorr(mNValidTrays_Run8,(TH2D*)0)
94  , mHitCorrModule(mNValidTrays_Run8,(TH2D*)0)
95 {
96 
97  // set default values
98  mEventCounter = 0;
99  mAcceptedEventCounter = 0;
100  mTofEventCounter = 0;
101  mTofStrobeEventCounter = 0;
102  mAcceptAndStrobe = 0;
103  mAcceptAndBeam = 0;
104 
105  mTofrGeom = 0;
106  mDaqMap = 0;
107  mTofINLCorr = 0;
108 
109  mWidthPad = 3.45;
110 
111  setValidAdcRange(1,1200);
112  setValidTdcRange(1,2047);
113  setOuterTrackGeometry();
114  setMinHitsPerTrack(15);
115  setMinFitPointsPerTrack(15);
116  setMinFitPointsOverMax(0.52);
117  setMaxDCA(9999.);
118 
119  setCreateHistoFlag(kFALSE);
120  setHistoFileName("tofana.root");
121  setCreateTreeFlag(kFALSE);
122  setSaveGeometry(kFALSE);
123  doPrintMemoryInfo = kFALSE;
124  doPrintCpuInfo = kFALSE;
125 }
126 
127 StTofrMatchMaker::~StTofrMatchMaker(){ /* nope */}
128 
129 //void StTofrMatchMaker::Clear(Option_t *opt){StMaker::Clear();}
130 
131 //---------------------------------------------------------------------------
132 Int_t StTofrMatchMaker::Init(){
133  gMessMgr->Info("StTofrMatchMaker -- initializing ...","OS");
134  LOG_INFO << "Valid TDC range: " << mMinValidTdc << " " << mMaxValidTdc << endm;
135  LOG_INFO << "Valid ADC range: " << mMinValidAdc << " " << mMaxValidAdc << endm;
136  LOG_INFO << "Minimum hits per track: " << mMinHitsPerTrack << endm;
137  LOG_INFO << "Minimum fitpoints per track: " << mMinFitPointsPerTrack << endm;
138  LOG_INFO << "Maximum DCA: " << mMaxDCA << endm;
139  if (!mOuterTrackGeometry)
140  gMessMgr->Warning("Warning: using standard trackgeometry()","OS");
141 
142  // m_Mode can be set by SetMode() method
143  if(m_Mode) {
144 // setHistoFileName("tofana.root");
145  } else {
146  setHistoFileName("");
147  }
148 
149  if (mHisto){
150  bookHistograms();
151  LOG_INFO << "Histograms are booked" << endm;
152  if (mHistoFileName!="")
153  LOG_INFO << "Histograms will be stored in " << mHistoFileName.c_str() << endm;
154  }
155 
156  // reset event counters
157  mEventCounter = 0;
158  mAcceptedEventCounter = 0;
159  mTofEventCounter = 0;
160  mTofStrobeEventCounter = 0;
161  mAcceptAndStrobe = 0;
162  mAcceptAndBeam = 0;
163 
164  return kStOK;
165 }
166 
167 
168 //---------------------------------------------------------------------------
169 Int_t StTofrMatchMaker::InitRun(Int_t runnumber){
170 
171  // determine TOF configuration from run#
172  mYear2 = (runnumber<4000000);
173  mYear3 = (runnumber>4000000&&runnumber<5000000);
174  mYear4 = (runnumber>5000000&&runnumber<6000000);
175  mYear5 = (runnumber>6000000&&runnumber<7000000);
176  mYear8 = (runnumber>9000000&&runnumber<10000000);
177  mYearX = (runnumber>10000000);
178 
179  if (mYearX){
180  Error(":InitRun","Wrong BFC configuration for run %d. Use StBTofHitMaker for Run9+ data.",runnumber);
181  return 0;
182  }
183 
184  gMessMgr->Info("StTofrMatchMaker -- Initializing TofGeometry (InitRun)","OS");
186  // TOFr geometry initializtion -- from GEANT geometry directly
187  // need St_geant_Maker be loaded before
189  mTofrGeom = new StTofrGeometry("tofrGeom","tofrGeom in MatchMaker");
190  if(!mTofrGeom->IsInitDone()) {
191  gMessMgr->Info("TofrGemetry initialization..." ,"OS");
192  TVolume *starHall = (TVolume *)GetDataSet("HALL");
193  mTofrGeom->Init(starHall);
194  }
195  // other makers can get this geometry
196  if(mGeometrySave) {
197  if(TDataSet *geom = GetDataSet("tofrGeometry")) delete geom;
198  AddConst(new TObjectSet("tofrGeometry",mTofrGeom));
199  }
200 
201  if(Debug()) {
202  gMessMgr->Info(" Test the TofrGeometry","OS");
203  // Test
204  TVolumeView* topNode = mTofrGeom->GetTopNode();
205  if(!topNode) {
206  gMessMgr->Warning("Warning: no Top node!","OS");
207  return kStWarn;
208  }
209  mTofrGeom->Print();
210 
211  topNode->ls(9);
212  // TShape *topshape = topNode->GetShape();
213  LOG_INFO << " # of trays = " << topNode->GetListSize() << endm;
214  TList *list = topNode->Nodes();
215  Int_t ibtoh =0;
216  TVolumeView *sectorVolume = 0;
217  for(Int_t i=0;i<list->GetSize();i++) {
218  sectorVolume = dynamic_cast<TVolumeView*> (list->At(i));
219  if ( i>=60 ) ibtoh = 1;
220  LOG_INFO << " test sector size = " <<sectorVolume->GetListSize() << endm;
221  StTofrGeomTray *tray = new StTofrGeomTray(ibtoh, sectorVolume, topNode);
222  tray->Print();
223  TVolumeView *trayVolume = tray->GetfView();
224  TList *list1 = trayVolume->Nodes();
225  LOG_INFO << " # of modules in tray " << tray->Index() << " = " << trayVolume->GetListSize() << endm;
226  if (!list1 ) continue;
227  TVolumeView *sensorVolume = 0;
228  for(Int_t j=0;j<list1->GetSize();j++) {
229  sensorVolume = dynamic_cast<TVolumeView*> (list1->At(j));
230  StTofrGeomSensor *sensor = new StTofrGeomSensor(sensorVolume, topNode);
231  sensor->Print();
232  delete sensor;
233  }
234  delete tray;
235  }
236  }
237 
239  // TOFr Daq map initialization -- load from StTofUtil
241  mDaqMap = new StTofrDaqMap();
242  mTofINLCorr = new StTofINLCorr();
243  if(mYear2||mYear3||mYear4) {
244  if(mYear3) mDaqMap->setNValidTrays(mNValidTrays_Run3);
245  if(mYear4) mDaqMap->setNValidTrays(mNValidTrays_Run4);
246  mDaqMap->init(this);
247  // AddConst(new TObjectSet("tofrDaqMap",mDaqMap));
248  LOG_INFO << " Initialize Daq map for run 2,3,4 ... " << endm;
249 
250  if(Debug()) {
251  for(Int_t i=0;i<mNTOFR;i++) {
252  LOG_INFO << i << " -- adc=" << mDaqMap->DaqChan2ADCChan(i) << " -- tdc=" << mDaqMap->DaqChan2TDCChan(i) << endm;
253  IntVec map = mDaqMap->DaqChan2Cell(i);
254  LOG_INFO << " tray=" << map[0] << " module=" << map[1] << " cell=" << map[2] << endm;
255  }
256  }
257 
259  // Load configuration parameters from dbase
260  // need "[shell] setenv Calibrations_tof reconV0"
262  gMessMgr->Info(" -- retrieving run parameters from Calibrations_tof","OS");
263  TDataSet *mDbDataSet = GetDataBase("Calibrations/tof/pvpdStrobeDef");
264  if (!mDbDataSet){
265  gMessMgr->Error("unable to get TOF pvpdStrobeDef table","OS");
266  // assert(mDbDataSet);
267  return kStErr;
268  }
269  St_pvpdStrobeDef* pvpdStrobeDef = static_cast<St_pvpdStrobeDef*>(mDbDataSet->Find("pvpdStrobeDef"));
270  if (!pvpdStrobeDef){
271  gMessMgr->Error("unable to find pVPD strobe def parameters","OS");
272  // assert(pvpdStrobeDef);
273  return kStErr;
274  }
275  pvpdStrobeDef_st *strobeDef = static_cast<pvpdStrobeDef_st*>(pvpdStrobeDef->GetArray());
276  Int_t numRows = pvpdStrobeDef->GetNRows();
277  if (mNPVPD != numRows) gMessMgr->Warning("#tubes inconsistency in dbase");
278  for (Int_t i=0;i<mNPVPD;i++){
279  Int_t ii = strobeDef[i].id - 1;
280  mStrobeTdcMin[ii] = strobeDef[i].strobeTdcMin;
281  mStrobeTdcMax[ii] = strobeDef[i].strobeTdcMax;
282  if (Debug())
283  LOG_INFO << "tube " << strobeDef[i].id << " min:"<< strobeDef[i].strobeTdcMin
284  <<" max:"<< strobeDef[i].strobeTdcMax<< endm;
285  }
286 
287  //TOFr pedestals
288  mDbDataSet = GetDataBase("Calibrations/tof/tofPedestal");
289  if (!mDbDataSet){
290  gMessMgr->Error("unable to access TOF tofPedestal table","OS");
291  // assert(mDbDataSet);
292  return kStErr;
293  }
294  St_tofPedestal* tofPed = static_cast<St_tofPedestal*>(mDbDataSet->Find("tofPedestal"));
295  if (!tofPed){
296  gMessMgr->Error("unable to find TOF pedestal table","OS");
297  return kStErr;
298  // assert(tofPed);
299  }
300  tofPedestal_st *pedestal = static_cast<tofPedestal_st*>(tofPed->GetArray());
301  Int_t entries = (Int_t)(pedestal[0].entries);
302  if (mNTOFR>entries) gMessMgr->Warning("#TOFr channels inconsistency in dbase");
303  for (Int_t i=0;i<entries;i++) {
304  Int_t ii = pedestal[0].daqChannel[i];
305  Int_t adcChan = pedestal[0].adcChan[i];
306  if(adcChan>=60)
307  mPedTOFr[ii] = (Double_t)(pedestal[0].adcPedestal[i]);
308  if (Debug())
309  LOG_INFO << "daqChannel" << ii << " ped:" << pedestal[0].adcPedestal[i] <<endm;
310  }
311 
312  /* only test
313  //TOFr configuration
314  St_tofTrayConfig* toftrayconfig = static_cast<St_tofTrayConfig*>(mDbDataSet->Find("tofTrayConfig"));
315  if(!toftrayconfig) {
316  gMessMgr->Error("unable to find TOF tray configuration parameters","OS");
317  // assert(toftrayconfig);
318  return kStErr;
319  }
320  tofTrayConfig_st *toftrayconf = static_cast<tofTrayConfig_st*>(toftrayconfig->GetArray());
321  entries = toftrayconf[0].entries;
322  for(Int_t i=0;i<entries;i++) {
323  Int_t iTray = toftrayconf[0].iTray[i];
324  Int_t nModules = toftrayconf[0].nModules[i];
325  if(Debug())
326  LOG_INFO << " " << iTray << " tray has " << nModules << " modules " << endm;
327  }
328  */
329 
330  } else if(mYear5) {
331  mDaqMap->setNValidTrays(mNValidTrays_Run5);
332  mDaqMap->initFromDbaseY5(this);
333  LOG_INFO << " Initialize Daq map for run 5 ... " << endm;
334 
335  if(Debug()) {
336  // if(0) {
337  for(Int_t i=0;i<mNTOFR5;i++) {
338  IntVec map = mDaqMap->Tofr5TDCChan2Cell(i);
339  LOG_INFO << "InitRun():i="<<i<<" tray=" << map[0] << " module=" << map[1] << " cell=" << map[2] << endm;
340  }
341  }
342  } else if(mYear8) {
343  mDaqMap->setNValidTrays(mNValidTrays_Run8);
344  mDaqMap->initFromDbaseGeneral(this);
345  LOG_INFO << " Initialize Daq map for run 8 ... " << endm;
346 
347  mTofINLCorr->setNValidTrays(mNValidTrays_Run8);
348  mTofINLCorr->initFromDbase(this);
349  LOG_INFO << " Initialize INL table for run 8 ... " << endm;
350 
351  }
352 
353  return kStOK;
354 }
355 
356 //----------------------------------------------------------------------------
357 Int_t StTofrMatchMaker::FinishRun(Int_t runnumber){
358 
359  LOG_INFO << "StTofrMatchMaker -- cleaning up geometry (FinishRun)" << endm;
360  if (mTofrGeom) delete mTofrGeom;
361  mTofrGeom=0;
362 
363  if(mDaqMap) delete mDaqMap;
364  mDaqMap = 0;
365 
366  if(mTofINLCorr) delete mTofINLCorr;
367  mTofINLCorr = 0;
368 
369  return kStOK;
370 }
371 
372 
373 //---------------------------------------------------------------------------
374 Int_t StTofrMatchMaker::Finish(){
375 
376  LOG_INFO << "StTofrMatchMaker ----- RUN SUMMARY ----- (Finish)\n"
377  << "\tProcessed " << mEventCounter << " events."
378  << " Accepted " << mAcceptedEventCounter << " events."
379  << " Rejected " << mEventCounter - mAcceptedEventCounter << " events\n"
380  << "\tTOF events " << mTofEventCounter
381  << ". Beam " << mTofEventCounter - mTofStrobeEventCounter
382  << " Strobe " << mTofStrobeEventCounter
383  << "\n\t Accept & Strobe " << mAcceptAndStrobe << " events\n"
384  << "\t Accept & Beam " << mAcceptAndBeam << " events" << endm;
385 
386  //if (mHisto) writeHistogramsToFile();
387  if (mHistoFileName!="") writeHistogramsToFile();
388  return kStOK;
389 }
390 
391 
392 //---------------------------------------------------------------------------
393 Int_t StTofrMatchMaker::Make(){
394  gMessMgr->Info("StTofrMatchMaker -- welcome","OS");
395 
396  Int_t iret = kStOK;
397  if(mYear2||mYear3||mYear4) {
398  iret = processEventYear2to4();
399  } else if(mYear5) {
400  iret = processEventYear5();
401  } else if(mYear8) {
402  iret = processEventYear8();
403  }
404  return iret;
405 }
406 
407 //---------------------------------------------------------------------------
408 Int_t StTofrMatchMaker::processEventYear2to4(){
409  if(mHisto) mEventCounterHisto->Fill(0);
410  // event selection ...
411  mEvent = (StEvent *) GetInputDS("StEvent");
412  if (!validEvent(mEvent)){
413  gMessMgr->Info("StTofrMatchMaker -- nothing to do ... bye-bye","OS");
414  return kStOK;
415  }
416 
417  if (mYear2) {
418  gMessMgr->Info("StTofrMatchMaker -- no TOFr in Year2!","OS");
419  return kStOK;
420  }
421 
422  // timing & memory info -only when requested-
423  StTimer timer;
424  if (doPrintCpuInfo) timer.start();
425  if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
426 
427  //.........................................................................
428  // check for tofCollection and fill local copy with ADC and TDC data
429  StTofCollection *theTof = mEvent->tofCollection();
430  getTofData(theTof);
431 
432  //.........................................................................
433  // update Tofr cell counters
434  Float_t sumAdcTofr(0.); Int_t nAdcTofr(0), nTdcTofr(0), nAdcTdcTofr(0);
435  for (Int_t i=0;i<mNTOFR;i++){
436  sumAdcTofr += mTofrAdc[i];
437  bool tdcValid = validTdc(mTofrTdc[i]);
438  bool adcValid = validAdc(mTofrAdc[i]);
439  if (tdcValid) nTdcTofr++;
440  if (adcValid) nAdcTofr++;
441  if (adcValid && tdcValid) nAdcTdcTofr++;
442  }
443  if (Debug())
444  LOG_INFO << " Tofr #Adc:" << nAdcTofr << " #Tdc:" << nTdcTofr << endm;
445 
446  //..................................
447  // update pVPD tubes counters
448  Float_t sumAdcPvpd=0; Int_t nAdcPvpd=0, nTdcPvpd=0;
449  for (Int_t i=0;i<mNPVPD;i++){
450  sumAdcPvpd += mPvpdAdc[i];
451  bool tdcValid = validTdc(mPvpdTdc[i]);
452  bool adcValid = validAdc(mPvpdAdc[i]);
453  if (tdcValid) nTdcPvpd++;
454  if (adcValid) nAdcPvpd++;
455  }
456  if (Debug())
457  LOG_INFO << " pVPD #Adc:" << nAdcPvpd << " #Tdc:" << nTdcPvpd << endm;
458 
459 
460 
461  // number of primary tracks
462  // (note: different meaning in event.root and richtof.root)
463  Int_t refmult(0);
464  bool richTofMuDST = (mEvent->summary()->numberOfExoticTracks() == -999);
465  if (richTofMuDST)
466  refmult = mEvent->summary()->numberOfGoodTracks();
467  else
468  refmult = uncorrectedNumberOfPrimaries(*mEvent);
469 
470  if (Debug()){
471  LOG_INFO << " #Tracks :" << mEvent->summary()->numberOfTracks()
472  << "\n #goodPrimaryTracks:" << mEvent->summary()->numberOfGoodPrimaryTracks()
473  << "\n #uncorr.prim.tracks :" << refmult << endm;
474  if (!richTofMuDST)
475  LOG_INFO << " #goodTracks (global):" << mEvent->summary()->numberOfGoodTracks() << endm;
476  }
477 
478  //.........................................................................
479  // A. build vector of candidate cells with valid ADC signals
480  // idVector validCellIdVec;
481  tofCellHitVector daqCellsHitVec;
482  // daqCellsHitVec.clear();
483  idVector validModuleVec;
484 
485  for (Int_t i=0;i<mNTOFR;i++){
486  Float_t rawAdc = mTofrAdc[i];
487  Float_t rawTdc = mTofrTdc[i];
488  if (mHisto) mDaqOccupancy->Fill(i);
489 
490  if (validAdc(rawAdc) && validTdc(rawTdc) && rawAdc>mPedTOFr[i] ){
491  IntVec map = mDaqMap->DaqChan2Cell(i);
492  Int_t IDtray = map[0];
493  Int_t IDmodule = map[1];
494  Int_t IDcell = map[2];
495 
496  StructCellHit aDaqCellHit;
497  aDaqCellHit.channel = i;
498  aDaqCellHit.tray = IDtray;
499  aDaqCellHit.module = IDmodule;
500  aDaqCellHit.cell = IDcell;
501  daqCellsHitVec.push_back(aDaqCellHit);
502  // additional valid number configuration
503  int id = IDtray*100+IDmodule;
504  bool ifind = kFALSE;
505  for(size_t iv=0;iv<validModuleVec.size();iv++) {
506  if(id==validModuleVec[iv]) {
507  ifind = kTRUE;
508  break;
509  }
510  }
511  if(!ifind) validModuleVec.push_back(id);
512 
513  if (mHisto) {
514  mDaqOccupancyValid->Fill(i);
515  Int_t adcchan = mDaqMap->DaqChan2ADCChan(i);
516  Int_t tdcchan = mDaqMap->DaqChan2TDCChan(i);
517  mADCTDCCorelation->Fill(tdcchan, adcchan);
518  }
519  if(Debug()) {
520  LOG_INFO <<"A: daqId=" << i << " rawAdc= " << rawAdc << " rawTdc="<< rawTdc <<endm;
521  }
522  }
523  }
524  // end of Sect.A
525  if(Debug()) {
526  LOG_INFO << " total # of cells = " << daqCellsHitVec.size() << endm;
527  for(size_t iv = 0;iv<validModuleVec.size();iv++) {
528  LOG_INFO << " module # " << validModuleVec[iv] << " Valid! " << endm;
529  }
530  }
531  if(mHisto) {
532  mCellsMultInEvent->Fill(daqCellsHitVec.size());
533  if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
534  }
535  if(!daqCellsHitVec.size()) return kStOK;
536 
537  //.........................................................................
538  // B. loop over global tracks and determine all cell-track matches
539  //
540  tofCellHitVector allCellsHitVec;
541  // allCellsHitVec.clear();
542  StructCellHit cellHit;
543 
544  StSPtrVecTrackNode& nodes = mEvent->trackNodes();
545  Int_t nAllTracks=0;
546  for (unsigned int iNode=0; iNode<nodes.size(); iNode++){
547  tofCellHitVector cellHitVec;
548  // cellHitVec.clear();
549  StTrack *theTrack = nodes[iNode]->track(global);
550 
551  // make sure we have a track, a miniDST might have removed it...
552  if (validTrack(theTrack)){
553  nAllTracks++;
554  StPhysicalHelixD theHelix = trackGeometry(theTrack)->helix();
555 
556  IntVec projTrayVec;
557  if(!mTofrGeom->projTrayVector(theHelix, projTrayVec)) continue;
558 
559  IntVec idVec;
560  DoubleVec pathVec;
561  PointVec crossVec;
562 
563 // idVec.clear();
564 // pathVec.clear();
565 // crossVec.clear();
566 
567  Int_t ncells = 0;
568  // if(mTofrGeom->HelixCrossCellIds(theHelix,idVec,pathVec,crossVec) ) {
569  if(mTofrGeom->HelixCrossCellIds(theHelix, validModuleVec, projTrayVec, idVec, pathVec, crossVec)) {
570  Int_t cells = idVec.size();
571  for (Int_t i=0; i<cells; i++) {
572  Int_t icell,imodule,itray;
573  Double_t local[3],global[3];
574  for(Int_t i2=0;i2<3;i2++){
575  local[i2]=0;
576  }
577  global[0]=crossVec[i].x();
578  global[1]=crossVec[i].y();
579  global[2]=crossVec[i].z();
580  mTofrGeom->DecodeCellId(idVec[i], icell, imodule, itray);
581  // LOG_INFO << " decode " << idVec[i] << " to tray#" << itray << " module#" << imodule << " cell#" << icell << endm;
582  StTofrGeomSensor* sensor =
583  mTofrGeom->GetGeomSensor(imodule,itray);
584  if(!sensor) {
585  gMessMgr->Warning("","OS") << " No sensitive module in the projection??? -- Something weird!!! " << endm;
586  continue;
587  }
588  sensor->Master2Local(&global[0],&local[0]);
589  icell = sensor->FindCellIndex(local);
590  // StThreeVectorD glo=sensor->GetCenterPosition();
591  StThreeVectorD glo(global[0], global[1], global[2]);
592  StThreeVectorD hitPos(local[0], local[1], local[2]);
593  delete sensor;
594  if (local[2]<=3.4&&local[2]>=-2.7) {
595  Int_t Iarray = mDaqMap->Cell2DaqChan(itray, imodule, icell);
596  if(Iarray>=mDAQOVERFLOW) continue;
597  ncells++;
598  cellHit.channel = Iarray;
599  cellHit.tray = itray;
600  cellHit.module = imodule;
601  cellHit.cell = icell;
602  cellHit.trackIdVec.push_back(iNode);
603  cellHit.hitPosition = glo; // global position
604  cellHit.zhit = (Float_t)hitPos.z();
605  cellHit.yhit = (Float_t)hitPos.y();
606  cellHitVec.push_back(cellHit);
607  allCellsHitVec.push_back(cellHit);
608  if(mHisto) {
609  mDaqOccupancyProj->Fill(Iarray);
610  mHitsPosition->Fill(hitPos.y(), hitPos.z());
611  }
612 
613  if(Debug()) {
614  LOG_INFO <<"B: nodeid=" << iNode << " projected in " << " tray="<< itray << " module="<<imodule<<" cell="<<icell<<endm;
615  LOG_INFO <<" hit position " << hitPos << endm;
616  }
617  }
618  } // for (Int_t i=0...)
619  } // endif(helixcross...)
620  if(ncells>0&&mHisto) mHitsMultPerTrack->Fill(ncells);
621 
622  } // if(ValidTrack)..
623  } // loop over nodes
624  if(Debug())
625  LOG_INFO << "B: matched/available/total #tracknodes: " <<allCellsHitVec.size() << "/" <<nAllTracks << "/" << nodes.size() << endm;
626  if(mHisto) {
627  mHitsMultInEvent->Fill(allCellsHitVec.size());
628  if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
629  }
630  // end of Sect.B
631 
632  //.........................................................................
633  // C. Match find Neighbours -- identify crosstalk
634  //
635  tofCellHitVector matchHitCellsVec;
636  // matchHitCellsVec.clear();
637 
638  tofCellHitVectorIter daqIter = daqCellsHitVec.begin();
639  for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
640  tofCellHitVectorIter proIter = allCellsHitVec.begin();
641  for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
642  if( (daqIter->tray==proIter->tray)&&
643  (daqIter->module==proIter->module) &&
644  ( ( (proIter->cell==6)&&((proIter->cell==daqIter->cell) ||
645  (proIter->cell==daqIter->cell+1)) )
646  || ( (proIter->cell==1)&&((proIter->cell==daqIter->cell) ||
647  (proIter->cell==daqIter->cell-1)) )
648  || ( (proIter->cell>=2&&proIter->cell<=6) &&
649  ( (proIter->cell==daqIter->cell) ||
650  (proIter->cell==daqIter->cell-1) ||
651  (proIter->cell==daqIter->cell+1) ) ) ) ) {
652  cellHit.channel = daqIter->channel;
653  cellHit.tray = daqIter->tray;
654  cellHit.module = daqIter->module;
655  cellHit.cell = daqIter->cell;
656  cellHit.hitPosition = proIter->hitPosition;
657  cellHit.trackIdVec = proIter->trackIdVec;
658  cellHit.zhit = proIter->zhit;
659  cellHit.yhit = proIter->yhit;
660  matchHitCellsVec.push_back(cellHit);
661  }
662  }
663  } //end {sec. C}
664  if(Debug()) {
665  LOG_INFO << "C: before/after: " << allCellsHitVec.size() << "/" << matchHitCellsVec.size() << endm;
666  }
667  if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
668 
669  //.........................................................................
670  // D. sort hit vectors and deal with (discard) cells matched by multiple tracks
671  //
672  Int_t nSingleHitCells(0);
673  Int_t nMultiHitsCells(0);
674 
675  tofCellHitVector singleHitCellsVec;
676  tofCellHitVector multiHitsCellsVec;
677 // singleHitCellsVec.clear();
678 // multiHitsCellsVec.clear();
679 
680  tofCellHitVector tempVec = matchHitCellsVec;
681  tofCellHitVector erasedVec = tempVec;
682  while (tempVec.size() != 0) {
683  Int_t nTracks = 0;
684  idVector trackIdVec;
685 
686  tofCellHitVectorIter tempIter=tempVec.begin();
687  tofCellHitVectorIter erasedIter=erasedVec.begin();
688  while(erasedIter!= erasedVec.end()) {
689  if(tempIter->tray == erasedIter->tray &&
690  tempIter->module == erasedIter->module &&
691  tempIter->cell == erasedIter->cell) {
692  nTracks++;
693  trackIdVec.push_back(erasedIter->trackIdVec.back()); // merge
694  erasedVec.erase(erasedIter);
695  erasedIter--;
696  }
697  erasedIter++;
698  }
699 
700  cellHit.channel = tempIter->channel;
701  cellHit.cell = tempIter->cell;
702  cellHit.module = tempIter->module;
703  cellHit.tray = tempIter->tray;
704  cellHit.hitPosition = tempIter->hitPosition;
705  cellHit.trackIdVec = trackIdVec;
706  cellHit.zhit = tempIter->zhit;
707  cellHit.yhit = tempIter->yhit;
708 
709  Float_t ycenter = (tempIter->cell-1-2.5)*mWidthPad;
710  Float_t dy = tempIter->yhit - ycenter;
711  Float_t dz = tempIter->zhit;
712 
713  if(mHisto) {
714  mTracksPerCellMatch1->Fill(trackIdVec.size());
715  mDaqOccupancyMatch1->Fill(tempIter->channel);
716  mDeltaHitMatch1->Fill(dy, dz);
717  }
718 
719  if (nTracks==1){
720  nSingleHitCells++;
721  singleHitCellsVec.push_back(cellHit);
722  } else if (nTracks>1){
723  nMultiHitsCells++;
724  multiHitsCellsVec.push_back(cellHit);
725  // for multiple hit cells either discard (yes) or
726  // find the most likely candidate.
727  } else {
728  LOG_INFO << "D: no tracks extrapolate to matched cell ... should not happen!" << endm;
729  }
730 
731  if (Debug()) {
732  LOG_INFO << "D: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
733  idVectorIter ij=trackIdVec.begin();
734  while (ij != trackIdVec.end()) { LOG_INFO << " " << *ij; ij++; }
735  LOG_INFO <<endm;
736  }
737 
738  tempVec = erasedVec;
739  }
740  if(Debug())
741  LOG_INFO << "D: before/after: " << matchHitCellsVec.size() << "/" << singleHitCellsVec.size() << endm;
742  //end of Sect.C
743  if(mHisto) {
744  mCellsPerEventMatch1->Fill(singleHitCellsVec.size()+multiHitsCellsVec.size());
745  if(singleHitCellsVec.size()) mEventCounterHisto->Fill(9);
746  }
747 
748  //.........................................................................
749  // E. sort and deal singleHitCellsVector for multiple cells associated to single tracks
750  //
751  tofCellHitVector FinalMatchedCellsVec;
752  // FinalMatchedCellsVec.clear();
753  tempVec = singleHitCellsVec;
754  if(mHisto) {
755  mCellsPerEventMatch2->Fill(tempVec.size());
756  for(unsigned int ii=0;ii<tempVec.size();ii++) {
757  mTracksPerCellMatch2->Fill(tempVec[ii].trackIdVec.size());
758  mDaqOccupancyMatch2->Fill(tempVec[ii].channel);
759  Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
760  Float_t dy = tempVec[ii].yhit-ycenter;
761  Float_t dz = tempVec[ii].zhit;
762  mDeltaHitMatch2->Fill(dy, dz);
763  }
764  }
765 
766  erasedVec = tempVec;
767  while (tempVec.size() != 0) {
768  StructCellHit cellHit;
769  Int_t nCells = 0;
770  idVector vTrackId;
771  vector<StThreeVectorD> vPosition;
772  vector<Int_t> vchannel, vtray, vmodule, vcell;
773  vector<Float_t> vzhit, vyhit;
774 
775  tofCellHitVectorIter tempIter=tempVec.begin();
776  tofCellHitVectorIter erasedIter=erasedVec.begin();
777  while(erasedIter!= erasedVec.end()) {
778  if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
779  nCells++;
780  vchannel.push_back(erasedIter->channel);
781  vtray.push_back(erasedIter->tray);
782  vmodule.push_back(erasedIter->module);
783  vcell.push_back(erasedIter->cell);
784  vPosition.push_back(erasedIter->hitPosition);
785  vTrackId.push_back(erasedIter->trackIdVec.back());
786  vzhit.push_back(erasedIter->zhit);
787  vyhit.push_back(erasedIter->yhit);
788 
789  erasedVec.erase(erasedIter);
790  erasedIter--;
791  }
792  erasedIter++;
793  }
794 
795 
796  if (nCells==1){
797  // for singly hit cell, copy data in singleHitCellsVec
798  cellHit.channel = vchannel[0];
799  cellHit.tray = vtray[0];
800  cellHit.module = vmodule[0];
801  cellHit.cell = vcell[0];
802  cellHit.trackIdVec.push_back(vTrackId[0]);
803  cellHit.hitPosition = vPosition[0];
804  cellHit.matchFlag = 0;
805  cellHit.zhit = vzhit[0];
806  cellHit.yhit = vyhit[0];
807 
808  FinalMatchedCellsVec.push_back(cellHit);
809 
810  // debugging output
811  if (Debug()) {
812  LOG_INFO << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
813  idVectorIter ij=vTrackId.begin();
814  while (ij != vTrackId.end()) { LOG_INFO << " " << *ij; ij++; }
815  LOG_INFO <<endm;
816  }
817  }
818  else if (nCells>1){ // for multiple hit cells find the most likely candidate.
819  Int_t thiscandidate(-99);
820  Int_t thisMatchFlag(0);
821 
822  // 1. sort on adc weight
823  Float_t weight(0);
824  vector<Int_t> weightCandidates;
825  thisMatchFlag = 1;
826  if (Debug()) LOG_INFO << "E: find ... weight ";
827  for (Int_t i=0;i<nCells;i++){
828  if(vchannel[i]>=0&&vchannel[i]<mNTOFR) {
829  Float_t adcwt = mTofrAdc[vchannel[i]]-mPedTOFr[vchannel[i]];
830  if(adcwt>weight) {
831  weight = adcwt;
832  weightCandidates.clear();
833  weightCandidates.push_back(i);
834  } else if (adcwt == weight) {
835  weightCandidates.push_back(i);
836  }
837  }
838  }
839  if (weightCandidates.size()==1){
840  thiscandidate = weightCandidates[0];
841  Int_t daqId = vchannel[thiscandidate];
842  if (Debug()) LOG_INFO << "candidate =" << daqId << endm;
843  }
844 
845  // 2. if still undecided check on hitposition
846  if (weightCandidates.size()>1){
847  Float_t ss(99.);
848  vector<Int_t> ssCandidates;
849  thisMatchFlag = 2;
850  if (Debug()) LOG_INFO << " ss ";
851  for (unsigned int i=0;i<weightCandidates.size();i++){
852  Int_t ii=weightCandidates[i];
853  // Int_t daqId = vchannel[ii];
854  StThreeVectorD ahit = vPosition[ii];
855  Float_t yy = ahit.y();
856  Float_t ycell = (vcell[ii]-1-2.5)*mWidthPad;
857  Float_t ll = fabs(yy-ycell);
858  if(ll<ss) {
859  ss = ll;
860  ssCandidates.clear();
861  ssCandidates.push_back(ii);
862  }else if (ll==ss)
863  ssCandidates.push_back(ii);
864  }
865  if (ssCandidates.size()==1){
866  thiscandidate = ssCandidates[0];
867  Int_t daqId = vchannel[thiscandidate];
868  if (Debug()) LOG_INFO << "candidate =" << daqId << endm;
869  }
870 
871  }
872 
873  if (thiscandidate>=0) {
874  cellHit.channel = vchannel[thiscandidate];
875  cellHit.tray = vtray[thiscandidate];
876  cellHit.module = vmodule[thiscandidate];
877  cellHit.cell = vcell[thiscandidate];
878  cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
879  cellHit.hitPosition = vPosition[thiscandidate];
880  cellHit.matchFlag = thisMatchFlag;
881  cellHit.zhit = vzhit[thiscandidate];
882  cellHit.yhit = vyhit[thiscandidate];
883 
884  FinalMatchedCellsVec.push_back(cellHit);
885 
886  // debugging output
887  if (Debug()) {
888  LOG_INFO << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:" << vTrackId[thiscandidate] << endm;
889  }
890  }
891 
892  } else {
893  LOG_INFO << "E: no cells belong to this track ... should not happen!" << endm;
894  }
895 
896  tempVec = erasedVec;
897  }
898 
899  LOG_INFO << "E: before/after: " << singleHitCellsVec.size() << "/" << FinalMatchedCellsVec.size() << endm;
900  // end of Sect.E
901 
902  //.........................................................................
903  // F. perform further selection and
904  // fill valid track histograms, ntuples and CellCollection
905  //
906  tempVec.clear();
907  tempVec = FinalMatchedCellsVec;
908  if(mHisto) {
909  if(FinalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
910  mCellsPerEventMatch3->Fill(tempVec.size());
911  for(unsigned int ii=0;ii<tempVec.size();ii++) {
912  mTracksPerCellMatch3->Fill(tempVec[ii].trackIdVec.size());
913  mDaqOccupancyMatch3->Fill(tempVec[ii].channel);
914  Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
915  Float_t dy = tempVec[ii].yhit - ycenter;
916  Float_t dz = tempVec[ii].zhit;
917  mDeltaHitMatch3->Fill(dy, dz);
918  }
919  }
920 
921  StTofCellCollection *mCellCollection = new StTofCellCollection;
922  Int_t nValidSingleHitCells(0), nValidSinglePrimHitCells(0);
923 
924  for (size_t ii=0; ii < FinalMatchedCellsVec.size(); ii++){
925  Int_t daqId = FinalMatchedCellsVec[ii].channel;
926  Int_t jj = daqId;
927  Int_t tray = FinalMatchedCellsVec[ii].tray;
928  Int_t module = FinalMatchedCellsVec[ii].module;
929  Int_t cell = FinalMatchedCellsVec[ii].cell;
930 
931  Float_t ycenter = (cell-1-2.5)*mWidthPad;
932  Float_t dy = FinalMatchedCellsVec[ii].yhit - ycenter;
933  if (FinalMatchedCellsVec[ii].trackIdVec.size()!=1)
934  LOG_INFO << "F: WHAT!?! mult.matched cell in single cell list " << daqId << endm;
935 
936  // 1. fill valid single track AND valid tdc histograms
937  if (validTdc(mTofrTdc[jj])) nValidSingleHitCells++;
938 
939  // get track-id from cell hit vector
940  unsigned int trackNode = FinalMatchedCellsVec[ii].trackIdVec[0];
941  StTrack *theTrack = nodes[trackNode]->track(primary);
942  StTrack *globalTrack = nodes[trackNode]->track(global);
943 
944  // 2. continue only if the (primary) track exists
945  if (validTofTrack(theTrack) && fabs(dy)<1.9 ){
946  nValidSinglePrimHitCells++;
947 
948  //--- store number of hits per track
949  Int_t nHitsPerTrack = theTrack->topologyMap().numberOfHits(kTpcId);
950 
951  // select the apropriate track geometry
952  StTrackGeometry *theTrackGeometry = trackGeometry(theTrack);
953 
954  //--- get momentum from track
955  const StThreeVectorF momentum = theTrackGeometry->momentum();
956 
957  //--- calculate flight path
958  Double_t pathLength = tofPathLength(&mEvent->primaryVertex()->position(),
959  &FinalMatchedCellsVec[ii].hitPosition,
960  theTrackGeometry->helix().curvature());
961 
962  //--- calculate local hit position on cell based first, last and middle plane
963  // (middle plane is the average of first and last plane, which is mathematically
964  // the same as CellHitVec.hitPosition ... )
965 // StThreeVectorD *pInnerLayer, *pOuterLayer;
966 // pInnerLayer = FinalMatchedCellsVec[ii].layerHitPositions.begin();
967 // pOuterLayer = FinalMatchedCellsVec[ii].layerHitPositions.end() - 1;
968 
969  //--- dig out from the dedx and rich pid traits
970  Float_t dedx(0.), cherang(0);
971  Int_t dedx_np(0), cherang_nph(0);
972  StSPtrVecTrackPidTraits& traits = theTrack->pidTraits();
973  for (unsigned int it=0;it<traits.size();it++){
974  if (traits[it]->detector() == kTpcId){
975  StDedxPidTraits *pid = dynamic_cast<StDedxPidTraits*>(traits[it]);
976  if (pid && pid->method() ==kTruncatedMeanId){
977  dedx = pid->mean();
978  dedx_np = pid->numberOfPoints();
979  }
980  } else if (traits[it]->detector() == kRichId){
981  StRichPidTraits *pid = dynamic_cast<StRichPidTraits*>(traits[it]);
982  if (pid){
983  StRichSpectra* pidinfo = pid->getRichSpectra();
984  if (pidinfo && pidinfo->getCherenkovPhotons()>2){
985  cherang = pidinfo->getCherenkovAngle();
986  cherang_nph = pidinfo->getCherenkovPhotons();
987  }
988  }
989  }
990  }
991 
992  //--- calculate local hit position on cell based on average hitposition
993  // Float_t localHitPos = mTofGeom->cellHitPosition(&allMatchedCellsVec[ii].hitPosition);
994 
995  // Fill TOF Cell Collection
996  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);
997  mCellCollection->push_back(tofCell);
998 
999  // dump debug data
1000  if (Debug()){
1001  LOG_INFO << "F: itray=" << tray << " imodule=" << module << " icell=" << cell << "\tnodeid:";
1002  idVectorIter ij=FinalMatchedCellsVec[ii].trackIdVec.begin();
1003  while (ij != FinalMatchedCellsVec[ii].trackIdVec.end()) { LOG_INFO << " " << *ij; ij++; }
1004  LOG_INFO << "\tR=" << 1/(theTrackGeometry->helix().curvature())
1005  << "\tpT=" << momentum.perp() << "\tp=" << momentum.mag()
1006  << "\thits="<< nHitsPerTrack << "\ts="<< pathLength
1007  << "\t#fitp=" <<theTrack->fitTraits().numberOfFitPoints(kTpcId)
1008  // << "\t#trkp=" <<theTrack->detectorInfo()->numberOfPoints(kTpcId)
1009  << " \tdedx=" << dedx
1010  << " \tdca="<< globalTrack->geometry()->helix().distance(mEvent->primaryVertex()->position())<<" and "<<theTrackGeometry->helix().distance(mEvent->primaryVertex()->position());
1011  if (cherang!=0) LOG_INFO << " \trich="<< cherang << " (" << cherang_nph << ")";
1012  LOG_INFO << endm;
1013  }
1014 
1015  } // track exists
1016  }
1017 
1018  storeMatchData(mCellCollection,theTof);
1019  delete mCellCollection;
1020 
1021  //check StEvent collections --
1022  if (theTof->dataPresent())
1023  LOG_INFO << " TofCollection: raw data container present" << endm;
1024  if (theTof->cellsPresent()){
1025  LOG_INFO << " TofCollection: cell container present."<<endm;
1026  if (Debug()){
1027  StSPtrVecTofCell& tmpCellTofVec = theTof->tofCells();
1028  for (size_t i = 0; i < tmpCellTofVec.size(); i++) {
1029  StTofCell* p = tmpCellTofVec[i];
1030  LOG_INFO << p->trayIndex() << " " << p->moduleIndex() << " " << p->cellIndex() << " " << p->adc() << " " << p->tdc() << " " << p->associatedTrack() << " " << p->matchFlag() << " " << p->position() << endm;
1031  }
1032  }
1033  }
1034  //-- end check
1035 
1036  LOG_INFO << "F: before/after" << FinalMatchedCellsVec.size() << "/" <<nValidSinglePrimHitCells << endm;
1037  // end of Sect.F
1038 
1039  LOG_INFO << "#(cell tracks): " << allCellsHitVec.size()
1040  << " #(hit cells): " << FinalMatchedCellsVec.size()
1041  << " #cells (valid tdc): " << nTdcTofr
1042  << "\n#(single hits): " << nSingleHitCells
1043  << " #(single valid hits): " << nValidSingleHitCells
1044  << " #(single prim valid hits): " << nValidSinglePrimHitCells
1045  << endm;
1046 
1047 
1048 
1049  if (doPrintMemoryInfo) {
1050  StMemoryInfo::instance()->snapshot();
1051  StMemoryInfo::instance()->print();
1052  }
1053  if (doPrintCpuInfo) {
1054  timer.stop();
1055  LOG_INFO << "CPU time for StTofrMatchMaker::Make(): "
1056  << timer.elapsedTime() << " sec\n" << endm;
1057  }
1058 
1059  LOG_INFO << "StTofrMatchMaker -- bye-bye" << endm;
1060  return kStOK;
1061 }
1062 
1063 //---------------------------------------------------------------------------
1064 Int_t StTofrMatchMaker::processEventYear5(){
1065  // leave as empty now
1066 
1067  if(mHisto) mEventCounterHisto->Fill(0);
1068  // event selection ...
1069  mEvent = (StEvent *) GetInputDS("StEvent");
1070  if (!validEvent(mEvent)){
1071  gMessMgr->Info("StTofrMatchMaker -- nothing to do ... bye-bye","OS");
1072  return kStOK;
1073  }
1074 
1075  // timing & memory info -only when requested-
1076  StTimer timer;
1077  if (doPrintCpuInfo) timer.start();
1078  if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
1079 
1080  //.........................................................................
1081  // check for tofCollection and fill local copy with ADC and TDC data
1082  StTofCollection *theTof = mEvent->tofCollection();
1083  // getTofData(theTof);
1084  // year5
1085  mSortTofRawData = new StSortTofRawData(theTof);
1086 
1087  IntVec validchannel = mSortTofRawData->GetValidChannel();
1088  // Any Hits in TOF+PVPD ?
1089  if(validchannel.size()<=0) {
1090  LOG_INFO << " No hits in TOF or pVPD! " << endm;
1091  // return kStOK;
1092  } else {
1093  LOG_INFO << " Number of TOF+pVPD fired hits: " << validchannel.size() << endm;
1094  }
1095  if(Debug()) {
1096  for(size_t iv=0;iv<validchannel.size();iv++) {
1097  //if(validchannel[iv]<0||validchannel[iv]>=mNTOFR5) continue; // skip pvpd
1098  LOG_INFO << " channel = " << validchannel[iv]<<endm;
1099  IntVec leTdc = mSortTofRawData->GetLeadingTdc(validchannel[iv]);
1100  IntVec teTdc = mSortTofRawData->GetTrailingTdc(validchannel[iv]);
1101  for(size_t iv1=0;iv1<leTdc.size();iv1++) {
1102  LOG_INFO << " leading Tdc = " << leTdc[iv1]<<endm;
1103  }
1104  for(size_t iv2=0;iv2<teTdc.size();iv2++) {
1105  LOG_INFO << " trailing Tdc = " << teTdc[iv2] << endm;
1106  }
1107  }
1108  } // end mdebug
1109 
1110  // number of primary tracks
1111  // (note: different meaning in event.root and richtof.root)
1112  Int_t refmult(0);
1113  refmult = uncorrectedNumberOfPrimaries(*mEvent);
1114 
1115  if (Debug()){
1116  LOG_INFO << " #Tracks :" << mEvent->summary()->numberOfTracks()
1117  << "\n #goodPrimaryTracks:" << mEvent->summary()->numberOfGoodPrimaryTracks()
1118  << "\n #uncorr.prim.tracks :" << refmult << endm;
1119  }
1120 
1121  //.........................................................................
1122  // A. build vector of candidate cells with valid ADC signals
1123  // idVector validCellIdVec;
1124  tofCellHitVector daqCellsHitVec;
1125  // daqCellsHitVec.clear();
1126  idVector validModuleVec;
1127 
1128  for(size_t ich=0;ich<validchannel.size();ich++) {
1129  int ichan = validchannel[ich];
1130  if(ichan<0||ichan>=mNTOFR5) continue; // tray
1131  IntVec map = mDaqMap->Tofr5TDCChan2Cell(ichan);
1132  if(map.size()!=3) continue;
1133 
1134  Int_t IDtray = map[0];
1135  Int_t IDmodule = map[1];
1136  Int_t IDcell = map[2];
1137 
1138  StructCellHit aDaqCellHit;
1139  aDaqCellHit.channel = ichan;
1140  aDaqCellHit.tray = IDtray;
1141  aDaqCellHit.module = IDmodule;
1142  aDaqCellHit.cell = IDcell;
1143  daqCellsHitVec.push_back(aDaqCellHit);
1144 
1145  // additional valid number configuration
1146  int id = IDtray*100+IDmodule;
1147  bool ifind = kFALSE;
1148  for(size_t iv=0;iv<validModuleVec.size();iv++) {
1149  if(id==validModuleVec[iv]) {
1150  ifind = kTRUE;
1151  break;
1152  }
1153  }
1154  if(!ifind) validModuleVec.push_back(id);
1155  }
1156 
1157  // end of Sect.A
1158  if(Debug()) {
1159  LOG_INFO << " total # of cells = " << daqCellsHitVec.size() << endm;
1160  for(size_t iv = 0;iv<validModuleVec.size();iv++) {
1161  LOG_INFO << " module # " << validModuleVec[iv] << " Valid! " << endm;
1162  }
1163  }
1164  if(mHisto) {
1165  mCellsMultInEvent->Fill(daqCellsHitVec.size());
1166  if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
1167  }
1168  if(!daqCellsHitVec.size()) return kStOK;
1169 
1170  //.........................................................................
1171  // B. loop over global tracks and determine all cell-track matches
1172  //
1173  tofCellHitVector allCellsHitVec;
1174  // allCellsHitVec.clear();
1175  StructCellHit cellHit;
1176 
1177  StSPtrVecTrackNode& nodes = mEvent->trackNodes();
1178  Int_t nAllTracks=0;
1179  for (unsigned int iNode=0; iNode<nodes.size(); iNode++){
1180  tofCellHitVector cellHitVec;
1181  // cellHitVec.clear();
1182  StTrack *theTrack = nodes[iNode]->track(global);
1183 
1184  // make sure we have a track, a miniDST might have removed it...
1185  if (validTrack(theTrack)){
1186  nAllTracks++;
1187  StPhysicalHelixD theHelix = trackGeometry(theTrack)->helix();
1188 
1189  IntVec projTrayVec;
1190  if(!mTofrGeom->projTrayVector(theHelix, projTrayVec)) continue;
1191 
1192  IntVec idVec;
1193  DoubleVec pathVec;
1194  PointVec crossVec;
1195 
1196 // idVec.clear();
1197 // pathVec.clear();
1198 // crossVec.clear();
1199 
1200  Int_t ncells = 0;
1201  // if(mTofrGeom->HelixCrossCellIds(theHelix,idVec,pathVec,crossVec) ) {
1202  if(mTofrGeom->HelixCrossCellIds(theHelix, validModuleVec, projTrayVec, idVec, pathVec, crossVec)) {
1203  Int_t cells = idVec.size();
1204  for (Int_t i=0; i<cells; i++) {
1205  Int_t icell,imodule,itray;
1206  Double_t local[3],global[3];
1207  for(Int_t i2=0;i2<3;i2++){
1208  local[i2]=0;
1209  }
1210  global[0]=crossVec[i].x();
1211  global[1]=crossVec[i].y();
1212  global[2]=crossVec[i].z();
1213  mTofrGeom->DecodeCellId(idVec[i], icell, imodule, itray);
1214  LOG_INFO << " decode " << idVec[i] << " to tray#" << itray << " module#" << imodule << " cell#" << icell << endm;
1215  StTofrGeomSensor* sensor =
1216  mTofrGeom->GetGeomSensor(imodule,itray);
1217  if(!sensor) {
1218  gMessMgr->Warning("","OS") << " No sensitive module in the projection??? -- Something weird!!! " << endm;
1219  continue;
1220  }
1221  sensor->Master2Local(&global[0],&local[0]);
1222  icell = sensor->FindCellIndex(local);
1223  // StThreeVectorD glo=sensor->GetCenterPosition();
1224  StThreeVectorD glo(global[0], global[1], global[2]);
1225  StThreeVectorD hitPos(local[0], local[1], local[2]);
1226  delete sensor;
1227  // if (local[2]<=3.4&&local[2]>=-2.7) {
1228  Int_t Iarray = mDaqMap->Tofr5Cell2TDCChan(itray, imodule, icell);
1229  if(Iarray>=mDAQOVERFLOW||Iarray<0) continue;
1230  ncells++;
1231  cellHit.channel = Iarray;
1232  cellHit.tray = itray;
1233  cellHit.module = imodule;
1234  cellHit.cell = icell;
1235  cellHit.trackIdVec.push_back(iNode);
1236  cellHit.hitPosition = glo; // global position
1237  cellHit.zhit = (Float_t)hitPos.z();
1238  cellHit.yhit = (Float_t)hitPos.y();
1239  cellHitVec.push_back(cellHit);
1240  allCellsHitVec.push_back(cellHit);
1241  if(mHisto) {
1242  mDaqOccupancyProj->Fill(Iarray);
1243  mHitsPosition->Fill(hitPos.y(), hitPos.z());
1244  }
1245 
1246  if(Debug()) {
1247  LOG_INFO <<"B: nodeid=" << iNode << " projected in " << " tray="<< itray << " module="<<imodule<<" cell="<<icell<<endm;
1248  LOG_INFO <<" hit position " << hitPos << endm;
1249  }
1250  // }
1251  } // for (Int_t i=0...)
1252  } // endif(helixcross...)
1253  if(ncells>0&&mHisto) mHitsMultPerTrack->Fill(ncells);
1254 
1255  } // if(ValidTrack)..
1256  } // loop over nodes
1257  if(Debug())
1258  LOG_INFO << "B: matched/available/total #tracknodes: " <<allCellsHitVec.size() << "/" <<nAllTracks << "/" << nodes.size() << endm;
1259  if(mHisto) {
1260  mHitsMultInEvent->Fill(allCellsHitVec.size());
1261  if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
1262  }
1263  // end of Sect.B
1264 
1265  //.........................................................................
1266  // C. Match find Neighbours -- identify crosstalk
1267  //
1268  tofCellHitVector matchHitCellsVec;
1269  // matchHitCellsVec.clear();
1270 
1271  tofCellHitVectorIter daqIter = daqCellsHitVec.begin();
1272  for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
1273  tofCellHitVectorIter proIter = allCellsHitVec.begin();
1274  for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
1275  if( (daqIter->tray==proIter->tray)&&
1276  (daqIter->module==proIter->module) &&
1277  ( ( (proIter->cell==6)&&((proIter->cell==daqIter->cell) ||
1278  (proIter->cell==daqIter->cell+1)) )
1279  || ( (proIter->cell==1)&&((proIter->cell==daqIter->cell) ||
1280  (proIter->cell==daqIter->cell-1)) )
1281  || ( (proIter->cell>=2&&proIter->cell<=6) &&
1282  ( (proIter->cell==daqIter->cell) ||
1283  (proIter->cell==daqIter->cell-1) ||
1284  (proIter->cell==daqIter->cell+1) ) ) ) ) {
1285  cellHit.channel = daqIter->channel;
1286  cellHit.tray = daqIter->tray;
1287  cellHit.module = daqIter->module;
1288  cellHit.cell = daqIter->cell;
1289  cellHit.hitPosition = proIter->hitPosition;
1290  cellHit.trackIdVec = proIter->trackIdVec;
1291  cellHit.zhit = proIter->zhit;
1292  cellHit.yhit = proIter->yhit;
1293  matchHitCellsVec.push_back(cellHit);
1294  }
1295  }
1296  } //end {sec. C}
1297  if(Debug()) {
1298  LOG_INFO << "C: before/after: " << allCellsHitVec.size() << "/" << matchHitCellsVec.size() << endm;
1299  }
1300  if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
1301 
1302  //.........................................................................
1303  // D. sort hit vectors and deal with (discard) cells matched by multiple tracks
1304  //
1305  Int_t nSingleHitCells(0);
1306  Int_t nMultiHitsCells(0);
1307 
1308  tofCellHitVector singleHitCellsVec;
1309  tofCellHitVector multiHitsCellsVec;
1310 // singleHitCellsVec.clear();
1311 // multiHitsCellsVec.clear();
1312 
1313  tofCellHitVector tempVec = matchHitCellsVec;
1314  tofCellHitVector erasedVec = tempVec;
1315  while (tempVec.size() != 0) {
1316  Int_t nTracks = 0;
1317  idVector trackIdVec;
1318 
1319  tofCellHitVectorIter tempIter=tempVec.begin();
1320  tofCellHitVectorIter erasedIter=erasedVec.begin();
1321  while(erasedIter!= erasedVec.end()) {
1322  if(tempIter->tray == erasedIter->tray &&
1323  tempIter->module == erasedIter->module &&
1324  tempIter->cell == erasedIter->cell) {
1325  nTracks++;
1326  trackIdVec.push_back(erasedIter->trackIdVec.back()); // merge
1327  erasedVec.erase(erasedIter);
1328  erasedIter--;
1329  }
1330  erasedIter++;
1331  }
1332 
1333  cellHit.channel = tempIter->channel;
1334  cellHit.cell = tempIter->cell;
1335  cellHit.module = tempIter->module;
1336  cellHit.tray = tempIter->tray;
1337  cellHit.hitPosition = tempIter->hitPosition;
1338  cellHit.trackIdVec = trackIdVec;
1339  cellHit.zhit = tempIter->zhit;
1340  cellHit.yhit = tempIter->yhit;
1341 
1342  Float_t ycenter = (tempIter->cell-1-2.5)*mWidthPad;
1343  Float_t dy = tempIter->yhit - ycenter;
1344  Float_t dz = tempIter->zhit;
1345 
1346  if(mHisto) {
1347  mTracksPerCellMatch1->Fill(trackIdVec.size());
1348  mDaqOccupancyMatch1->Fill(tempIter->channel);
1349  mDeltaHitMatch1->Fill(dy, dz);
1350  }
1351 
1352  if (nTracks==1){
1353  nSingleHitCells++;
1354  singleHitCellsVec.push_back(cellHit);
1355  } else if (nTracks>1){
1356  nMultiHitsCells++;
1357  multiHitsCellsVec.push_back(cellHit);
1358  // for multiple hit cells either discard (yes) or
1359  // find the most likely candidate.
1360  } else {
1361  LOG_INFO << "D: no tracks extrapolate to matched cell ... should not happen!" << endm;
1362  }
1363 
1364  if (Debug()) {
1365  LOG_INFO << "D: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
1366  idVectorIter ij=trackIdVec.begin();
1367  while (ij != trackIdVec.end()) { LOG_INFO << " " << *ij; ij++; }
1368  LOG_INFO <<endm;
1369  }
1370 
1371  tempVec = erasedVec;
1372  }
1373  if(Debug())
1374  LOG_INFO << "D: before/after: " << matchHitCellsVec.size() << "/" << singleHitCellsVec.size() << endm;
1375  //end of Sect.C
1376  if(mHisto) {
1377  mCellsPerEventMatch1->Fill(singleHitCellsVec.size()+multiHitsCellsVec.size());
1378  if(singleHitCellsVec.size()) mEventCounterHisto->Fill(9);
1379  }
1380 
1381  //.........................................................................
1382  // E. sort and deal singleHitCellsVector for multiple cells associated to single tracks
1383  //
1384  tofCellHitVector FinalMatchedCellsVec;
1385  // FinalMatchedCellsVec.clear();
1386  tempVec = singleHitCellsVec;
1387  if(mHisto) {
1388  mCellsPerEventMatch2->Fill(tempVec.size());
1389  for(unsigned int ii=0;ii<tempVec.size();ii++) {
1390  mTracksPerCellMatch2->Fill(tempVec[ii].trackIdVec.size());
1391  mDaqOccupancyMatch2->Fill(tempVec[ii].channel);
1392  Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
1393  Float_t dy = tempVec[ii].yhit-ycenter;
1394  Float_t dz = tempVec[ii].zhit;
1395  mDeltaHitMatch2->Fill(dy, dz);
1396  }
1397  }
1398 
1399  erasedVec = tempVec;
1400  while (tempVec.size() != 0) {
1401  StructCellHit cellHit;
1402  Int_t nCells = 0;
1403  idVector vTrackId;
1404  vector<StThreeVectorD> vPosition;
1405  vector<Int_t> vchannel, vtray, vmodule, vcell;
1406  vector<Float_t> vzhit, vyhit;
1407 
1408  tofCellHitVectorIter tempIter=tempVec.begin();
1409  tofCellHitVectorIter erasedIter=erasedVec.begin();
1410  while(erasedIter!= erasedVec.end()) {
1411  if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
1412  nCells++;
1413  vchannel.push_back(erasedIter->channel);
1414  vtray.push_back(erasedIter->tray);
1415  vmodule.push_back(erasedIter->module);
1416  vcell.push_back(erasedIter->cell);
1417  vPosition.push_back(erasedIter->hitPosition);
1418  vTrackId.push_back(erasedIter->trackIdVec.back());
1419  vzhit.push_back(erasedIter->zhit);
1420  vyhit.push_back(erasedIter->yhit);
1421 
1422  erasedVec.erase(erasedIter);
1423  erasedIter--;
1424  }
1425  erasedIter++;
1426  }
1427 
1428  if (nCells==1){
1429  // for singly hit cell, copy data in singleHitCellsVec
1430  cellHit.channel = vchannel[0];
1431  cellHit.tray = vtray[0];
1432  cellHit.module = vmodule[0];
1433  cellHit.cell = vcell[0];
1434  cellHit.trackIdVec.push_back(vTrackId[0]);
1435  cellHit.hitPosition = vPosition[0];
1436  cellHit.matchFlag = 0;
1437  cellHit.zhit = vzhit[0];
1438  cellHit.yhit = vyhit[0];
1439 
1440  FinalMatchedCellsVec.push_back(cellHit);
1441 
1442  // debugging output
1443  if (Debug()) {
1444  LOG_INFO << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
1445  idVectorIter ij=vTrackId.begin();
1446  while (ij != vTrackId.end()) { LOG_INFO << " " << *ij; ij++; }
1447  LOG_INFO <<endm;
1448  }
1449  }
1450  else if (nCells>1){ // for multiple hit cells find the most likely candidate.
1451  Int_t thiscandidate(-99);
1452  Int_t thisMatchFlag(0);
1453 
1454  // sort on hitposition
1455  Float_t ss(99.);
1456  vector<Int_t> ssCandidates;
1457  thisMatchFlag = 2;
1458  if (Debug()) LOG_INFO << " ss " << endm;
1459  for (Int_t i=0;i<nCells;i++){
1460  Float_t yy = vyhit[i];
1461  Float_t ycell = (vcell[i]-1-2.5)*mWidthPad;
1462  Float_t ll = fabs(yy-ycell);
1463  if(ll<ss) {
1464  ss = ll;
1465  ssCandidates.clear();
1466  ssCandidates.push_back(i);
1467  }else if (ll==ss)
1468  ssCandidates.push_back(i);
1469  }
1470  if (ssCandidates.size()==1){
1471  thiscandidate = ssCandidates[0];
1472  Int_t daqId = vchannel[thiscandidate];
1473  if (Debug()) LOG_INFO << "candidate =" << daqId << endm;
1474  }
1475 
1476 
1477  if (thiscandidate>=0) {
1478  cellHit.channel = vchannel[thiscandidate];
1479  cellHit.tray = vtray[thiscandidate];
1480  cellHit.module = vmodule[thiscandidate];
1481  cellHit.cell = vcell[thiscandidate];
1482  cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
1483  cellHit.hitPosition = vPosition[thiscandidate];
1484  cellHit.matchFlag = thisMatchFlag;
1485  cellHit.zhit = vzhit[thiscandidate];
1486  cellHit.yhit = vyhit[thiscandidate];
1487 
1488  FinalMatchedCellsVec.push_back(cellHit);
1489 
1490  // debugging output
1491  if (Debug()) {
1492  LOG_INFO << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:" << vTrackId[thiscandidate] << endm;
1493  }
1494  }
1495 
1496  } else {
1497  LOG_INFO << "E: no cells belong to this track ... should not happen!" << endm;
1498  }
1499 
1500  tempVec = erasedVec;
1501  }
1502 
1503  LOG_INFO << "E: before/after: " << singleHitCellsVec.size() << "/" << FinalMatchedCellsVec.size() << endm;
1504  // end of Sect.E
1505 
1506  //.........................................................................
1507  // F. perform further selection and
1508  // fill valid track histograms, ntuples and CellCollection
1509  //
1510  tempVec.clear();
1511  tempVec = FinalMatchedCellsVec;
1512  if(mHisto) {
1513  if(FinalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
1514  mCellsPerEventMatch3->Fill(tempVec.size());
1515  for(unsigned int ii=0;ii<tempVec.size();ii++) {
1516  mTracksPerCellMatch3->Fill(tempVec[ii].trackIdVec.size());
1517  mDaqOccupancyMatch3->Fill(tempVec[ii].channel);
1518  Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
1519  Float_t dy = tempVec[ii].yhit - ycenter;
1520  Float_t dz = tempVec[ii].zhit;
1521  mDeltaHitMatch3->Fill(dy, dz);
1522  }
1523  }
1524 
1525  StTofCellCollection *mCellCollection = new StTofCellCollection;
1526  Int_t nValidSingleHitCells(0), nValidSinglePrimHitCells(0);
1527 
1528  for (size_t ii=0; ii < FinalMatchedCellsVec.size(); ii++){
1529  Int_t daqId = FinalMatchedCellsVec[ii].channel;
1530  Int_t jj = daqId;
1531  Int_t tray = FinalMatchedCellsVec[ii].tray;
1532  Int_t module = FinalMatchedCellsVec[ii].module;
1533  Int_t cell = FinalMatchedCellsVec[ii].cell;
1534 
1535  Float_t ycenter = (cell-1-2.5)*mWidthPad;
1536  Float_t dy = FinalMatchedCellsVec[ii].yhit - ycenter;
1537  if (FinalMatchedCellsVec[ii].trackIdVec.size()!=1)
1538  LOG_INFO << "F: WHAT!?! mult.matched cell in single cell list " << daqId << endm;
1539 
1540 
1541  /* move INL correction to calibration maker
1542  // Read in Leading and Trailing edge TDC, apply on INL correction
1543  int tmptdc = (mSortTofRawData->GetLeadingTdc(jj))[0];
1544  // do inl correction
1545  int bin = int(tmptdc)&0x03ff;
1546  float tmptdc_f = tmptdc + GetINLcorr(4, jj, bin);
1547  float letime = tmptdc_f*25./1024; // ns
1548 
1549  tmptdc=(mSortTofRawData->GetTrailingTdc(jj))[0];
1550  // do inl correction
1551  bin = int(tmptdc)&0x0ff;
1552  tmptdc_f = tmptdc + GetINLcorr(5, jj, bin);
1553  float tetime = tmptdc_f*100./1024; // ns
1554  float tot = tetime - letime; // ns
1555  */
1556 
1557  int letdc = (mSortTofRawData->GetLeadingTdc(jj))[0];
1558  int tetdc=(mSortTofRawData->GetTrailingTdc(jj))[0];
1559 
1560  // get track-id from cell hit vector
1561  unsigned int trackNode = FinalMatchedCellsVec[ii].trackIdVec[0];
1562  StTrack *theTrack = nodes[trackNode]->track(primary);
1563  StTrack *globalTrack = nodes[trackNode]->track(global);
1564 
1565  // 2. continue only if the (primary) track exists
1566  if (validTofTrack(theTrack) && fabs(dy)<1.9 ){
1567  nValidSinglePrimHitCells++;
1568 
1569  //--- store number of hits per track
1570  Int_t nHitsPerTrack = theTrack->topologyMap().numberOfHits(kTpcId);
1571 
1572  // select the apropriate track geometry
1573  StTrackGeometry *theTrackGeometry = trackGeometry(theTrack);
1574 
1575  //--- get momentum from track
1576  const StThreeVectorF momentum = theTrackGeometry->momentum();
1577 
1578  //--- calculate flight path
1579  Double_t pathLength = tofPathLength(&mEvent->primaryVertex()->position(),
1580  &FinalMatchedCellsVec[ii].hitPosition,
1581  theTrackGeometry->helix().curvature());
1582 
1583  //--- calculate local hit position on cell based first, last and middle plane
1584  // (middle plane is the average of first and last plane, which is mathematically
1585  // the same as CellHitVec.hitPosition ... )
1586 // StThreeVectorD *pInnerLayer, *pOuterLayer;
1587 // pInnerLayer = FinalMatchedCellsVec[ii].layerHitPositions.begin();
1588 // pOuterLayer = FinalMatchedCellsVec[ii].layerHitPositions.end() - 1;
1589 
1590  //--- dig out from the dedx and rich pid traits
1591  Float_t dedx(0.), cherang(0);
1592  Int_t dedx_np(0), cherang_nph(0);
1593  StSPtrVecTrackPidTraits& traits = theTrack->pidTraits();
1594  for (unsigned int it=0;it<traits.size();it++){
1595  if (traits[it]->detector() == kTpcId){
1596  StDedxPidTraits *pid = dynamic_cast<StDedxPidTraits*>(traits[it]);
1597  if (pid && pid->method() ==kTruncatedMeanId){
1598  dedx = pid->mean();
1599  dedx_np = pid->numberOfPoints();
1600  }
1601  } else if (traits[it]->detector() == kRichId){
1602  StRichPidTraits *pid = dynamic_cast<StRichPidTraits*>(traits[it]);
1603  if (pid){
1604  StRichSpectra* pidinfo = pid->getRichSpectra();
1605  if (pidinfo && pidinfo->getCherenkovPhotons()>2){
1606  cherang = pidinfo->getCherenkovAngle();
1607  cherang_nph = pidinfo->getCherenkovPhotons();
1608  }
1609  }
1610  }
1611  }
1612 
1613  //--- calculate local hit position on cell based on average hitposition
1614  // Float_t localHitPos = mTofGeom->cellHitPosition(&allMatchedCellsVec[ii].hitPosition);
1615 
1616  // Fill TOF Cell Collection
1617  // year5, raw letdc and tetdc stored as "tdc" and "adc" in StTofCell
1618  //
1619  // 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);
1620  StTofCell *tofCell = new StTofCell(tray, module, cell, daqId, tetdc,letdc,theTrack,FinalMatchedCellsVec[ii].zhit,FinalMatchedCellsVec[ii].matchFlag,FinalMatchedCellsVec[ii].hitPosition);
1621  mCellCollection->push_back(tofCell);
1622 
1623  // dump debug data
1624  if (Debug()){
1625  LOG_INFO << "F: itray=" << tray << " imodule=" << module << " icell=" << cell << "\tnodeid:";
1626  idVectorIter ij=FinalMatchedCellsVec[ii].trackIdVec.begin();
1627  while (ij != FinalMatchedCellsVec[ii].trackIdVec.end()) { LOG_INFO << " " << *ij; ij++; }
1628  LOG_INFO << "\tR=" << 1/(theTrackGeometry->helix().curvature())
1629  << "\tpT=" << momentum.perp() << "\tp=" << momentum.mag()
1630  << "\thits="<< nHitsPerTrack << "\ts="<< pathLength
1631  << "\t#fitp=" <<theTrack->fitTraits().numberOfFitPoints(kTpcId)
1632  // << "\t#trkp=" <<theTrack->detectorInfo()->numberOfPoints(kTpcId)
1633  << " \tdedx=" << dedx
1634  << " \tdca="<< globalTrack->geometry()->helix().distance(mEvent->primaryVertex()->position())<<" and "<<theTrackGeometry->helix().distance(mEvent->primaryVertex()->position());
1635  if (cherang!=0) LOG_INFO << " \trich="<< cherang << " (" << cherang_nph << ")";
1636  LOG_INFO << endm;
1637  }
1638 
1639  } // track exists
1640  } // end final matched cells
1641 
1642  storeMatchData(mCellCollection,theTof);
1643  delete mCellCollection;
1644 
1645  delete mSortTofRawData;
1646  mSortTofRawData = 0;
1647 
1648  //check StEvent collections --
1649  if (theTof->dataPresent())
1650  LOG_INFO << " TofCollection: raw data container present" << endm;
1651  if (theTof->cellsPresent()){
1652  LOG_INFO << " TofCollection: cell container present."<<endm;
1653  if (Debug()){
1654  StSPtrVecTofCell& tmpCellTofVec = theTof->tofCells();
1655  LOG_INFO << " # of matched cells " << tmpCellTofVec.size() << endm;
1656  for (size_t i = 0; i < tmpCellTofVec.size(); i++) {
1657  StTofCell* p = tmpCellTofVec[i];
1658  LOG_INFO << p->trayIndex() << " " << p->moduleIndex() << " " << p->cellIndex() << " " << p->trailingEdgeTime() << " " << p->leadingEdgeTime() << " " << p->associatedTrack() << " " << p->matchFlag() << " " << p->position() << endm;
1659  }
1660  }
1661  }
1662  //-- end check
1663 
1664  LOG_INFO << "F: before/after" << FinalMatchedCellsVec.size() << "/" <<nValidSinglePrimHitCells << endm;
1665  // end of Sect.F
1666 
1667  LOG_INFO << "#(cell tracks): " << allCellsHitVec.size()
1668  << " #(hit cells): " << FinalMatchedCellsVec.size()
1669  << "\n#(single hits): " << nSingleHitCells
1670  << " #(single valid hits): " << nValidSingleHitCells
1671  << " #(single prim valid hits): " << nValidSinglePrimHitCells
1672  << endm;
1673 
1674 
1675 
1676  if (doPrintMemoryInfo) {
1677  StMemoryInfo::instance()->snapshot();
1678  StMemoryInfo::instance()->print();
1679  }
1680  if (doPrintCpuInfo) {
1681  timer.stop();
1682  LOG_INFO << "CPU time for StTofrMatchMaker::Make(): "
1683  << timer.elapsedTime() << " sec\n" << endm;
1684  }
1685 
1686  LOG_INFO << "StTofrMatchMaker -- bye-bye" << endm;
1687 
1688 
1689 
1690  return kStOK;
1691 }
1692 
1693 //---------------------------------------------------------------------------
1694 Int_t StTofrMatchMaker::processEventYear8(){
1695  // leave as empty now
1696 
1697  if(Debug()) LOG_INFO << " processing event in run 8 " << endm;
1698  if(mHisto) mEventCounterHisto->Fill(0);
1699  // event selection ...
1700  mEvent = (StEvent *) GetInputDS("StEvent");
1701 // if (!validEvent(mEvent)){
1702  if(!mEvent || !(mEvent->tofCollection()) || !(mEvent->tofCollection()->rawdataPresent()) ) {
1703  gMessMgr->Info("StTofrMatchMaker -- nothing to do ... bye-bye","OS");
1704  return kStOK;
1705  }
1706  if(mHisto) mEventCounterHisto->Fill(1);
1707 
1708  // timing & memory info -only when requested-
1709  StTimer timer;
1710  if (doPrintCpuInfo) timer.start();
1711  if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
1712 
1713  //.........................................................................
1714  // check for tofCollection and fill local copy with ADC and TDC data
1715  StTofCollection *theTof = mEvent->tofCollection();
1716 
1717  //.........................................................................
1718  // push raw data into StTofData for QA:
1719  // Be careful, the vpd trayId in TofData is 121 and 122 (**)
1720  // and
1721  //.........................................................................
1722  // A. build vector of candidate cells with valid ADC signals
1723  // idVector validCellIdVec;
1724  tofCellHitVector daqCellsHitVec;
1725  // daqCellsHitVec.clear();
1726  idVector validModuleVec;
1727 
1728  mSortTofRawData = new StSortTofRawData(theTof, mDaqMap);
1729 
1730  // multi-tray system
1731  IntVec validtray = mDaqMap->ValidTrays();
1732  for(size_t i=0;i<validtray.size();i++) {
1733  int trayId = validtray[i];
1734  IntVec validchannel = mSortTofRawData->GetValidChannel(trayId);
1735  if(Debug()) LOG_INFO << " Number of fired hits on tray " << trayId << " = " << validchannel.size() << endm;
1736 
1737  for(size_t iv=0;iv<validchannel.size();iv++) {
1738  IntVec leTdc = mSortTofRawData->GetLeadingTdc(trayId, validchannel[iv], kTRUE);
1739  IntVec teTdc = mSortTofRawData->GetTrailingTdc(trayId, validchannel[iv], kTRUE);
1740 
1741  if(!leTdc.size() || !teTdc.size()) continue;
1742 
1743  int chan = validchannel[iv];
1744  IntVec map = mDaqMap->TDIGChan2Cell(chan);
1745  int moduleId = map[0];
1746  int cellId = map[1];
1747 
1748  StructCellHit aDaqCellHit;
1749  aDaqCellHit.channel = chan;
1750  aDaqCellHit.tray = trayId;
1751  aDaqCellHit.module = moduleId;
1752  aDaqCellHit.cell = cellId;
1753  daqCellsHitVec.push_back(aDaqCellHit);
1754 
1755  // additional valid number configuration
1756  int id = trayId*100+moduleId;
1757  bool ifind = kFALSE;
1758  for(size_t im=0;im<validModuleVec.size();im++) {
1759  if(id==validModuleVec[im]) {
1760  ifind = kTRUE;
1761  break;
1762  }
1763  }
1764  if(!ifind) validModuleVec.push_back(id);
1765 
1766  if(mHisto) {
1767  // mDaqOccupancyValid->Fill(chan);
1768  mDaqOccupancyValid->Fill((moduleId-1)*mNCell+(cellId-1));
1769  mDaqOccupancyValidAll->Fill((trayId-76)*mNTOF+(moduleId-1)*mNCell+(cellId-1));
1770  }
1771  //
1772  // store data from trays and vpds into StTofData
1773  //
1774  int dataIndex = (trayId-1)*mNTOF + (moduleId-1)*mNCell + (cellId-1);
1775  StTofData *aData = new StTofData(dataIndex,0,0,0,0,leTdc[0],teTdc[0]);
1776  theTof->addData(aData);
1777 
1778  if(Debug()) {
1779  for(size_t iv1=0;iv1<leTdc.size();iv1++) {
1780  LOG_INFO << " leading Tdc = " << leTdc[iv1]<<endm;
1781  }
1782  for(size_t iv2=0;iv2<teTdc.size();iv2++) {
1783  LOG_INFO << " trailing Tdc = " << teTdc[iv2] << endm;
1784  }
1785  } // end debug
1786  } // end channel
1787 
1788  } // end tray
1789 
1790  // vpd -> StTofData
1791  for(int ivpd=0;ivpd<2;ivpd++) { // west and east sides
1792  int trayId = (ivpd==0) ? mWestVpdTrayId : mEastVpdTrayId;
1793  IntVec validtube = mSortTofRawData->GetValidChannel(trayId);
1794  if(Debug()) LOG_INFO << " Number of fired hits on tray(vpd) " << trayId << " = " << validtube.size() << endm;
1795 
1796  if(!validtube.size()) continue;
1797  for(int i=0;i<mNVPD;i++) {
1798  int tubeId = i+1;
1799  int lechan = (ivpd==0) ? mDaqMap->WestPMT2TDIGLeChan(tubeId) : mDaqMap->EastPMT2TDIGLeChan(tubeId);
1800  IntVec leTdc = mSortTofRawData->GetLeadingTdc(trayId, lechan, kTRUE);
1801  IntVec teTdc = mSortTofRawData->GetTrailingTdc(trayId, lechan, kTRUE); // channel number should be le, sorted in StSortTofRawData
1802 
1803  if(leTdc.size()&&mHisto) mDaqOccupancyVpd->Fill(ivpd*mNVPD+i);
1804 
1805  if(leTdc.size() && teTdc.size()) {
1806  int dataIndex = (ivpd+120)*mNTOF + (tubeId-1);
1807  StTofData *aData = new StTofData(dataIndex,0,0,0,0,leTdc[0],teTdc[0]);
1808  theTof->addData(aData);
1809 
1810  if(mHisto) mDaqOccupancyValidVpd->Fill(ivpd*mNVPD+i);
1811  }
1812  }
1813  }
1814  //
1815 
1816  // end of Sect.A
1817  if(Debug()) {
1818  LOG_INFO << " total # of cells = " << daqCellsHitVec.size() << endm;
1819  for(size_t iv = 0;iv<validModuleVec.size();iv++) {
1820  LOG_INFO << " module # " << validModuleVec[iv] << " Valid! " << endm;
1821  }
1822  }
1823  if(mHisto) {
1824  mCellsMultInEvent->Fill(daqCellsHitVec.size());
1825  if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
1826  }
1827 // if(!daqCellsHitVec.size()) return kStOK;
1828 
1829  //.........................................................................
1830  // B. loop over global tracks and determine all cell-track matches
1831  //
1832  tofCellHitVector allCellsHitVec;
1833  // allCellsHitVec.clear();
1834  StructCellHit cellHit;
1835 
1836  StSPtrVecTrackNode& nodes = mEvent->trackNodes();
1837  Int_t nAllTracks=0;
1838  for (unsigned int iNode=0; iNode<nodes.size(); iNode++){
1839  tofCellHitVector cellHitVec;
1840  // cellHitVec.clear();
1841  StGlobalTrack *theTrack = dynamic_cast<StGlobalTrack*>(nodes[iNode]->track(global));
1842  if(!theTrack) continue;
1843 
1844  StThreeVectorF mom = theTrack->geometry()->momentum();
1845  float pt = mom.perp();
1846  float eta = mom.pseudoRapidity();
1847  float phi = mom.phi();
1848  if (phi<0.) phi += 2.*3.14159;
1849 
1850  float dEdx = -999.;
1851  int ndEdxpts = 0;
1852  float nSigmaPion = -999.;
1853  static StTpcDedxPidAlgorithm PidAlgorithm;
1854  static StPionPlus* Pion = StPionPlus::instance();
1855 
1856  const StParticleDefinition* pd = theTrack->pidTraits(PidAlgorithm);
1857  if (pd) {
1858  nSigmaPion = PidAlgorithm.numberOfSigma(Pion);
1859  }
1860  if ( PidAlgorithm.traits() ) {
1861  dEdx = PidAlgorithm.traits()->mean();
1862  ndEdxpts = PidAlgorithm.traits()->numberOfPoints();
1863  }
1864  int nfitpts = theTrack->fitTraits().numberOfFitPoints(kTpcId);
1865 // cout << " dEdx = " << dEdx << endl;
1866 // cout << " nSigmaPi = " << nSigmaPion << endl;
1867 // if (pt<0.2) continue;
1868 
1869  // make sure we have a track, a miniDST might have removed it...
1870  if (validTrackRun8(theTrack)){
1871  if(mHisto) {
1872  mTrackPtEta->Fill(pt, eta);
1873  mTrackPtPhi->Fill(pt, phi);
1874  mTrackNFitPts->Fill(nfitpts);
1875  if(dEdx>0.) mTrackdEdxvsp->Fill(mom.mag(), dEdx*1.e6);
1876  if(fabs(nSigmaPion)<5.) mNSigmaPivsPt->Fill(pt, nSigmaPion+5.*theTrack->geometry()->charge());
1877  }
1878 
1879  if(mSaveTree) {
1880  trackTree.pt = pt;
1881  trackTree.eta = eta;
1882  trackTree.phi = phi;
1883  trackTree.nfitpts = nfitpts;
1884  trackTree.dEdx = dEdx*1.e6;
1885  trackTree.ndEdxpts = ndEdxpts;
1886  trackTree.charge = theTrack->geometry()->charge();
1887  trackTree.projTrayId = 0;
1888  trackTree.projCellChan = -1;
1889  trackTree.projY = -999.;
1890  trackTree.projZ = -999.;
1891  }
1892 
1893  nAllTracks++;
1894  StPhysicalHelixD theHelix = trackGeometry(theTrack)->helix();
1895 
1896 // IntVec projTrayVec;
1897 // if(!mTofrGeom->projTrayVector(theHelix, projTrayVec)) continue;
1898 
1899  IntVec idVec;
1900  DoubleVec pathVec;
1901  PointVec crossVec;
1902 
1903 // idVec.clear();
1904 // pathVec.clear();
1905 // crossVec.clear();
1906 
1907  Int_t ncells = 0;
1908  if(mTofrGeom->HelixCrossCellIds(theHelix,idVec,pathVec,crossVec) ) {
1909 // if(mTofrGeom->HelixCrossCellIds(theHelix, validModuleVec, projTrayVec, idVec, pathVec, crossVec)) {
1910  Int_t cells = idVec.size();
1911  for (Int_t i=0; i<cells; i++) {
1912  Int_t icell,imodule,itray;
1913  Double_t local[3],global[3];
1914  for(Int_t i2=0;i2<3;i2++){
1915  local[i2]=0;
1916  }
1917  global[0]=crossVec[i].x();
1918  global[1]=crossVec[i].y();
1919  global[2]=crossVec[i].z();
1920  mTofrGeom->DecodeCellId(idVec[i], icell, imodule, itray);
1921  LOG_INFO << " decode " << idVec[i] << " to tray#" << itray << " module#" << imodule << " cell#" << icell << endm;
1922  StTofrGeomSensor* sensor =
1923  mTofrGeom->GetGeomSensor(imodule,itray);
1924  if(!sensor) {
1925  gMessMgr->Warning("","OS") << " No sensitive module in the projection??? -- Something weird!!! " << endm;
1926  continue;
1927  }
1928  sensor->Master2Local(&global[0],&local[0]);
1929  icell = sensor->FindCellIndex(local);
1930  // StThreeVectorD glo=sensor->GetCenterPosition();
1931  StThreeVectorD glo(global[0], global[1], global[2]);
1932  StThreeVectorD hitPos(local[0], local[1], local[2]);
1933  delete sensor;
1934  if (local[2]<=2.7&&local[2]>=-3.4) {
1935  Int_t Iarray = mDaqMap->Cell2TDIGChan(imodule, icell);
1936  if(Iarray>=mDAQOVERFLOW||Iarray<0) continue;
1937  ncells++;
1938  cellHit.channel = Iarray;
1939  cellHit.tray = itray;
1940  cellHit.module = imodule;
1941  cellHit.cell = icell;
1942  cellHit.trackIdVec.push_back(iNode);
1943  cellHit.hitPosition = glo; // global position
1944  cellHit.zhit = (Float_t)hitPos.z();
1945  cellHit.yhit = (Float_t)hitPos.y();
1946  cellHitVec.push_back(cellHit);
1947  allCellsHitVec.push_back(cellHit);
1948  if(mHisto) {
1949 // mDaqOccupancyProj->Fill(Iarray);
1950  if(itray>=76&&itray<=80) {
1951  mDaqOccupancyProj->Fill((imodule-1)*mNCell+(icell-1));
1952  mDaqOccupancyProjAll->Fill((itray-76)*mNTOF+(imodule-1)*mNCell+(icell-1));
1953  mHitsPosition->Fill(hitPos.y(), hitPos.z());
1954  }
1955  }
1956 
1957  if(mSaveTree) {
1958  trackTree.projTrayId = itray;
1959  trackTree.projCellChan = (imodule-1)*mNCell+(icell-1);
1960  trackTree.projY = local[1];
1961  trackTree.projZ = local[2];
1962  }
1963 
1964  if(Debug()) {
1965  LOG_INFO <<"B: nodeid=" << iNode << " projected in " << " tray="<< itray << " module="<<imodule<<" cell="<<icell<<endm;
1966  LOG_INFO <<" hit position " << hitPos << endm;
1967  }
1968  }
1969  } // for (Int_t i=0...)
1970  } // endif(helixcross...)
1971  if(ncells>0&&mHisto) mHitsMultPerTrack->Fill(ncells);
1972 
1973  if(mHisto && mSaveTree) mTrackTree->Fill();
1974 
1975  } // if(ValidTrack)..
1976  } // loop over nodes
1977  if(Debug())
1978  LOG_INFO << "B: matched/available/total #tracknodes: " <<allCellsHitVec.size() << "/" <<nAllTracks << "/" << nodes.size() << endm;
1979  if(mHisto) {
1980  mHitsMultInEvent->Fill(allCellsHitVec.size());
1981  if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
1982  }
1983  // end of Sect.B
1984 
1985  //.........................................................................
1986  // C. Match find Neighbours -- identify crosstalk
1987  //
1988  tofCellHitVector matchHitCellsVec;
1989  // matchHitCellsVec.clear();
1990 
1991  tofCellHitVectorIter daqIter = daqCellsHitVec.begin();
1992  for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
1993  tofCellHitVectorIter proIter = allCellsHitVec.begin();
1994  for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
1995 
1996  int daqIndex = (daqIter->module-1)*6 + (daqIter->cell-1);
1997  int proIndex = (proIter->module-1)*6 + (proIter->cell-1);
1998  int hisIndex = daqIter->tray - 76;
1999  int daqAllIndex = (daqIter->tray - 76)*192 + daqIndex;
2000  int proAllIndex = (proIter->tray - 76)*192 + proIndex;
2001  if (mHisto) mHitCorrAll->Fill(proAllIndex,daqAllIndex);
2002  if(daqIter->tray==proIter->tray) {
2003  if (mHisto) {
2004  if(hisIndex>=0&&hisIndex<5) {
2005  mHitCorr[hisIndex]->Fill(proIndex,daqIndex);
2006  mHitCorrModule[hisIndex]->Fill(proIter->module-1,daqIter->module-1);
2007  } else {
2008  cout << " weird tray # " << daqIter->tray << endl;
2009  }
2010  }
2011  }
2012  if( (daqIter->tray==proIter->tray)&&
2013  (daqIter->module==proIter->module) &&
2014  ( ( (proIter->cell==6)&&((proIter->cell==daqIter->cell) ||
2015  (proIter->cell==daqIter->cell+1)) )
2016  || ( (proIter->cell==1)&&((proIter->cell==daqIter->cell) ||
2017  (proIter->cell==daqIter->cell-1)) )
2018  || ( (proIter->cell>=2&&proIter->cell<=6) &&
2019  ( (proIter->cell==daqIter->cell) ||
2020  (proIter->cell==daqIter->cell-1) ||
2021  (proIter->cell==daqIter->cell+1) ) ) ) ) {
2022  cellHit.channel = daqIter->channel;
2023  cellHit.tray = daqIter->tray;
2024  cellHit.module = daqIter->module;
2025  cellHit.cell = daqIter->cell;
2026  cellHit.hitPosition = proIter->hitPosition;
2027  cellHit.trackIdVec = proIter->trackIdVec;
2028  cellHit.zhit = proIter->zhit;
2029  cellHit.yhit = proIter->yhit;
2030  matchHitCellsVec.push_back(cellHit);
2031  }
2032  }
2033  } //end {sec. C}
2034  if(Debug()) {
2035  LOG_INFO << "C: before/after: " << allCellsHitVec.size() << "/" << matchHitCellsVec.size() << endm;
2036  }
2037  if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
2038 
2039  //.........................................................................
2040  // D. sort hit vectors and deal with (discard) cells matched by multiple tracks
2041  //
2042  Int_t nSingleHitCells(0);
2043  Int_t nMultiHitsCells(0);
2044 
2045  tofCellHitVector singleHitCellsVec;
2046  tofCellHitVector multiHitsCellsVec;
2047 // singleHitCellsVec.clear();
2048 // multiHitsCellsVec.clear();
2049 
2050  tofCellHitVector tempVec = matchHitCellsVec;
2051  tofCellHitVector erasedVec = tempVec;
2052  while (tempVec.size() != 0) {
2053  Int_t nTracks = 0;
2054  idVector trackIdVec;
2055 
2056  tofCellHitVectorIter tempIter=tempVec.begin();
2057  tofCellHitVectorIter erasedIter=erasedVec.begin();
2058  while(erasedIter!= erasedVec.end()) {
2059  if(tempIter->tray == erasedIter->tray &&
2060  tempIter->module == erasedIter->module &&
2061  tempIter->cell == erasedIter->cell) {
2062  nTracks++;
2063  trackIdVec.push_back(erasedIter->trackIdVec.back()); // merge
2064  erasedVec.erase(erasedIter);
2065  erasedIter--;
2066  }
2067  erasedIter++;
2068  }
2069 
2070  cellHit.channel = tempIter->channel;
2071  cellHit.cell = tempIter->cell;
2072  cellHit.module = tempIter->module;
2073  cellHit.tray = tempIter->tray;
2074  cellHit.hitPosition = tempIter->hitPosition;
2075  cellHit.trackIdVec = trackIdVec;
2076  cellHit.zhit = tempIter->zhit;
2077  cellHit.yhit = tempIter->yhit;
2078 
2079  Float_t ycenter = (tempIter->cell-1-2.5)*mWidthPad;
2080  Float_t dy = tempIter->yhit - ycenter;
2081  Float_t dz = tempIter->zhit;
2082 
2083  if(mHisto) {
2084  mTracksPerCellMatch1->Fill(trackIdVec.size());
2085 // mDaqOccupancyMatch1->Fill(tempIter->channel);
2086  mDaqOccupancyMatch1->Fill((tempIter->module-1)*mNCell+(tempIter->cell-1));
2087  mDeltaHitMatch1->Fill(dy, dz);
2088  }
2089 
2090  if (nTracks==1){
2091  nSingleHitCells++;
2092  singleHitCellsVec.push_back(cellHit);
2093  } else if (nTracks>1){
2094  nMultiHitsCells++;
2095  multiHitsCellsVec.push_back(cellHit);
2096  // for multiple hit cells either discard (yes) or
2097  // find the most likely candidate.
2098  } else {
2099  LOG_INFO << "D: no tracks extrapolate to matched cell ... should not happen!" << endm;
2100  }
2101 
2102  if (Debug()) {
2103  LOG_INFO << "D: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
2104  idVectorIter ij=trackIdVec.begin();
2105  while (ij != trackIdVec.end()) { LOG_INFO << " " << *ij; ij++; }
2106  LOG_INFO <<endm;
2107  }
2108 
2109  tempVec = erasedVec;
2110  }
2111  if(Debug())
2112  LOG_INFO << "D: before/after: " << matchHitCellsVec.size() << "/" << singleHitCellsVec.size() << endm;
2113  //end of Sect.C
2114  if(mHisto) {
2115  mCellsPerEventMatch1->Fill(singleHitCellsVec.size()+multiHitsCellsVec.size());
2116  if(singleHitCellsVec.size()) mEventCounterHisto->Fill(9);
2117  }
2118 
2119  //.........................................................................
2120  // E. sort and deal singleHitCellsVector for multiple cells associated to single tracks
2121  //
2122  tofCellHitVector FinalMatchedCellsVec;
2123  // FinalMatchedCellsVec.clear();
2124  tempVec = singleHitCellsVec;
2125  if(mHisto) {
2126  mCellsPerEventMatch2->Fill(tempVec.size());
2127  for(unsigned int ii=0;ii<tempVec.size();ii++) {
2128  mTracksPerCellMatch2->Fill(tempVec[ii].trackIdVec.size());
2129  // mDaqOccupancyMatch2->Fill(tempVec[ii].channel);
2130  mDaqOccupancyMatch2->Fill((tempVec[ii].module-1)*mNCell+(tempVec[ii].cell-1));
2131  Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
2132  Float_t dy = tempVec[ii].yhit-ycenter;
2133  Float_t dz = tempVec[ii].zhit;
2134  mDeltaHitMatch2->Fill(dy, dz);
2135  }
2136  }
2137 
2138  erasedVec = tempVec;
2139  while (tempVec.size() != 0) {
2140  StructCellHit cellHit;
2141  Int_t nCells = 0;
2142  idVector vTrackId;
2143  vector<StThreeVectorD> vPosition;
2144  vector<Int_t> vchannel, vtray, vmodule, vcell;
2145  vector<Float_t> vzhit, vyhit;
2146 
2147  tofCellHitVectorIter tempIter=tempVec.begin();
2148  tofCellHitVectorIter erasedIter=erasedVec.begin();
2149  while(erasedIter!= erasedVec.end()) {
2150  if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
2151  nCells++;
2152  vchannel.push_back(erasedIter->channel);
2153  vtray.push_back(erasedIter->tray);
2154  vmodule.push_back(erasedIter->module);
2155  vcell.push_back(erasedIter->cell);
2156  vPosition.push_back(erasedIter->hitPosition);
2157  vTrackId.push_back(erasedIter->trackIdVec.back());
2158  vzhit.push_back(erasedIter->zhit);
2159  vyhit.push_back(erasedIter->yhit);
2160 
2161  erasedVec.erase(erasedIter);
2162  erasedIter--;
2163  }
2164  erasedIter++;
2165  }
2166 
2167  if (nCells==1){
2168  // for singly hit cell, copy data in singleHitCellsVec
2169  cellHit.channel = vchannel[0];
2170  cellHit.tray = vtray[0];
2171  cellHit.module = vmodule[0];
2172  cellHit.cell = vcell[0];
2173  cellHit.trackIdVec.push_back(vTrackId[0]);
2174  cellHit.hitPosition = vPosition[0];
2175  cellHit.matchFlag = 0;
2176  cellHit.zhit = vzhit[0];
2177  cellHit.yhit = vyhit[0];
2178 
2179  FinalMatchedCellsVec.push_back(cellHit);
2180 
2181  // debugging output
2182  if (Debug()) {
2183  LOG_INFO << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
2184  idVectorIter ij=vTrackId.begin();
2185  while (ij != vTrackId.end()) { LOG_INFO << " " << *ij; ij++; }
2186  LOG_INFO <<endm;
2187  }
2188  }
2189  else if (nCells>1){ // for multiple hit cells find the most likely candidate.
2190  Int_t thiscandidate(-99);
2191  Int_t thisMatchFlag(0);
2192 
2193  // sort on hitposition
2194  Float_t ss(99.);
2195  vector<Int_t> ssCandidates;
2196  thisMatchFlag = 2;
2197  if (Debug()) LOG_INFO << " ss " << endm;
2198  for (Int_t i=0;i<nCells;i++){
2199  Float_t yy = vyhit[i];
2200  Float_t ycell = (vcell[i]-1-2.5)*mWidthPad;
2201  Float_t ll = fabs(yy-ycell);
2202  if(ll<ss) {
2203  ss = ll;
2204  ssCandidates.clear();
2205  ssCandidates.push_back(i);
2206  }else if (ll==ss)
2207  ssCandidates.push_back(i);
2208  }
2209  if (ssCandidates.size()==1){
2210  thiscandidate = ssCandidates[0];
2211  Int_t daqId = vchannel[thiscandidate];
2212  if (Debug()) LOG_INFO << "candidate =" << daqId << endm;
2213  }
2214 
2215 
2216  if (thiscandidate>=0) {
2217  cellHit.channel = vchannel[thiscandidate];
2218  cellHit.tray = vtray[thiscandidate];
2219  cellHit.module = vmodule[thiscandidate];
2220  cellHit.cell = vcell[thiscandidate];
2221  cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
2222  cellHit.hitPosition = vPosition[thiscandidate];
2223  cellHit.matchFlag = thisMatchFlag;
2224  cellHit.zhit = vzhit[thiscandidate];
2225  cellHit.yhit = vyhit[thiscandidate];
2226 
2227  FinalMatchedCellsVec.push_back(cellHit);
2228 
2229  // debugging output
2230  if (Debug()) {
2231  LOG_INFO << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:" << vTrackId[thiscandidate] << endm;
2232  }
2233  }
2234 
2235  } else {
2236  LOG_INFO << "E: no cells belong to this track ... should not happen!" << endm;
2237  }
2238 
2239  tempVec = erasedVec;
2240  }
2241 
2242  LOG_INFO << "E: before/after: " << singleHitCellsVec.size() << "/" << FinalMatchedCellsVec.size() << endm;
2243  // end of Sect.E
2244 
2245  //.........................................................................
2246  // F. perform further selection and
2247  // fill valid track histograms, ntuples and CellCollection
2248  //
2249  tempVec.clear();
2250  tempVec = FinalMatchedCellsVec;
2251  if(mHisto) {
2252  if(FinalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
2253  mCellsPerEventMatch3->Fill(tempVec.size());
2254  for(unsigned int ii=0;ii<tempVec.size();ii++) {
2255  mTracksPerCellMatch3->Fill(tempVec[ii].trackIdVec.size());
2256 // mDaqOccupancyMatch3->Fill(tempVec[ii].channel);
2257  mDaqOccupancyMatch3->Fill((tempVec[ii].module-1)*mNCell+(tempVec[ii].cell-1));
2258  Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
2259  Float_t dy = tempVec[ii].yhit - ycenter;
2260  Float_t dz = tempVec[ii].zhit;
2261  mDeltaHitMatch3->Fill(dy, dz);
2262  }
2263  }
2264 
2265  StTofCellCollection *mCellCollection = new StTofCellCollection;
2266  Int_t nValidSingleHitCells(0), nValidSinglePrimHitCells(0);
2267 
2268  for (size_t ii=0; ii < FinalMatchedCellsVec.size(); ii++){
2269  Int_t daqId = FinalMatchedCellsVec[ii].channel;
2270  Int_t jj = daqId;
2271  Int_t tray = FinalMatchedCellsVec[ii].tray;
2272  Int_t module = FinalMatchedCellsVec[ii].module;
2273  Int_t cell = FinalMatchedCellsVec[ii].cell;
2274 
2275 // Float_t ycenter = (cell-1-2.5)*mWidthPad;
2276 // Float_t dy = FinalMatchedCellsVec[ii].yhit - ycenter;
2277  if (FinalMatchedCellsVec[ii].trackIdVec.size()!=1)
2278  LOG_INFO << "F: WHAT!?! mult.matched cell in single cell list " << daqId << endm;
2279 
2280 
2281  // Read in Leading and Trailing edge TDC, apply on INL correction
2282  //
2283  int tmptdc = (mSortTofRawData->GetLeadingTdc(tray,jj,kTRUE))[0];
2284  int bin = (int)tmptdc&0x3ff;
2285  double tmptdc_f = tmptdc + mTofINLCorr->getTrayINLCorr(tray, jj, bin);
2286  double letime = tmptdc_f*VHRBIN2PS / 1000.;
2287 
2288  tmptdc=(mSortTofRawData->GetTrailingTdc(tray,jj,kTRUE))[0];
2289  bin = (int)tmptdc&0x3ff;
2290  tmptdc_f = tmptdc + mTofINLCorr->getTrayINLCorr(tray, jj, bin);
2291  double tetime = tmptdc_f*VHRBIN2PS / 1000.;
2292 
2293  // get track-id from cell hit vector
2294  unsigned int trackNode = FinalMatchedCellsVec[ii].trackIdVec[0];
2295 // StTrack *theTrack = nodes[trackNode]->track(primary);
2296  StTrack *globalTrack = nodes[trackNode]->track(global);
2297 
2298  // 2. continue only if the (primary) track exists
2299 // if (validTofTrack(theTrack) && fabs(dy)<1.9 ){
2300  nValidSinglePrimHitCells++;
2301 
2302  //--- store number of hits per track
2303  // Int_t nHitsPerTrack = theTrack->topologyMap().numberOfHits(kTpcId);
2304 
2305  // select the apropriate track geometry
2306  // StTrackGeometry *theTrackGeometry = trackGeometry(globalTrack);
2307 
2308  //--- calculate local hit position on cell based on average hitposition
2309  // Float_t localHitPos = mTofGeom->cellHitPosition(&allMatchedCellsVec[ii].hitPosition);
2310 
2311  // Fill TOF Cell Collection
2312  StTofCell *tofCell = new StTofCell(tray, module, cell, daqId, globalTrack, FinalMatchedCellsVec[ii].zhit, FinalMatchedCellsVec[ii].matchFlag, FinalMatchedCellsVec[ii].hitPosition);
2313  tofCell->setLeadingEdgeTime(letime);
2314  tofCell->setTrailingEdgeTime(tetime);
2315  mCellCollection->push_back(tofCell);
2316 
2317  // dump debug data
2318  if (Debug()){
2319  LOG_INFO << "F: itray=" << tray << " imodule=" << module << " icell=" << cell << "\tnodeid:";
2320  idVectorIter ij=FinalMatchedCellsVec[ii].trackIdVec.begin();
2321  while (ij != FinalMatchedCellsVec[ii].trackIdVec.end()) { LOG_INFO << " " << *ij; ij++; }
2322  LOG_INFO << endm;
2323  }
2324 
2325 // } // track exists
2326  } // end final matched cells
2327 
2328 
2329  // put INL corrected vpd information in StTofCell as well
2330  StSPtrVecTofData &tofData = theTof->tofData();
2331  for(size_t ii = 0; ii<tofData.size(); ii++) {
2332  StTofData *aData = dynamic_cast<StTofData *>(tofData[ii]);
2333  if(!aData) continue;
2334 
2335  int dataIndex = aData->dataIndex();
2336  int trayId = dataIndex / mNTOF;
2337  if(trayId<120) continue; // only vpd selected
2338  int ewId = trayId - 120 + 1; // 1: west, 2: east
2339 
2340  int tubeId = dataIndex % mNTOF + 1;
2341  int lechan = (ewId==1) ? mDaqMap->WestPMT2TDIGLeChan(tubeId) : mDaqMap->EastPMT2TDIGLeChan(tubeId);
2342  int techan = (ewId==1) ? mDaqMap->WestPMT2TDIGTeChan(tubeId) : mDaqMap->EastPMT2TDIGTeChan(tubeId);
2343 
2344  int tmptdc = aData->leadingTdc();
2345  int bin = (int)tmptdc&0x3ff;
2346  double tmptdc_f = tmptdc + mTofINLCorr->getVpdINLCorr(ewId, lechan, bin);
2347  double letime = tmptdc_f*VHRBIN2PS / 1000.;
2348 
2349  tmptdc = aData->trailingTdc();
2350  bin = (int)tmptdc&0x3ff;
2351  tmptdc_f = tmptdc + mTofINLCorr->getVpdINLCorr(ewId, techan, bin);
2352  double tetime = tmptdc_f*VHRBIN2PS / 1000.;
2353 
2354  StThreeVectorF zero(0.,0.,0.);
2355  StTofCell *tofCell = new StTofCell(120+ewId, 0, tubeId, lechan, 0, 0, 0, zero);
2356  tofCell->setLeadingEdgeTime(letime);
2357  tofCell->setTrailingEdgeTime(tetime);
2358  mCellCollection->push_back(tofCell);
2359 
2360  if (Debug()){
2361  LOG_INFO << "F: itray=" << trayId << " imodule=" << 0 << " itube=" << tubeId << " letime=" << letime << " tetime=" << tetime << endm;
2362  }
2363  }
2364 
2365 
2366  storeMatchData(mCellCollection,theTof);
2367  delete mCellCollection;
2368 
2369  delete mSortTofRawData;
2370  mSortTofRawData = 0;
2371 
2372  //check StEvent collections --
2373  if (theTof->dataPresent())
2374  LOG_INFO << " TofCollection: raw data container present" << endm;
2375  if (theTof->cellsPresent()){
2376  LOG_INFO << " TofCollection: cell container present."<<endm;
2377  if (Debug()){
2378  StSPtrVecTofCell& tmpCellTofVec = theTof->tofCells();
2379  LOG_INFO << " # of matched cells " << tmpCellTofVec.size() << endm;
2380  for (size_t i = 0; i < tmpCellTofVec.size(); i++) {
2381  StTofCell* p = tmpCellTofVec[i];
2382  LOG_INFO << p->trayIndex() << " " << p->moduleIndex() << " " << p->cellIndex() << " " << p->trailingEdgeTime() << " " << p->leadingEdgeTime() << " " << p->associatedTrack() << " " << p->matchFlag() << " " << p->position() << endm;
2383  }
2384  }
2385  }
2386  //-- end check
2387 
2388  LOG_INFO << "F: before/after" << FinalMatchedCellsVec.size() << "/" <<nValidSinglePrimHitCells << endm;
2389  // end of Sect.F
2390 
2391  LOG_INFO << "#(cell tracks): " << allCellsHitVec.size()
2392  << " #(hit cells): " << FinalMatchedCellsVec.size()
2393  << "\n#(single hits): " << nSingleHitCells
2394  << " #(single valid hits): " << nValidSingleHitCells
2395  << " #(single prim valid hits): " << nValidSinglePrimHitCells
2396  << endm;
2397 
2398 
2399 
2400  if (doPrintMemoryInfo) {
2401  StMemoryInfo::instance()->snapshot();
2402  StMemoryInfo::instance()->print();
2403  }
2404  if (doPrintCpuInfo) {
2405  timer.stop();
2406  LOG_INFO << "CPU time for StTofrMatchMaker::Make(): "
2407  << timer.elapsedTime() << " sec\n" << endm;
2408  }
2409 
2410  LOG_INFO << "StTofrMatchMaker -- bye-bye" << endm;
2411 
2412 
2413 
2414  return kStOK;
2415 }
2416 
2417 //---------------------------------------------------------------------------
2418 // store local slat collection in StEvent's tofCollection
2419 Int_t StTofrMatchMaker::storeMatchData(StTofCellCollection *cellCollection,
2420  StTofCollection* tofCollection){
2421  if(!tofCollection){
2422  LOG_INFO << "Error: no TofCollection -- returning" << endm;
2423  return kStErr;
2424  }
2425 
2426  for (size_t j=0;j<cellCollection->size();j++){
2427  tofCollection->addCell(cellCollection->getCell(j));
2428  if (Debug())
2429  LOG_INFO << "storing " << j << " " << " tray:"
2430  << cellCollection->getCell(j)->trayIndex() << " module:"
2431  << cellCollection->getCell(j)->moduleIndex() << " cell:"
2432  << cellCollection->getCell(j)->cellIndex() << endm;
2433  }
2434  return kStOK;
2435 }
2436 
2437 
2438 //---------------------------------------------------------------------------
2439 // create a local copy of the raw tofp data tofData in StEvent's tofCollection
2440 Int_t StTofrMatchMaker::getTofData(StTofCollection* tofCollection){
2441  if (!tofCollection) return kStERR;
2442  LOG_INFO << " Read in tof data ... " << endm;
2443  StSPtrVecTofData &tofData = tofCollection->tofData();
2444 
2445  // perform consistency check
2446  bool dataOK(true);
2447  LOG_INFO << "TOF raw data consistency test ...";
2448 
2449  Int_t tofsize = tofData.size();
2450  Int_t nTOFr = 0;
2451  if(mYear4) {
2452  if(tofsize<184) {
2453  gMessMgr->Warning("The size of tofData is NOT 184!","OS");
2454  nTOFr = 72; // DAQ not updated
2455  } else {
2456  nTOFr = mNTOFR;
2457  }
2458  } else if(mYear3) {
2459  nTOFr = 72;
2460  } else if(mYear2) {
2461  nTOFr = 0;
2462  }
2463 
2464  for (Int_t i=0;i<nTOFr;i++){
2465  Int_t iAdc = mDaqMap->DaqChan2ADCChan(i);
2466  mTofrAdc[i] = tofData[iAdc]->adc();
2467  Int_t iTdc = mDaqMap->DaqChan2TDCChan(i);
2468  mTofrTdc[i] = tofData[iTdc]->tdc();
2469  }
2470 
2471  for (Int_t i=0;i<6;i++){
2472  mPvpdAdc[i] = tofData[42+i]->adc();
2473  mPvpdTdc[i] = tofData[42+i]->tdc();
2474  if (mYear3||mYear4)
2475  mPvpdAdcLoRes[i] = tofData[54+i]->adc();
2476  }
2477 
2478  if (!dataOK) return kStWarn;
2479 
2480  return kStOK;
2481 }
2482 
2483 
2484 //---------------------------------------------------------------------------
2485 // Book histograms and create ordered collection for easy manipulation
2486 void StTofrMatchMaker::bookHistograms(void){
2487 
2488  mEventCounterHisto = new TH1D("eventCounter","eventCounter",20,0,20);
2489 
2490  mADCTDCCorelation = new TH2D("ADCTDCCorelation","ADCTDCCorelation",200,0,200,200,0,200);
2491  mCellsMultInEvent = new TH1D("cellsPerEvent","cellsPerEvent",100,0,100);
2492  mHitsMultInEvent = new TH1D("hitsPerEvent","hitsPerEvent",100,0,100);
2493  mHitsMultPerTrack = new TH1D("hitsPerTrack","hitsPerTrack",10,0,10);
2494  mDaqOccupancy = new TH1D("daqOccupancy","daqOccupancy",192,0,192);
2495  mDaqOccupancyValid= new TH1D("daqOccupancyValid","daqOccupancyValid",192,0,192);
2496  mDaqOccupancyProj = new TH1D("daqOccupancyProj","daqOccupancyProj",192,0,192);
2497  mHitsPosition = new TH2D("hitsPosition","hitsPositions",300,-15.,15.,200,-5.,5.);
2498 
2499  mDaqOccupancyValidAll = new TH1D("daqOccupancyValidAll","",960,0,960);
2500  mDaqOccupancyProjAll = new TH1D("daqOccupancyProjAll","",960,0,960);
2501 
2502  mDaqOccupancyVpd = new TH1D("daqOccupancyVpd","",38,0,38);
2503  mDaqOccupancyValidVpd = new TH1D("daqOccupancyValidVpd","",38,0,38);
2504 
2505  // correlation
2506  for(int i=0;i<mNValidTrays_Run8;i++) {
2507  char hisname[100];
2508  sprintf(hisname,"Tray_%d",i+76);
2509  mHitCorr[i] = new TH2D(hisname,"",192,0,192,192,0,192);
2510  sprintf(hisname,"Tray_%d_module",i+76);
2511  mHitCorrModule[i] = new TH2D(hisname,"",32,0,32,32,0,32);
2512  }
2513  mHitCorrAll = new TH2D("Tray_All","",960,0,960,960,0,960);
2514 
2515  // TPC tracks
2516  if(mSaveTree) {
2517  mTrackTree = new TTree("track","track");
2518  mTrackTree->Branch("pt",&trackTree.pt,"pt/F");
2519  mTrackTree->Branch("eta",&trackTree.eta,"eta/F");
2520  mTrackTree->Branch("phi",&trackTree.phi,"phi/F");
2521  mTrackTree->Branch("nfitpts",&trackTree.nfitpts,"nfitpts/I");
2522  mTrackTree->Branch("dEdx",&trackTree.dEdx,"dEdx/F");
2523  mTrackTree->Branch("ndEdxpts",&trackTree.ndEdxpts,"ndEdxpts/I");
2524  mTrackTree->Branch("charge",&trackTree.charge,"charge/I");
2525  mTrackTree->Branch("projTrayId",&trackTree.projTrayId,"projTrayId/I");
2526  mTrackTree->Branch("projCellChan",&trackTree.projCellChan,"projCellChan/I");
2527  mTrackTree->Branch("projY",&trackTree.projY,"projY/F");
2528  mTrackTree->Branch("projZ",&trackTree.projZ,"projZ/F");
2529  }
2530 
2531  mTrackPtEta = new TH2D("trackPtEta","",100,0.,5.,60,-1.5,1.5);
2532  mTrackPtPhi = new TH2D("trackPtPhi","",100,0.,5.,120,0.,2*3.14159);
2533  mTrackNFitPts = new TH1D("trackNFitPts","",50,0.,50.);
2534  mTrackdEdxvsp = new TH2D("trackdEdxvsp","",500,0.,5.,1000,0.,10.);
2535  mNSigmaPivsPt = new TH2D("nSigmaPivsPt","",500,0.,5.,1000,-10.,10.);
2536 
2537  // primary association
2538  mCellsPerEventMatch1 = new TH1D("cellsPerEventMatch1","cellPerEventMatch1",100,0,100);
2539  mHitsPerEventMatch1 = new TH1D("hitsPerEventMatch1","hitsPerEventMatch1",100,0,100);
2540  mCellsPerTrackMatch1 = new TH1D("cellsPerTrackMatch1","cellsPerTrackMatch1",100,0,100);
2541  mTracksPerCellMatch1 = new TH1D("tracksPerCellMatch1","tracksPerCellMatch1",100,0,100);
2542  mDaqOccupancyMatch1 = new TH1D("daqOccupancyMatch1","daqOccupancyMatch1",192,0,192);
2543  mDeltaHitMatch1 = new TH2D("deltaHitMatch1","deltaHitMatch1",300,-15,15,200,-5.,5.);
2544 
2545  // kick out multi-hit
2546  mCellsPerEventMatch2 = new TH1D("cellsPerEventMatch2","cellPerEventMatch2",100,0,100);
2547  mHitsPerEventMatch2 = new TH1D("hitsPerEventMatch2","hitsPerEventMatch2",100,0,100);
2548  mCellsPerTrackMatch2 = new TH1D("cellsPerTrackMatch2","cellsPerTrackMatch2",100,0,100);
2549  mTracksPerCellMatch2 = new TH1D("tracksPerCellMatch2","tracksPerCellMatch2",100,0,100);
2550  mDaqOccupancyMatch2 = new TH1D("daqOccupancyMatch2","daqOccupancyMatch2",192,0,192);
2551  mDeltaHitMatch2 = new TH2D("deltaHitMatch2","deltaHitMatch2",300,-15,15,200,-5.,5.);
2552 
2553  // sort out multi matched cells
2554  mCellsPerEventMatch3 = new TH1D("cellsPerEventMatch3","cellPerEventMatch3",100,0,100);
2555  mHitsPerEventMatch3 = new TH1D("hitsPerEventMatch3","hitsPerEventMatch3",100,0,100);
2556  mCellsPerTrackMatch3 = new TH1D("cellsPerTrackMatch3","cellsPerTrackMatch3",100,0,100);
2557  mTracksPerCellMatch3 = new TH1D("tracksPerCellMatch3","tracksPerCellMatch3",100,0,100);
2558  mDaqOccupancyMatch3 = new TH1D("daqOccupancyMatch3","daqOccupancyMatch3",192,0,192);
2559  mDeltaHitMatch3 = new TH2D("deltaHitMatch3","deltaHitMatch3",300,-15,15,200,-5.,5.);
2560 
2561  return;
2562 }
2563 
2564 
2565 //---------------------------------------------------------------------------
2566 // store histograms in a seperate root file
2567 void StTofrMatchMaker::writeHistogramsToFile(){
2568  // Output file
2569  TFile *theHistoFile = new TFile(mHistoFileName.c_str(), "RECREATE");
2570  LOG_INFO << "StTofrMatchMaker::writeHistogramsToFile()"
2571  << " histogram file " << mHistoFileName << endm;
2572 
2573  theHistoFile->cd();
2574 
2575  if(mHisto) {
2576 
2577  mEventCounterHisto->Write();
2578  mADCTDCCorelation->Write();
2579  mCellsMultInEvent->Write();
2580  mHitsMultInEvent->Write();
2581  mHitsMultPerTrack->Write();
2582  mDaqOccupancy->Write();
2583  mDaqOccupancyValid->Write();
2584  mDaqOccupancyProj->Write();
2585  mHitsPosition->Write();
2586 
2587  mDaqOccupancyValidAll->Write();
2588  mDaqOccupancyProjAll->Write();
2589  mDaqOccupancyVpd->Write();
2590  mDaqOccupancyValidVpd->Write();
2591 
2592  for(int i=0;i<mNValidTrays_Run8;i++) {
2593  mHitCorr[i]->Write();
2594  mHitCorrModule[i]->Write();
2595  }
2596  mHitCorrAll->Write();
2597 
2598  mTrackPtEta->Write();
2599  mTrackPtPhi->Write();
2600  mTrackNFitPts->Write();
2601  mTrackdEdxvsp->Write();
2602  mNSigmaPivsPt->Write();
2603 
2604  mCellsPerEventMatch1->Write();
2605  mHitsPerEventMatch1->Write();
2606  mCellsPerTrackMatch1->Write();
2607  mTracksPerCellMatch1->Write();
2608  mDaqOccupancyMatch1->Write();
2609  mDeltaHitMatch1->Write();
2610 
2611  mCellsPerEventMatch2->Write();
2612  mHitsPerEventMatch2->Write();
2613  mCellsPerTrackMatch2->Write();
2614  mTracksPerCellMatch2->Write();
2615  mDaqOccupancyMatch2->Write();
2616  mDeltaHitMatch2->Write();
2617 
2618  mCellsPerEventMatch3->Write();
2619  mHitsPerEventMatch3->Write();
2620  mCellsPerTrackMatch3->Write();
2621  mTracksPerCellMatch3->Write();
2622  mDaqOccupancyMatch3->Write();
2623  mDeltaHitMatch3->Write();
2624 
2625  theHistoFile->Write();
2626  theHistoFile->Close();
2627 
2628  }
2629 
2630  if(mSaveTree) {
2631  TFile *theTreeFile = new TFile((mHistoFileName+".tree.root").c_str(), "RECREATE");
2632  theTreeFile->cd();
2633  mTrackTree->Write();
2634  theTreeFile->Write();
2635  theTreeFile->Close();
2636  }
2637 
2638  return;
2639 }
2640 
2641 
2642 //---------------------------------------------------------------------------
2643 // determine pVPD event type (strobe or beam)
2644 bool StTofrMatchMaker::strobeEvent(StSPtrVecTofData& tofData){
2645  // determine strobe event from pVPD TDC data
2646 
2647  Int_t nStrobedPvpdTdcs=0;
2648  for(Int_t i=0;i<mNPVPD;i++)
2649  if((tofData[42+i]->tdc()>mStrobeTdcMin[i]) &&
2650  (tofData[42+i]->tdc()<mStrobeTdcMax[i]))
2651  nStrobedPvpdTdcs++;
2652 
2653  if (nStrobedPvpdTdcs==mNPVPD) return true;
2654 
2655  return false;
2656 }
2657 
2658 
2659 //---------------------------------------------------------------------------
2660 // determine whether this is a valid TOF beam event
2661 bool StTofrMatchMaker::validEvent(StEvent *event){
2662  mEventCounter++;
2663  // 1. must have non-zero pointer
2664  if (!event) return false;
2665  if(mHisto) mEventCounterHisto->Fill(1);
2666 
2667  // 2. must have a valid primary vertex
2668 // if (!event->primaryVertex()) return false;
2669  mAcceptedEventCounter++;
2670  if(mHisto) mEventCounterHisto->Fill(2);
2671 
2672  // 3a. must have TOF collection
2673  if (!event->tofCollection()){
2674  LOG_INFO << "TOF is not present" << endm;
2675  return false;
2676  }
2677  if(mHisto) mEventCounterHisto->Fill(3);
2678 
2679  // 3b. must have TOF raw data available
2680  if (!(event->tofCollection()->dataPresent()||event->tofCollection()->rawdataPresent())){
2681  LOG_INFO << "TOF is present but no Raw Data" << endm;
2682  if (!(event->tofCollection()->cellsPresent())){
2683  LOG_INFO << " and no Cell Data" << endm;
2684  }
2685  return false;
2686  }
2687  mTofEventCounter++;
2688  if(mHisto) mEventCounterHisto->Fill(4);
2689 
2690 
2691 
2692  // 4. must be a TOF beam event, i.e. a non-strobe event
2693  if(event->tofCollection()->dataPresent()) {
2694  StSPtrVecTofData &tofData = event->tofCollection()->tofData();
2695  LOG_INFO << " tofData size = " << tofData.size() << endm;
2696  if (strobeEvent(tofData)){
2697  mTofStrobeEventCounter++;
2698  if (event->primaryVertex()) mAcceptAndStrobe++; // keep track of #valid strobed evts
2699  gMessMgr->Info("strobe event","OTS");
2700  return false;
2701  }
2702  }
2703  if(mHisto) mEventCounterHisto->Fill(5);
2704 
2705  mAcceptAndBeam++;
2706 
2707  // and we have a winner!
2708  gMessMgr->Info("TOF present ... and valid beam event","OTS");
2709 
2710  return true;
2711 }
2712 
2713 
2714 //---------------------------------------------------------------------------
2715 // determine whether this is a valid TPC track
2716 bool StTofrMatchMaker::validTrack(StTrack *track){
2717  // 1. no track, no go.
2718  if (!track) return false;
2719 
2720  // 2. track quality flag, should be >0
2721  if (track->flag()<=0) return false;
2722 
2723  // 3. minimum #hits per track
2724  if (track->topologyMap().numberOfHits(kTpcId) < mMinHitsPerTrack) return false;
2725  // 4. minimum #fit points per track
2726  if (track->fitTraits().numberOfFitPoints(kTpcId) < mMinFitPointsPerTrack) return false;
2727  // 5. minimum #fit points over #maximum points
2728  float ratio = (1.0*track->fitTraits().numberOfFitPoints(kTpcId)) / (1.0*track->numberOfPossiblePoints(kTpcId));
2729  if (ratio < mMinFitPointsOverMax) return false;
2730 
2731  return true;
2732 }
2733 
2734 //---------------------------------------------------------------------------
2735 // determine whether this is a valid TPC track
2736 bool StTofrMatchMaker::validTrackRun8(StGlobalTrack *track){
2737  // 1. no track, no go.
2738  if (!track) return false;
2739 
2740  // 2. track quality flag, should be >0
2741  if (track->flag()<=0) return false;
2742 
2743  // 3. minimum #hits per track
2744  if (track->topologyMap().numberOfHits(kTpcId) < mMinHitsPerTrack) return false;
2745  // 4. minimum #fit points per track
2746  if (track->fitTraits().numberOfFitPoints(kTpcId) < mMinFitPointsPerTrack) return false;
2747  // 5. chi2 per track
2748  float chi2 = track->fitTraits().chi2(0);
2749  if ( chi2 > 10. ) return false;
2750 
2751  // 6. error
2752  StDcaGeometry *dcaGeom = track->dcaGeometry();
2753  if(!dcaGeom) return false;
2754  float pt = dcaGeom->pt();
2755  float sigma_pT = sqrt(dcaGeom->errMatrix()[9])*pt*pt;
2756  if ( sigma_pT/pt > 0.3 ) return false;
2757 
2758  return true;
2759 }
2760 
2761 
2762 //---------------------------------------------------------------------------
2763 // determine whether this is a valid TOF track
2764 bool StTofrMatchMaker::validTofTrack(StTrack *track){
2765  // select valid tracks for time-of-flight calculations
2766 
2767  // 1. track must exist
2768  if (!track) return false;
2769 
2770  // 2. track quality flag, should be >0
2771  if (track->flag()<=0) return false;
2772 
2773  // 3. track must be a primary track
2774  if (!dynamic_cast<StPrimaryTrack*>(track)) return false;
2775 
2776  // 4. DCA cut (obsolete?)
2777  Double_t DCA= track->impactParameter();
2778  Int_t charge = track->geometry()->charge();
2779  if (DCA > mMaxDCA) {LOG_INFO << "dca>max:" << DCA<< endm; return false;}
2780  if (charge==0) { LOG_INFO << " neutral charge" << endm; return false; }
2781 
2782  return true;
2783 }
2784 
2785 
2786 //---------------------------------------------------------------------------
2787 // returns the proper track geometry, based on a global user setting
2788 StTrackGeometry* StTofrMatchMaker::trackGeometry(StTrack* track){
2789  // returns apropriate StTrackGeometry (standard or outerGeometry)
2790  if (!track) return 0;
2791  StTrackGeometry *thisTrackGeometry;
2792  if (mOuterTrackGeometry)
2793  thisTrackGeometry = track->outerGeometry();
2794  else
2795  thisTrackGeometry = track->geometry();
2796  return thisTrackGeometry;
2797 }
2798 
2799 /*****************************************************************
2800  *
2801  * $Log: StTofrMatchMaker.cxx,v $
2802  * Revision 1.34 2018/02/26 23:26:52 smirnovd
2803  * StTof: Remove outdated ClassImp macro
2804  *
2805  * Revision 1.33 2018/02/26 23:13:21 smirnovd
2806  * Move embedded CVS log messages to the end of file
2807  *
2808  * Revision 1.32 2012/12/17 22:57:29 geurts
2809  * bugfix (tickets #2456/#2457)
2810  *
2811  * Revision 1.31 2012/12/14 06:36:04 geurts
2812  * Changed global database calls to direct table access and/or removed deprecated database access code.
2813  *
2814  * Revision 1.30 2009/07/28 16:06:42 geurts
2815  * Bug Ticket #1591: explicit initialization of mStrobeTdcMin, mStrobeTdcMax, and mPedTOFr.
2816  *
2817  * Revision 1.29 2009/07/24 22:42:08 fine
2818  * replace the deprecated API of the STAR messenger
2819  *
2820  * Revision 1.28 2009/07/24 22:33:33 fine
2821  * Make the code C++ compliant
2822  *
2823  * Revision 1.27 2009/06/09 19:45:35 jeromel
2824  * Changes for BT#1428
2825  *
2826  * Revision 1.26 2008/07/23 19:22:03 dongx
2827  * New track quality cuts for Run8
2828  *
2829  * Revision 1.25 2008/06/24 21:58:13 dongx
2830  * fixed a bug of crashing due to potential empty track in trackNodes
2831  *
2832  * Revision 1.24 2008/06/19 16:11:54 dongx
2833  * fixed bug of filling an undefined histogram in case of production
2834  *
2835  * Revision 1.23 2008/05/12 17:16:37 dongx
2836  * letime and tetime in StTofCell stored in nano-seconds
2837  *
2838  * Revision 1.22 2008/05/06 18:41:39 dongx
2839  * - Fixed bug in ouput histogram filename switch
2840  * - Added switch for tpc track tree output
2841  *
2842  * Revision 1.21 2008/04/23 18:20:15 dongx
2843  * vpd letime and tetime stored in double precision
2844  *
2845  * Revision 1.20 2008/04/22 22:31:22 dongx
2846  * leadingEdgeTime and trailingEdgeTime stored as double precision in StTofCell
2847  *
2848  * Revision 1.19 2008/04/22 20:55:13 dongx
2849  * leadingEdgeTime and trailingEdgeTime stored as double precision in StTofCell
2850  *
2851  * Revision 1.18 2008/03/27 00:16:03 dongx
2852  * update for Run8 finished.
2853  *
2854  * Revision 1.17 2007/11/29 22:43:11 dongx
2855  * changed vpd trayId definition to 121 (East) and 122 (West)
2856  *
2857  * Revision 1.16 2007/11/28 02:17:08 dongx
2858  * trayId for vpd in StTofCell: 901 (E), 902 (W)
2859  * dataIndex for vpd in StTofData: use 121 and 122 as trayId
2860  *
2861  * Revision 1.15 2007/11/22 00:22:37 dongx
2862  * update for run8 - first version
2863  *
2864  * Revision 1.14 2007/04/17 23:02:20 dongx
2865  * replaced with standard STAR Loggers
2866  *
2867  * Revision 1.13 2007/02/28 23:31:59 dongx
2868  * completion for Run V matching
2869  * - trailing tdc and leading tdc stored as adc and tdc in StTofCell
2870  * - multi-hit association cases use hit position priority
2871  *
2872  * Revision 1.12 2005/04/12 17:31:56 dongx
2873  * update for year 5 data - not completed, leave as empty at present
2874  *
2875  * Revision 1.11 2004/08/11 18:57:44 dongx
2876  * deltay quality cut applied
2877  *
2878  * Revision 1.10 2004/08/10 19:25:13 dongx
2879  * updated to be compatible with RunII data
2880  *
2881  * Revision 1.9 2004/07/13 16:12:32 dongx
2882  * continuing update for the DAQ array in Run IV
2883  *
2884  * Revision 1.8 2004/07/12 17:21:19 dongx
2885  * fix the crashing when running events before day 030 (TofDaq not updated to 184 yet)
2886  *
2887  * Revision 1.7 2004/05/03 23:08:50 dongx
2888  * change according to the update of StTofrGeometry, save CPU time by 100 times
2889  *
2890  * Revision 1.5 2004/04/01 20:08:54 dongx
2891  * fix a bug about the hit position stored in TofCell
2892  *
2893  * Revision 1.4 2004/03/16 22:30:49 dongx
2894  * fix the warning message when compiling
2895  *
2896  * Revision 1.3 2004/03/11 22:30:34 dongx
2897  * -move m_Mode control to Init()
2898  * -clear up
2899  *
2900  * Revision 1.2 2004/03/09 17:44:56 dongx
2901  * first release
2902  *
2903  */
StMaker::m_Mode
Int_t m_Mode
counters
Definition: StMaker.h:81
Error
Definition: FJcore.h:367
StHelix::x
double x(double s) const
coordinates of helix at point s
Definition: StHelix.hh:182
StHelix::distance
double distance(const StThreeVector< double > &p, bool scanPeriods=true) const
minimal distance between point and helix
Definition: StHelix.cc:240
TDataSet::ls
virtual void ls(Option_t *option="") const
Definition: TDataSet.cxx:495
kStWarn
Definition: Stypes.h:42
kStOK
Definition: Stypes.h:40
kStErr
Definition: Stypes.h:44
TDataSet::Find
virtual TDataSet * Find(const char *path) const
Definition: TDataSet.cxx:362
kStERR
Definition: Stypes.h:45
Generated by   doxygen 1.8.5