StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StBTofMatchMaker.cxx
1 /*******************************************************************
2  *
3  * $Id: StBTofMatchMaker.cxx,v 1.21 2018/03/16 18:38:49 genevb Exp $
4  *
5  * Author: Xin Dong
6  *****************************************************************
7  *
8  * Description: BTof Match Maker to do the matching between the
9  * fired celles and TPC tracks
10  *
11  *****************************************************************
12  *
13  * $Log: StBTofMatchMaker.cxx,v $
14  * Revision 1.21 2018/03/16 18:38:49 genevb
15  * Use TGeo initializer for BTof geometry
16  *
17  * Revision 1.20 2017/10/20 17:50:33 smirnovd
18  * Squashed commit of the following:
19  *
20  * StBTof: Remove outdated ClassImp macro
21  *
22  * Prefer explicit namespace for std:: names in header files
23  *
24  * Removed unnecessary specification of default std::allocator
25  *
26  * Frank signed-off
27  *
28  * Revision 1.19 2012/05/07 14:11:16 fisyak
29  * Keep btofGeometry in const area for future use
30  *
31  * Revision 1.18 2011/08/04 19:14:02 geurts
32  * Bug fix: allow ideal geometry setting in the case an alignment file is used [Patrick]
33  *
34  * Revision 1.17 2011/07/27 16:13:58 geurts
35  * Alignment calibration modifications [Patrick Huck]:
36  * - modified to open the local Z window cut to determine the z offset
37  * - variables mZLocalCut, mCalculateAlign and mAlignFileName added
38  * - functions setCalculateAlign and setAlignFileName added
39  *
40  * Revision 1.16 2010/08/09 19:18:45 geurts
41  * Include local theta calculation in CellHit structure. Pass LocalTheta info on to TOF PID traits. [Masa]
42  *
43  * Revision 1.15 2010/07/16 04:25:16 geurts
44  * initialize mUseIdealGeometry to be kFALSE in ctor
45  *
46  * Revision 1.14 2010/07/14 20:35:21 geurts
47  * introduce switch to enable ideal MC geometry, without alignment updates. Default: disabled
48  *
49  * Revision 1.13 2010/05/25 22:09:38 geurts
50  * improved database handling and reduced log output
51  *
52  * Revision 1.12 2010/05/20 22:58:47 geurts
53  * Keep BTofMatchMaker from crashing ungracefully when no mEvent or BTOF Collection is found
54  *
55  * Revision 1.11 2010/03/22 19:40:28 dongx
56  * Fixed a bug in setting index2Primary in processMuDst
57  *
58  * Revision 1.10 2010/03/19 22:25:39 dongx
59  * - Added getBTofGeom() function for outside use
60  * - Remove AddConst(btofGeometry) to avoid crash due to duplication
61  * - TOT selection window opened to 40 ns
62  * - Added CPU timer printouts for processStEvent() funciton
63  *
64  * Revision 1.9 2010/03/04 21:59:27 dongx
65  * Further addition in the initial clean up for primary tracks too
66  *
67  * Revision 1.8 2010/03/04 21:40:54 dongx
68  * Added clean up in processMuDst to remove associations done before
69  *
70  * Revision 1.7 2010/03/04 03:54:28 dongx
71  * Further improvement on the processMuDst speed in Step F.
72  *
73  * Revision 1.6 2010/03/04 00:08:47 dongx
74  * Removed primary check for globals at projection in accessing MuDst function as it takes significant CPU time. Waiting for MuDst update
75  * Added some more debugging output for CPU time usage
76  *
77  * Revision 1.5 2010/02/25 05:17:10 dongx
78  * Geometry initalization moved from Init() to InitRun()
79  *
80  * Revision 1.4 2009/09/15 00:30:45 dongx
81  * 1) Added the functionality to perform the matching with MuDst directly.
82  * 2) Several updates on the track cuts used for matching
83  * - flag<1000 was added
84  * - nHits>15 cut was removed
85  * 3) Created a new StBTofPidTraits for any primary track
86  * 4) Local Z window cut set to symmetric (fabs(localz)<3.05)
87  * 5) Some small changes in the LOGGER output.
88  *
89  * Revision 1.3 2009/08/26 20:33:56 dongx
90  * Geometry init moved to Init() function, also allow reading in from others
91  *
92  * Revision 1.2 2009/07/24 18:52:53 dongx
93  * - Local Z window restricted in the projection
94  * - ToT selection is used firstly when more than one hits associated with a track
95  * - matchFlag updated
96  * 0: no matching
97  * 1: 1-1 matching
98  * 2: 1-2 matching, pick up the one with higher ToT value (<25ns)
99  * 3: 1-2 matching, pick up the one with closest projection posision along y
100  *
101  * Revision 1.1 2009/06/23 13:15:03 geurts
102  * *** empty log message ***
103  *
104  *
105  *******************************************************************/
106 #include <iostream>
107 #include "StEventTypes.h"
108 #include "Stypes.h"
109 #include "StThreeVectorD.hh"
110 #include "StThreeVectorF.hh"
111 #include "StHelix.hh"
112 #include "StTrackGeometry.h"
113 #include "StDcaGeometry.h"
114 #include "StDedxPidTraits.h"
115 #include "StTrackPidTraits.h"
116 #include "StBTofPidTraits.h"
117 #include "StarClassLibrary/StParticleTypes.hh"
118 #include "StarClassLibrary/StParticleDefinition.hh"
119 #include "StTpcDedxPidAlgorithm.h"
120 #include "StEventUtilities/StuRefMult.hh"
121 #include "PhysicalConstants.h"
122 #include "StPhysicalHelixD.hh"
123 #include "StHelix.hh"
124 #include "StBTofCollection.h"
125 #include "StBTofUtil/tofPathLength.hh"
126 #include "StBTofUtil/StBTofGeometry.h"
127 #include "StBTofUtil/StBTofHitCollection.h"
128 #include "tables/St_tofConfig_Table.h"
129 #include "tables/St_tofTrayConfig_Table.h"
130 #include "TH1.h"
131 #include "TH2.h"
132 #include "TFile.h"
133 #include "TTree.h"
134 #include "StMessMgr.h"
135 #include "StMemoryInfo.hh"
136 #include "StTimer.hh"
137 #include "TGeoManager.h"
138 
139 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
140 #include "StMuDSTMaker/COMMON/StMuDst.h"
141 #include "StMuDSTMaker/COMMON/StMuBTofHit.h"
142 #include "StMuDSTMaker/COMMON/StMuTrack.h"
143 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
144 #include "StMuDSTMaker/COMMON/StMuBTofPidTraits.h"
145 
146 #include "StBTofMatchMaker.h"
147 //#include "TMemStat.h"
148 
149 
150 //---------------------------------------------------------------------------
151 StBTofMatchMaker::StBTofMatchMaker(const Char_t *name): StMaker(name){
152  // set default values
153  mEventCounter = 0;
154  mAcceptedEventCounter = 0;
155  mTofEventCounter = 0;
156  mAcceptAndBeam = 0;
157 
158  mBTofGeom = 0;
159 
160  mWidthPad = 3.45;
161  mZLocalCut = 3.05;
162 
163  setOuterTrackGeometry();
164  setMinHitsPerTrack(15);
165  setMinFitPointsPerTrack(15);
166  setMinFitPointsOverMax(0.52);
167  setMaxDCA(9999.);
168 
169  setCreateHistoFlag(kFALSE);
170  setHistoFileName("tofana.root");
171  setCreateTreeFlag(kFALSE);
172  setSaveGeometry(kFALSE);
173  mInitFromOther = kFALSE;
174  mUseIdealGeometry = kFALSE;
175  mCalculateAlign = kFALSE;
176  setAlignFileName("");
177  doPrintMemoryInfo = kFALSE;
178  doPrintCpuInfo = kFALSE;
179 
180  mEvent = 0;
181  mMuDst = 0;
182  mMuDstIn = kFALSE;
183 }
184 
185 StBTofMatchMaker::~StBTofMatchMaker(){ /* nope */}
186 
187 //void StBTofMatchMaker::Clear(Option_t *opt){StMaker::Clear();}
188 
189 //---------------------------------------------------------------------------
190 Int_t StBTofMatchMaker::Init(){
191  LOG_INFO << "Initializing match settings:" << endm;
192  LOG_INFO << " Minimum fitpoints per track: " << mMinFitPointsPerTrack << endm;
193  LOG_INFO << " Maximum DCA: " << mMaxDCA << endm;
194 
195  if (!mOuterTrackGeometry) {
196  LOG_WARN << "using standard trackgeometry()" << endm;
197  }
198 
199  // m_Mode can be set by SetMode() method
200  if(m_Mode) {
201 // setHistoFileName("tofana.root");
202  } else {
203  setHistoFileName("");
204  }
205 
206  if (mHisto){
207  bookHistograms();
208  LOG_INFO << "Histograms are booked" << endm;
209  if (mHistoFileName!="") {
210  LOG_INFO << "Histograms will be stored in " << mHistoFileName.c_str() << endm;
211  }
212  }
213 
214  // reset event counters
215  mEventCounter = 0;
216  mAcceptedEventCounter = 0;
217  mTofEventCounter = 0;
218  mAcceptAndBeam = 0;
219 
220  // for alignment calculate, we start from ideal geometry
221  if(mCalculateAlign) {
222  if (mAlignFileName=="") mUseIdealGeometry = kTRUE;
223  mZLocalCut = 5.0;
224  }
225 
226  return kStOK;
227 }
228 
229 
230 //---------------------------------------------------------------------------
231 Int_t StBTofMatchMaker::InitRun(Int_t runnumber){
232 
233  // determine TOF configuration from run#
234  LOG_INFO << "Initializing BTOF Geometry:" << endm;
236  // TOF geometry initializtion -- from GEANT geometry directly
237  // need St_geant_Maker be loaded before
239  StTimer geomTimer;
240  geomTimer.start();
241  mBTofGeom = 0;
242  if(TDataSet *geom = GetDataSet("btofGeometry")) {
243  mBTofGeom = (StBTofGeometry *)geom->GetObject();
244  LOG_INFO << " Found btofGeometry ... " << endm;
245  mInitFromOther = kTRUE;
246  } else {
247  mBTofGeom = new StBTofGeometry("btofGeometry","btofGeometry in MatchMaker");
248  LOG_INFO << " Create a new btofGeometry ... " << endm;
249  AddConst(new TObjectSet("btofGeometry",mBTofGeom));
250  }
251  if(mBTofGeom && !mBTofGeom->IsInitDone()) {
252  LOG_INFO << " BTofGeometry initialization ... " << endm;
253  //fg if(runnumber<1000000) mBTofGeom->SetMCOn();
254  if (mUseIdealGeometry) mBTofGeom->SetMCOn();
255  else mBTofGeom->SetMCOff();
256  LOG_INFO << " Alignment file: " << mAlignFileName.c_str() << endm;
257  mBTofGeom->SetAlignFile(mAlignFileName.c_str());
258  TVolume *starHall = gGeoManager ? nullptr : (TVolume *) GetDataSet("HALL");
259  mBTofGeom->Init(this, starHall, gGeoManager);
260  }
261 
262  geomTimer.stop();
263  LOG_INFO << "CPU time for StBTofGeometry initialization : " << geomTimer.elapsedTime() << " sec" << endm;
264 
265  return kStOK;
266 }
267 
268 //----------------------------------------------------------------------------
269 Int_t StBTofMatchMaker::FinishRun(Int_t runnumber){
270 
271  LOG_INFO << "StBTofMatchMaker -- cleaning up geometry (Finish)" << endm;
272  if (!mInitFromOther && mBTofGeom) delete mBTofGeom;
273  mBTofGeom=0;
274 
275  return kStOK;
276 }
277 
278 
279 //---------------------------------------------------------------------------
280 Int_t StBTofMatchMaker::Finish(){
281 
282  LOG_DEBUG << "StBTofMatchMaker ----- RUN SUMMARY ----- (Finish)\n"
283  << "\tProcessed " << mEventCounter << " events."
284  << " Accepted " << mAcceptedEventCounter << " events."
285  << " Rejected " << mEventCounter - mAcceptedEventCounter << " events\n"
286  << "\tTOF events " << mTofEventCounter
287  << "\t Accept & Beam " << mAcceptAndBeam << " events" << endm;
288 
289  //if (mHisto) writeHistogramsToFile();
290  if (mHistoFileName!="") writeHistogramsToFile();
291 
292 
293  return kStOK;
294 }
295 
296 //---------------------------------------------------------------------------
297 Int_t StBTofMatchMaker::Make(){
298  LOG_INFO << "StBTofMatchMaker -- welcome" << endm;
299 
300  if(mMuDstIn) processMuDst();
301  else processStEvent();
302 
303  return kStOK;
304 }
305 
306 //---------------------------------------------------------------------------
307 void StBTofMatchMaker::processStEvent(){
308 
309  if(mHisto) mEventCounterHisto->Fill(0);
310  // event selection ...
311  mEvent = (StEvent *) GetInputDS("StEvent");
312  if(!mEvent || !(mEvent->btofCollection()) || !(mEvent->btofCollection()->hitsPresent()) ) {
313  if (!mEvent) {LOG_INFO << "no StEvent" << endm;}
314  else
315  if (!(mEvent->btofCollection())) {LOG_INFO << "no BTof Collection" << endm;}
316  else
317  if (!(mEvent->btofCollection()->hitsPresent()) ) LOG_INFO << "no BTOF hits present" << endm;
318  LOG_INFO << "StBTofMatchMaker -- nothing to do ... bye-bye" << endm;
319  return;
320  }
321  if(mHisto) mEventCounterHisto->Fill(1);
322 
323  // timing & memory info -only when requested-
324  StTimer timer;
325  if (doPrintCpuInfo) timer.start();
326  if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
327 
328  //.........................................................................
329  // check for tofCollection and fill local copy with ADC and TDC data
330  StBTofCollection *theTof = mEvent->btofCollection();
331 
332  //.........................................................................
333  // read data from StBTofHit
335  //
336  tofCellHitVector daqCellsHitVec;
337  idVector validModuleVec;
338 
339  // multi-tray system
340  StSPtrVecBTofHit& tofHits = theTof->tofHits();
341  LOG_INFO << " Number of BTOF Hits = " << tofHits.size() << endm;
342  for(size_t i=0;i<tofHits.size();i++) {
343  //fg StBTofHit* aHit = dynamic_cast<StBTofHit *>(tofHits[i]);
344  StBTofHit* aHit = tofHits[i];
345  if(!aHit) continue;
346  if(aHit->tray()<=0||aHit->tray()>mNTray) continue; // barrel tray hits
347 
348  int trayId = aHit->tray();
349  int moduleId = aHit->module();
350  int cellId = aHit->cell();
351 
352  if(Debug()) { LOG_INFO <<"A: fired hit in " << " tray="<< trayId << " module="<<moduleId<<" cell="<<cellId<<endm; }
353 
354  StructCellHit aDaqCellHit;
355  aDaqCellHit.tray = trayId;
356  aDaqCellHit.module = moduleId;
357  aDaqCellHit.cell = cellId;
358  aDaqCellHit.tot = aHit->tot();
359  aDaqCellHit.index2BTofHit = i;
360  daqCellsHitVec.push_back(aDaqCellHit);
361 
362  // additional valid number configuration
363  int id = trayId*100+moduleId;
364  //fg bool ifind = kFALSE;
365  //fg for(size_t im=0;im<validModuleVec.size();im++) {
366  //fg if(id==validModuleVec[im]) {
367  //fg ifind = kTRUE;
368  //fg break;
369  //fg }
370  //fg }
371  //fg if(!ifind) validModuleVec.push_back(id);
372  if (find(validModuleVec.begin(), validModuleVec.end(), id) == validModuleVec.end())
373  validModuleVec.push_back(id);
374 
375  if(mHisto) {
376  mDaqOccupancy[trayId-1]->Fill((moduleId-1)*mNCell+(cellId-1));
377  }
378  }
379 
380  // end of Sect.A
381  if(Debug()) {
382  LOG_INFO << " total # of cells = " << daqCellsHitVec.size() << endm;
383  for(size_t iv = 0;iv<validModuleVec.size();iv++) {
384  LOG_DEBUG << " module # " << validModuleVec[iv] << " Valid! " << endm;
385  }
386  }
387  if(mHisto) {
388  mCellsMultInEvent->Fill(daqCellsHitVec.size());
389  if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
390  }
391 // if(!daqCellsHitVec.size()) return;
392  if (doPrintCpuInfo) {
393  timer.stop();
394  LOG_INFO << "CPU time after Step A - loading hits : "
395  << timer.elapsedTime() << " sec" << endm;
396  timer.start();
397  }
398 
399  //.........................................................................
401  //
402  tofCellHitVector allCellsHitVec;
403  StructCellHit cellHit;
404 
405  StSPtrVecTrackNode& nodes = mEvent->trackNodes();
406  Int_t nAllTracks=0;
407  Int_t nPrimaryHits = 0;
408  for (unsigned int iNode=0; iNode<nodes.size(); iNode++){
409  tofCellHitVector cellHitVec;
410  // cellHitVec.clear();
411  StGlobalTrack *theTrack = dynamic_cast<StGlobalTrack*>(nodes[iNode]->track(global));
412  if(!theTrack) continue;
413 
414  bool isPrimary = kFALSE;
415  StPrimaryTrack *pTrack = dynamic_cast<StPrimaryTrack*>(theTrack->node()->track(primary));
416  if(pTrack) isPrimary = kTRUE;
417 
418  StThreeVectorF mom = theTrack->geometry()->momentum();
419  float pt = mom.perp();
420  float eta = mom.pseudoRapidity();
421  float phi = mom.phi();
422  //fg if (phi<0.) phi += 2.*3.14159;
423  if (phi<0.) phi += 2.*M_PI;
424 
425  float dEdx = -999.;
426  int ndEdxpts = 0;
427  float nSigmaPion = -999.;
428  static StTpcDedxPidAlgorithm PidAlgorithm;
429  static StPionPlus* Pion = StPionPlus::instance();
430 
431  const StParticleDefinition* pd = theTrack->pidTraits(PidAlgorithm);
432  if (pd) {
433  nSigmaPion = PidAlgorithm.numberOfSigma(Pion);
434  }
435  if ( PidAlgorithm.traits() ) {
436  dEdx = PidAlgorithm.traits()->mean();
437  ndEdxpts = PidAlgorithm.traits()->numberOfPoints();
438  }
439  int nfitpts = theTrack->fitTraits().numberOfFitPoints(kTpcId);
440 
441  // make sure we have a track, a miniDST might have removed it...
442  if (validTrack(theTrack)){
443  if(mHisto) {
444  mTrackPtEta->Fill(pt, eta);
445  mTrackPtPhi->Fill(pt, phi);
446  mTrackNFitPts->Fill(nfitpts);
447  if(dEdx>0.) mTrackdEdxvsp->Fill(mom.mag(), dEdx*1.e6);
448  if(fabs(nSigmaPion)<5.) mNSigmaPivsPt->Fill(pt, nSigmaPion+5.*theTrack->geometry()->charge());
449  }
450 
451  if(mSaveTree) {
452  trackTree.pt = pt;
453  trackTree.eta = eta;
454  trackTree.phi = phi;
455  trackTree.nfitpts = nfitpts;
456  trackTree.dEdx = dEdx*1.e6;
457  trackTree.ndEdxpts = ndEdxpts;
458  trackTree.charge = theTrack->geometry()->charge();
459  trackTree.projTrayId = 0;
460  trackTree.projCellChan = -1;
461  trackTree.projY = -999.;
462  trackTree.projZ = -999.;
463  }
464 
465  nAllTracks++;
466  StPhysicalHelixD theHelix = trackGeometry(theTrack)->helix();
467 
468 // IntVec projTrayVec;
469 // if(!mBTofGeom->projTrayVector(theHelix, projTrayVec)) continue;
470 
471  IntVec idVec;
472  DoubleVec pathVec, thetaVec;
473  PointVec crossVec;
474 
475 // idVec.clear();
476 // pathVec.clear();
477 // crossVec.clear();
478 
479  Int_t ncells = 0;
480  if(mBTofGeom->HelixCrossCellIds(theHelix,idVec,pathVec,crossVec,thetaVec) ) {
481 // if(mBTofGeom->HelixCrossCellIds(theHelix, validModuleVec, projTrayVec, idVec, pathVec, crossVec)) {
482  Int_t cells = idVec.size();
483  for (Int_t i=0; i<cells; i++) {
484  Int_t icell,imodule,itray;
485  Double_t local[3],global[3];
486  for(Int_t i2=0;i2<3;i2++){
487  local[i2]=0;
488  }
489  global[0]=crossVec[i].x();
490  global[1]=crossVec[i].y();
491  global[2]=crossVec[i].z();
492  mBTofGeom->DecodeCellId(idVec[i], icell, imodule, itray);
493 // LOG_INFO << " decode " << idVec[i] << " to tray#" << itray << " module#" << imodule << " cell#" << icell << endm;
494  StBTofGeomSensor* sensor = mBTofGeom->GetGeomSensor(imodule,itray);
495  if(!sensor) {
496  LOG_WARN << " No sensitive module in the projection??? -- Something weird!!! " << endm;
497  continue;
498  }
499  sensor->Master2Local(&global[0],&local[0]);
500  icell = sensor->FindCellIndex(local);
501  // StThreeVectorD glo=sensor->GetCenterPosition();
502  StThreeVectorD glo(global[0], global[1], global[2]);
503  StThreeVectorD hitPos(local[0], local[1], local[2]);
504 // delete sensor; /// function in StBTofGeometry modified. Don't delete here.
505  if (fabs(local[2])<mZLocalCut) { // to be consistent with GEANT geometry
506  // and the alignment calibration
507 // if (local[2]<=2.7&&local[2]>=-3.4) {
508 // if (fabs(local[2])<4.) { // loose local z cut at first step
509  ncells++;
510  cellHit.tray = itray;
511  cellHit.module = imodule;
512  cellHit.cell = icell;
513  cellHit.trackIdVec.push_back((Int_t)iNode);
514  cellHit.hitPosition = glo; // global position
515  cellHit.zhit = (Float_t)hitPos.z();
516  cellHit.yhit = (Float_t)hitPos.y();
517  cellHit.theta = (Double_t)thetaVec[i];
518  cellHitVec.push_back(cellHit);
519  allCellsHitVec.push_back(cellHit);
520 
521  if(isPrimary) nPrimaryHits++;
522 
523  if(mHisto) {
524  mDaqOccupancyProj[itray-1]->Fill((imodule-1)*mNCell+(icell-1));
525  mHitsPosition->Fill(hitPos.y(), hitPos.z());
526  }
527 
528  if(mSaveTree) {
529  trackTree.projTrayId = itray;
530  trackTree.projCellChan = (imodule-1)*mNCell+(icell-1);
531  trackTree.projY = local[1];
532  trackTree.projZ = local[2];
533  }
534 
535  if(Debug()) {
536  LOG_DEBUG <<"B: nodeid=" << iNode << " projected in " << " tray="<< itray << " module="<<imodule<<" cell="<<icell<<endm;
537  LOG_DEBUG <<" hit position " << hitPos << endm;
538  // LOG_DEBUG <<" momemtum= " << pt << " " << eta << " " << phi << endm;
539  }
540  }
541  } // for (Int_t i=0...)
542  } // endif(helixcross...)
543  if(ncells>0&&mHisto) mHitsMultPerTrack->Fill(ncells);
544 
545  if(mHisto && mSaveTree) mTrackTree->Fill();
546 
547  } // if(ValidTrack)..
548  } // loop over nodes
549  if(Debug()) { LOG_INFO << "B: matched/available/total #tracknodes: " <<allCellsHitVec.size() << "/" <<nAllTracks << "/" << nodes.size() << endm; }
550  if(mHisto) {
551  mHitsMultInEvent->Fill(allCellsHitVec.size());
552  mHitsPrimaryInEvent->Fill(nPrimaryHits);
553  if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
554  }
555  // end of Sect.B
556  if (doPrintCpuInfo) {
557  timer.stop();
558  LOG_INFO << "CPU time after Step B - projection : "
559  << timer.elapsedTime() << " sec" << endm;
560  timer.start();
561  }
562 
563  //.........................................................................
565  //
566  tofCellHitVector matchHitCellsVec;
567 
568  tofCellHitVectorIter daqIter = daqCellsHitVec.begin();
569  for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
570  tofCellHitVectorIter proIter = allCellsHitVec.begin();
571  for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
572 
573  int daqIndex = (daqIter->module-1)*6 + (daqIter->cell-1);
574  int proIndex = (proIter->module-1)*6 + (proIter->cell-1);
575  int hisIndex = daqIter->tray - 1;
576 // int daqAllIndex = (daqIter->tray - 1)*mNTOF + daqIndex;
577 // int proAllIndex = (proIter->tray - 1)*mNTOF + proIndex;
578  if(daqIter->tray==proIter->tray) {
579  if (mHisto) {
580  if(hisIndex>=0&&hisIndex<mNTray) {
581  mHitCorr[hisIndex]->Fill(proIndex,daqIndex);
582  mHitCorrModule[hisIndex]->Fill(proIter->module-1,daqIter->module-1);
583  } else {
584  LOG_WARN << " weird tray # " << daqIter->tray << endm;
585  }
586  }
587  }
588  if( (daqIter->tray==proIter->tray)&&
589  (daqIter->module==proIter->module) &&
590  ( ( (proIter->cell==6)&&((proIter->cell==daqIter->cell) ||
591  (proIter->cell==daqIter->cell+1)) )
592  || ( (proIter->cell==1)&&((proIter->cell==daqIter->cell) ||
593  (proIter->cell==daqIter->cell-1)) )
594  || ( (proIter->cell>=2&&proIter->cell<=6) &&
595  ( (proIter->cell==daqIter->cell) ||
596  (proIter->cell==daqIter->cell-1) ||
597  (proIter->cell==daqIter->cell+1) ) ) ) ) {
598  cellHit.tray = daqIter->tray;
599  cellHit.module = daqIter->module;
600  cellHit.cell = daqIter->cell;
601  cellHit.hitPosition = proIter->hitPosition;
602  cellHit.trackIdVec = proIter->trackIdVec;
603  cellHit.zhit = proIter->zhit;
604  cellHit.yhit = proIter->yhit;
605  cellHit.tot = daqIter->tot;
606  cellHit.index2BTofHit = daqIter->index2BTofHit;
607  cellHit.theta = proIter->theta;
608  matchHitCellsVec.push_back(cellHit);
609  }
610  }
611  } //end {sec. C}
612  if(Debug()) { LOG_INFO << "C: before/after: " << allCellsHitVec.size() << "/" << matchHitCellsVec.size() << endm; }
613  if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
614  if (doPrintCpuInfo) {
615  timer.stop();
616  LOG_INFO << "CPU time after Step C - matching : "
617  << timer.elapsedTime() << " sec" << endm;
618  timer.start();
619  }
620 
621  //.........................................................................
623  //
624  Int_t nSingleHitCells(0);
625  Int_t nMultiHitsCells(0);
626 
627  tofCellHitVector singleHitCellsVec;
628  tofCellHitVector multiHitsCellsVec;
629 
630  tofCellHitVector tempVec = matchHitCellsVec;
631  tofCellHitVector erasedVec = tempVec;
632  while (tempVec.size() != 0) {
633  Int_t nTracks = 0;
634  idVector trackIdVec;
635 
636  tofCellHitVectorIter tempIter=tempVec.begin();
637  tofCellHitVectorIter erasedIter=erasedVec.begin();
638  while(erasedIter!= erasedVec.end()) {
639  if(tempIter->tray == erasedIter->tray &&
640  tempIter->module == erasedIter->module &&
641  tempIter->cell == erasedIter->cell) {
642  nTracks++;
643  trackIdVec.push_back(erasedIter->trackIdVec.back()); // merge
644  erasedVec.erase(erasedIter);
645  erasedIter--;
646  }
647  erasedIter++;
648  }
649 
650  cellHit.cell = tempIter->cell;
651  cellHit.module = tempIter->module;
652  cellHit.tray = tempIter->tray;
653  cellHit.hitPosition = tempIter->hitPosition;
654  cellHit.trackIdVec = trackIdVec;
655  cellHit.zhit = tempIter->zhit;
656  cellHit.yhit = tempIter->yhit;
657  cellHit.tot = tempIter->tot;
658  cellHit.index2BTofHit = tempIter->index2BTofHit;
659  cellHit.theta = tempIter->theta;
660 
661  Float_t ycenter = (tempIter->cell-1-2.5)*mWidthPad;
662  Float_t dy = tempIter->yhit - ycenter;
663  Float_t dz = tempIter->zhit;
664 
665  if(mHisto) {
666  mTracksPerCellMatch1->Fill(trackIdVec.size());
667 // mDaqOccupancyMatch1->Fill((tempIter->module-1)*mNCell+(tempIter->cell-1));
668  mDeltaHitMatch1->Fill(dy, dz);
669  }
670 
671  if (nTracks==1){
672  nSingleHitCells++;
673  singleHitCellsVec.push_back(cellHit);
674  } else if (nTracks>1){
675  nMultiHitsCells++;
676  multiHitsCellsVec.push_back(cellHit);
677  // for multiple hit cells either discard (yes) or
678  // find the most likely candidate.
679  } else {
680  LOG_WARN << "D: no tracks extrapolate to matched cell ... should not happen!" << endm;
681  }
682 
683  if(Debug()) {
684  LOG_DEBUG << "D: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
685  idVectorIter ij=trackIdVec.begin();
686  while (ij != trackIdVec.end()) { LOG_DEBUG << " " << *ij; ij++; }
687  LOG_DEBUG << endm;
688  }
689 
690  tempVec = erasedVec;
691  }
692  if(Debug()) { LOG_INFO << "D: before/after: " << matchHitCellsVec.size() << "/" << singleHitCellsVec.size() << endm; }
693  //end of Sect.D
694 
695  if(mHisto) {
696  mCellsPerEventMatch1->Fill(singleHitCellsVec.size()+multiHitsCellsVec.size());
697  if(singleHitCellsVec.size()) mEventCounterHisto->Fill(9);
698  }
699  if (doPrintCpuInfo) {
700  timer.stop();
701  LOG_INFO << "CPU time after Step D - erasing : "
702  << timer.elapsedTime() << " sec" << endm;
703  timer.start();
704  }
705 
706  //.........................................................................
708  //
709  tofCellHitVector FinalMatchedCellsVec;
710  // FinalMatchedCellsVec.clear();
711  tempVec = singleHitCellsVec;
712  if(mHisto) {
713  mCellsPerEventMatch2->Fill(tempVec.size());
714  for(unsigned int ii=0;ii<tempVec.size();ii++) {
715  mTracksPerCellMatch2->Fill(tempVec[ii].trackIdVec.size());
716 // mDaqOccupancyMatch2->Fill((tempVec[ii].module-1)*mNCell+(tempVec[ii].cell-1));
717  Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
718  Float_t dy = tempVec[ii].yhit-ycenter;
719  Float_t dz = tempVec[ii].zhit;
720  mDeltaHitMatch2->Fill(dy, dz);
721  }
722  }
723 
724  erasedVec = tempVec;
725  while (tempVec.size() != 0) {
726  StructCellHit cellHit;
727  Int_t nCells = 0;
728  idVector vTrackId;
729  vector<StThreeVectorD> vPosition;
730  vector<Int_t> vchannel, vtray, vmodule, vcell;
731  vector<Float_t> vzhit, vyhit;
732  vector<Double_t> vtot, vtheta;
733  vector<Int_t> vindex2BTofHit;
734 
735  tofCellHitVectorIter tempIter=tempVec.begin();
736  tofCellHitVectorIter erasedIter=erasedVec.begin();
737  while(erasedIter!= erasedVec.end()) {
738  if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
739  nCells++;
740  vtray.push_back(erasedIter->tray);
741  vmodule.push_back(erasedIter->module);
742  vcell.push_back(erasedIter->cell);
743  vPosition.push_back(erasedIter->hitPosition);
744  vTrackId.push_back(erasedIter->trackIdVec.back());
745  vzhit.push_back(erasedIter->zhit);
746  vyhit.push_back(erasedIter->yhit);
747  vtot.push_back(erasedIter->tot);
748  vindex2BTofHit.push_back(erasedIter->index2BTofHit);
749  vtheta.push_back(erasedIter->theta);
750 
751  erasedVec.erase(erasedIter);
752  erasedIter--;
753  }
754  erasedIter++;
755  }
756 
757  if (nCells==1){
758  // for singly hit cell, copy data in singleHitCellsVec
759  cellHit.tray = vtray[0];
760  cellHit.module = vmodule[0];
761  cellHit.cell = vcell[0];
762  cellHit.trackIdVec.push_back(vTrackId[0]);
763  cellHit.hitPosition = vPosition[0];
764  cellHit.matchFlag = 1;
765  cellHit.zhit = vzhit[0];
766  cellHit.yhit = vyhit[0];
767  cellHit.tot = vtot[0];
768  cellHit.index2BTofHit = vindex2BTofHit[0];
769  cellHit.theta = vtheta[0];
770 
771  FinalMatchedCellsVec.push_back(cellHit);
772 
773  // debugging output
774  if(Debug()) {
775  LOG_DEBUG << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
776  idVectorIter ij=vTrackId.begin();
777  while (ij != vTrackId.end()) { LOG_DEBUG << " " << *ij; ij++; }
778  LOG_DEBUG << endm;
779  }
780 
781  }
782  else if (nCells>1){ // for multiple hit cells find the most likely candidate.
783  Int_t thiscandidate(-99);
784  Int_t thisMatchFlag(0);
785 
786  // sort on tot
787  Float_t tot(0.);
788  vector<Int_t> ttCandidates;
789  for (Int_t i=0;i<nCells;i++) {
790  Double_t tt = vtot[i];
791  if(tt<40.&&tt>tot) { // open the ToT cut to 40 ns
792  tot = tt;
793  ttCandidates.clear();
794  ttCandidates.push_back(i);
795  } else if (tt==tot) {
796  ttCandidates.push_back(i);
797  }
798  }
799  if (ttCandidates.size()==1) {
800  thiscandidate = ttCandidates[0];
801  thisMatchFlag = 2;
802  } else if (ttCandidates.size()>1) { // sort on hitposition
803  Float_t ss(99.);
804  vector<Int_t> ssCandidates;
805  for(size_t j=0;j<ttCandidates.size();j++) {
806  Float_t yy = vyhit[ttCandidates[j]];
807  Float_t ycell = (vcell[ttCandidates[j]]-1-2.5)*mWidthPad;
808  Float_t ll = fabs(yy-ycell);
809  if(ll<ss) {
810  ss = ll;
811  ssCandidates.clear();
812  ssCandidates.push_back(ttCandidates[j]);
813  }else if (ll==ss)
814  ssCandidates.push_back(ttCandidates[j]);
815  }
816  if (ssCandidates.size()==1){
817  thiscandidate = ssCandidates[0];
818  thisMatchFlag = 3;
819  }
820  }
821 
822  if (thiscandidate>=0) {
823  cellHit.tray = vtray[thiscandidate];
824  cellHit.module = vmodule[thiscandidate];
825  cellHit.cell = vcell[thiscandidate];
826  cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
827  cellHit.hitPosition = vPosition[thiscandidate];
828  cellHit.matchFlag = thisMatchFlag;
829  cellHit.zhit = vzhit[thiscandidate];
830  cellHit.yhit = vyhit[thiscandidate];
831  cellHit.tot = vtot[thiscandidate];
832  cellHit.index2BTofHit = vindex2BTofHit[thiscandidate];
833  cellHit.theta = vtheta[thiscandidate];
834 
835  FinalMatchedCellsVec.push_back(cellHit);
836 
837  // debugging output
838  if(Debug()) { LOG_DEBUG << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:" << vTrackId[thiscandidate] << endm; }
839  }
840 
841  } else {
842  LOG_WARN << "E: no cells belong to this track ... should not happen!" << endm;
843  }
844 
845  tempVec = erasedVec;
846  }
847 
848  if(Debug()) { LOG_INFO << "E: before/after: " << singleHitCellsVec.size() << "/" << FinalMatchedCellsVec.size() << endm; }
849  // end of Sect.E
850  if (doPrintCpuInfo) {
851  timer.stop();
852  LOG_INFO << "CPU time after Step E - sorting : "
853  << timer.elapsedTime() << " sec" << endm;
854  timer.start();
855  }
856 
857  //.........................................................................
859  //
860  if(mHisto) {
861  if(FinalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
862  mCellsPerEventMatch3->Fill(FinalMatchedCellsVec.size());
863  }
864 
865 // StSPtrVecBTofHit& tofHits = theTof->tofHits();
866  Int_t nValidSingleHitCells(0), nValidSinglePrimHitCells(0);
867 
868  for (size_t ii=0; ii < FinalMatchedCellsVec.size(); ii++){
869  Int_t tray = FinalMatchedCellsVec[ii].tray;
870  Int_t module = FinalMatchedCellsVec[ii].module;
871  Int_t cell = FinalMatchedCellsVec[ii].cell;
872 
873  Float_t ycenter = (cell-1-2.5)*mWidthPad;
874  Float_t dy = FinalMatchedCellsVec[ii].yhit - ycenter;
875  Float_t dz = FinalMatchedCellsVec[ii].zhit;
876  if (FinalMatchedCellsVec[ii].trackIdVec.size()!=1)
877  LOG_WARN << "F: WHAT!?! mult.matched cell in single cell list " << tray << " " << module << " " << cell << endm;
878 
879  if(mHisto) {
880  mTracksPerCellMatch3->Fill(FinalMatchedCellsVec[ii].trackIdVec.size());
881 // mDaqOccupancyMatch3->Fill((module-1)*mNCell+(cell-1));
882  mDeltaHitMatch3->Fill(dy, dz);
883  mDeltaHitFinal[tray-1]->Fill(dy,dz);
884  }
885 
886  // get track-id from cell hit vector
887  int trackNode = FinalMatchedCellsVec[ii].trackIdVec[0];
888  StGlobalTrack *globalTrack = dynamic_cast<StGlobalTrack*>(nodes[trackNode]->track(global));
889  if(!globalTrack) {
890  LOG_WARN << "Wrong global track!" << endm;
891  continue;
892  }
893 
894  // Fill association in TOF Hit Collection
895  //fg StBTofHit *tofHit = dynamic_cast<StBTofHit *>(tofHits[FinalMatchedCellsVec[ii].index2BTofHit]);
896  StBTofHit *tofHit = tofHits[FinalMatchedCellsVec[ii].index2BTofHit];
897  if(tofHit->tray()!=tray || tofHit->module()!=module || tofHit->cell()!=cell) {
898  LOG_WARN << "Wrong hit in the BTofHitCollection!" << endm;
899  continue;
900  }
901  nValidSingleHitCells++;
902 
904  tofHit->setAssociatedTrack(globalTrack);
905 
906  // Fill the matched data in StBTofPidTraits
907  StBTofPidTraits *pidTof = new StBTofPidTraits();
908  pidTof->setTofHit(tofHit);
909  pidTof->setMatchFlag(FinalMatchedCellsVec[ii].matchFlag);
910  pidTof->setYLocal(dy);
911  pidTof->setZLocal(FinalMatchedCellsVec[ii].zhit);
912  pidTof->setPosition(FinalMatchedCellsVec[ii].hitPosition);
913  pidTof->setThetaLocal(FinalMatchedCellsVec[ii].theta);
914  globalTrack->addPidTraits(pidTof);
915 
916  StPrimaryTrack *pTrack = dynamic_cast<StPrimaryTrack*>(nodes[trackNode]->track(primary));
917  if(pTrack) {
918  nValidSinglePrimHitCells++;
919  StBTofPidTraits *ppidTof = new StBTofPidTraits();
920  ppidTof->setTofHit(tofHit);
921  ppidTof->setMatchFlag(FinalMatchedCellsVec[ii].matchFlag);
922  ppidTof->setYLocal(dy);
923  ppidTof->setZLocal(FinalMatchedCellsVec[ii].zhit);
924  ppidTof->setPosition(FinalMatchedCellsVec[ii].hitPosition);
925  ppidTof->setThetaLocal(FinalMatchedCellsVec[ii].theta);
926  pTrack->addPidTraits(ppidTof);
927  }
928 
929  } // end final matched cells
930 
931  if(mHisto) {
932  mCellsPrimaryPerEventMatch3->Fill(nValidSinglePrimHitCells);
933  }
934 
935  if(Debug()) { LOG_INFO << "F: before/after" << FinalMatchedCellsVec.size() << "/" <<nValidSinglePrimHitCells << endm; }
936  // end of Sect.F
937  if (doPrintCpuInfo) {
938  timer.stop();
939  LOG_INFO << "CPU time after Step F - final : "
940  << timer.elapsedTime() << " sec" << endm;
941  timer.start();
942  }
943 
944  LOG_INFO << "#(daq): " << daqCellsHitVec.size()
945  << " #(proj): " << allCellsHitVec.size()
946  << " #(prim proj): " << nPrimaryHits
947  << " #(matched): " << FinalMatchedCellsVec.size()
948  << " #(single valid): " << nValidSingleHitCells
949  << " #(single prim valid): " << nValidSinglePrimHitCells
950  << endm;
951 
952  //check StEvent collections --
953  if (theTof->hitsPresent()){
954  LOG_DEBUG << " BTofCollection: hit container present."<<endm;
955  if (Debug()){
956  StSPtrVecBTofHit& tmpCellTofVec = theTof->tofHits();
957  LOG_INFO << " # of hits in this event:" << tmpCellTofVec.size() << endm;
958  for (size_t i = 0; i < tmpCellTofVec.size(); i++) {
959  StBTofHit* p = tmpCellTofVec[i];
960  LOG_INFO << (*p) << endm;
961  }
962  }
963  }
964  //-- end check
965 
966 
967  if (doPrintMemoryInfo) {
968  StMemoryInfo::instance()->snapshot();
969  StMemoryInfo::instance()->print();
970  }
971  if (doPrintCpuInfo) {
972  timer.stop();
973  LOG_INFO << "CPU time for StBTofMatchMaker::Make(): "
974  << timer.elapsedTime() << " sec" << endm;
975  }
976 
977  LOG_DEBUG << "StBTofMatchMaker -- bye-bye" << endm;
978 
979  return;
980 }
981 
982 //---------------------------------------------------------------------------
983 void StBTofMatchMaker::processMuDst(){
984 
985  if(mHisto) mEventCounterHisto->Fill(0);
986  // event selection ...
987  StMuDstMaker *mMuDstMaker = (StMuDstMaker *)GetMaker("MuDst");
988  if(!mMuDstMaker) {
989  LOG_WARN << " No MuDstMaker ... bye-bye ... " << endm;
990  return;
991  }
992  mMuDst = mMuDstMaker->muDst();
993  if(!mMuDst) {
994  LOG_WARN << " No MuDst ... bye-bye" << endm;
995  return;
996  }
997 
998  // timing & memory info -only when requested-
999  StTimer timer;
1000  if (doPrintCpuInfo) timer.start();
1001  if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
1002 
1003  //.........................................................................
1004  // read data from StBTofHit
1006  //
1007  tofCellHitVector daqCellsHitVec;
1008  idVector validModuleVec;
1009 
1010  // multi-tray system
1011  Int_t nhits = mMuDst->numberOfBTofHit();
1012  LOG_INFO << " Number of BTOF Hits = " << nhits << endm;
1013  if(mHisto&&nhits>0) mEventCounterHisto->Fill(1);
1014  for(int i=0;i<nhits;i++) {
1015  StMuBTofHit *aHit = (StMuBTofHit*)mMuDst->btofHit(i);
1016  if(!aHit) continue;
1017  if(aHit->tray()<=0||aHit->tray()>mNTray) continue; // barrel tray hits
1018 
1019  // clean up any association done before
1020  aHit->setIndex2Primary(-1);
1021  aHit->setIndex2Global(-1);
1022  aHit->setAssociatedTrackId(-1);
1023 
1024  int trayId = aHit->tray();
1025  int moduleId = aHit->module();
1026  int cellId = aHit->cell();
1027 
1028  if(Debug()) { LOG_INFO <<"A: fired hit in " << " tray="<< trayId << " module="<<moduleId<<" cell="<<cellId<<endm; }
1029 
1030  StructCellHit aDaqCellHit;
1031  aDaqCellHit.tray = trayId;
1032  aDaqCellHit.module = moduleId;
1033  aDaqCellHit.cell = cellId;
1034  aDaqCellHit.tot = aHit->tot();
1035  aDaqCellHit.index2BTofHit = i;
1036  daqCellsHitVec.push_back(aDaqCellHit);
1037 
1038  // additional valid number configuration
1039  int id = trayId*100+moduleId;
1040  //fg bool ifind = kFALSE;
1041  //fg for(size_t im=0;im<validModuleVec.size();im++) {
1042  //fg if(id==validModuleVec[im]) {
1043  //fg ifind = kTRUE;
1044  //fg break;
1045  //fg }
1046  //fg }
1047  //fg if(!ifind) validModuleVec.push_back(id);
1048  if (find(validModuleVec.begin(), validModuleVec.end(), id) == validModuleVec.end())
1049  validModuleVec.push_back(id);
1050 
1051  if(mHisto) {
1052  mDaqOccupancy[trayId-1]->Fill((moduleId-1)*mNCell+(cellId-1));
1053  }
1054  }
1055 
1056  // end of Sect.A
1057  if(Debug()) {
1058  LOG_DEBUG << " total # of cells = " << daqCellsHitVec.size() << endm;
1059  for(size_t iv = 0;iv<validModuleVec.size();iv++) {
1060  LOG_DEBUG << " module # " << validModuleVec[iv] << " Valid! " << endm;
1061  }
1062  }
1063  if(mHisto) {
1064  mCellsMultInEvent->Fill(daqCellsHitVec.size());
1065  if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
1066  }
1067 // if(!daqCellsHitVec.size()) return;
1068  if (doPrintCpuInfo) {
1069  timer.stop();
1070  LOG_INFO << "CPU time after Step A - loading hits : "
1071  << timer.elapsedTime() << " sec" << endm;
1072  timer.start();
1073  }
1074 
1075 
1076  //.........................................................................
1078  //
1079  Int_t index2Primary[50000]; // map
1080  memset(index2Primary, -1, sizeof(index2Primary));
1081  for(int ip=0;ip<(int)mMuDst->array(muPrimary)->GetEntries();ip++) {
1082  StMuTrack *pTrack = (StMuTrack *)mMuDst->array(muPrimary)->UncheckedAt(ip);
1083  if(!pTrack) continue;
1084  int gIndex = pTrack->index2Global();
1085  if(gIndex<0) continue;
1086  index2Primary[gIndex] = ip;
1087  }
1088 
1089  tofCellHitVector allCellsHitVec;
1090  StructCellHit cellHit;
1091 
1092  Int_t nGlobals = mMuDst->numberOfGlobalTracks();
1093  Int_t nAllTracks=0;
1094  Int_t nPrimaryHits = 0;
1095  for (int iNode=0; iNode<nGlobals; iNode++){
1096  tofCellHitVector cellHitVec;
1097  // cellHitVec.clear();
1098  StMuTrack *theTrack = mMuDst->globalTracks(iNode);
1099  if(!theTrack) continue;
1100 
1101  // clean up any association done before
1102  StMuBTofPidTraits pidTof;
1103  theTrack->setBTofPidTraits(pidTof);
1104  theTrack->setIndex2BTofHit(-1);
1105 
1106  bool isPrimary = kFALSE;
1107  int pIndex = index2Primary[iNode];
1108  if(pIndex>=0) {
1109  isPrimary = kTRUE;
1110  StMuTrack *pTrack = (StMuTrack *)mMuDst->array(muPrimary)->UncheckedAt(pIndex);
1111  if(pTrack) {
1112  pTrack->setBTofPidTraits(pidTof);
1113  pTrack->setIndex2BTofHit(-1);
1114  }
1115  }
1116 
1117  StThreeVectorF mom = theTrack->momentum();
1118  float pt = mom.perp();
1119  float eta = mom.pseudoRapidity();
1120  float phi = mom.phi();
1121  //fg if (phi<0.) phi += 2.*3.14159;
1122  if (phi<0.) phi += 2.*M_PI;
1123 
1124  float nSigmaPion = theTrack->nSigmaPion();
1125  float dEdx = theTrack->dEdx();
1126  int ndEdxpts = theTrack->nHitsDedx();
1127  int nfitpts = theTrack->nHitsFit(kTpcId);
1128 
1129  // make sure we have a track, a miniDST might have removed it...
1130  if (validTrack(theTrack)){
1131  if(mHisto) {
1132  mTrackPtEta->Fill(pt, eta);
1133  mTrackPtPhi->Fill(pt, phi);
1134  mTrackNFitPts->Fill(nfitpts);
1135  if(dEdx>0.) mTrackdEdxvsp->Fill(mom.mag(), dEdx*1.e6);
1136  if(fabs(nSigmaPion)<5.) mNSigmaPivsPt->Fill(pt, nSigmaPion+5.*theTrack->charge());
1137  }
1138 
1139  if(mSaveTree) {
1140  trackTree.pt = pt;
1141  trackTree.eta = eta;
1142  trackTree.phi = phi;
1143  trackTree.nfitpts = nfitpts;
1144  trackTree.dEdx = dEdx*1.e6;
1145  trackTree.ndEdxpts = ndEdxpts;
1146  trackTree.charge = theTrack->charge();
1147  trackTree.projTrayId = 0;
1148  trackTree.projCellChan = -1;
1149  trackTree.projY = -999.;
1150  trackTree.projZ = -999.;
1151  }
1152 
1153  nAllTracks++;
1154  StPhysicalHelixD theHelix = mOuterTrackGeometry ? theTrack->outerHelix() : theTrack->helix();
1155 
1156 // IntVec projTrayVec;
1157 // if(!mBTofGeom->projTrayVector(theHelix, projTrayVec)) continue;
1158 
1159  IntVec idVec;
1160  DoubleVec pathVec, thetaVec;
1161  PointVec crossVec;
1162 
1163 // idVec.clear();
1164 // pathVec.clear();
1165 // crossVec.clear();
1166 
1167  Int_t ncells = 0;
1168  if(mBTofGeom->HelixCrossCellIds(theHelix,idVec,pathVec,crossVec,thetaVec) ) {
1169 // if(mBTofGeom->HelixCrossCellIds(theHelix, validModuleVec, projTrayVec, idVec, pathVec, crossVec)) {
1170  Int_t cells = idVec.size();
1171  for (Int_t i=0; i<cells; i++) {
1172  Int_t icell,imodule,itray;
1173  Double_t local[3],global[3];
1174  for(Int_t i2=0;i2<3;i2++){
1175  local[i2]=0;
1176  }
1177  global[0]=crossVec[i].x();
1178  global[1]=crossVec[i].y();
1179  global[2]=crossVec[i].z();
1180  mBTofGeom->DecodeCellId(idVec[i], icell, imodule, itray);
1181 // LOG_INFO << " decode " << idVec[i] << " to tray#" << itray << " module#" << imodule << " cell#" << icell << endm;
1182  StBTofGeomSensor* sensor = mBTofGeom->GetGeomSensor(imodule,itray);
1183  if(!sensor) {
1184  LOG_WARN << " No sensitive module in the projection??? -- Something weird!!! " << endm;
1185  continue;
1186  }
1187  sensor->Master2Local(&global[0],&local[0]);
1188  icell = sensor->FindCellIndex(local);
1189  // StThreeVectorD glo=sensor->GetCenterPosition();
1190  StThreeVectorD glo(global[0], global[1], global[2]);
1191  StThreeVectorD hitPos(local[0], local[1], local[2]);
1192 // delete sensor; /// function in StBTofGeometry modified. Don't delete here.
1193  if (fabs(local[2])<mZLocalCut) { // to be consistent with GEANT geometry
1194  // and the alignment calibration
1195 // if (local[2]<=2.7&&local[2]>=-3.4) {
1196 // if (fabs(local[2])<4.) { // loose local z cut at first step
1197  ncells++;
1198  cellHit.tray = itray;
1199  cellHit.module = imodule;
1200  cellHit.cell = icell;
1201  cellHit.trackIdVec.push_back(iNode);
1202  cellHit.hitPosition = glo; // global position
1203  cellHit.zhit = (Float_t)hitPos.z();
1204  cellHit.yhit = (Float_t)hitPos.y();
1205  cellHit.theta = (Double_t)thetaVec[i];
1206  cellHitVec.push_back(cellHit);
1207  allCellsHitVec.push_back(cellHit);
1208 
1209  if(isPrimary) nPrimaryHits++;
1210 
1211  if(mHisto) {
1212  mDaqOccupancyProj[itray-1]->Fill((imodule-1)*mNCell+(icell-1));
1213  mHitsPosition->Fill(hitPos.y(), hitPos.z());
1214  }
1215 
1216  if(mSaveTree) {
1217  trackTree.projTrayId = itray;
1218  trackTree.projCellChan = (imodule-1)*mNCell+(icell-1);
1219  trackTree.projY = local[1];
1220  trackTree.projZ = local[2];
1221  }
1222 
1223  if(Debug()) {
1224  LOG_DEBUG <<"B: nodeid=" << iNode << " projected in " << " tray="<< itray << " module="<<imodule<<" cell="<<icell<<endm;
1225  LOG_DEBUG <<" hit position " << hitPos << endm;
1226  // LOG_DEBUG <<" momemtum= " << pt << " " << eta << " " << phi << endm;
1227  }
1228  }
1229  } // for (Int_t i=0...)
1230  } // endif(helixcross...)
1231  if(ncells>0&&mHisto) mHitsMultPerTrack->Fill(ncells);
1232 
1233  if(mHisto && mSaveTree) mTrackTree->Fill();
1234 
1235  } // if(ValidTrack)..
1236  } // loop over nodes
1237  if(Debug()) { LOG_INFO << "B: matched/available/total #tracknodes: " <<allCellsHitVec.size() << "/" <<nAllTracks << "/" << nGlobals << endm; }
1238  if(mHisto) {
1239  mHitsMultInEvent->Fill(allCellsHitVec.size());
1240  mHitsPrimaryInEvent->Fill(nPrimaryHits);
1241  if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
1242  }
1243  // end of Sect.B
1244  if (doPrintCpuInfo) {
1245  timer.stop();
1246  LOG_INFO << "CPU time after Step B - projection : "
1247  << timer.elapsedTime() << " sec" << endm;
1248  timer.start();
1249  }
1250 
1251  //.........................................................................
1253  //
1254  tofCellHitVector matchHitCellsVec;
1255 
1256  tofCellHitVectorIter daqIter = daqCellsHitVec.begin();
1257  for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
1258  tofCellHitVectorIter proIter = allCellsHitVec.begin();
1259  for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
1260 
1261  int daqIndex = (daqIter->module-1)*6 + (daqIter->cell-1);
1262  int proIndex = (proIter->module-1)*6 + (proIter->cell-1);
1263  int hisIndex = daqIter->tray - 1;
1264 // int daqAllIndex = (daqIter->tray - 1)*mNTOF + daqIndex;
1265 // int proAllIndex = (proIter->tray - 1)*mNTOF + proIndex;
1266  if(daqIter->tray==proIter->tray) {
1267  if (mHisto) {
1268  if(hisIndex>=0&&hisIndex<mNTray) {
1269  mHitCorr[hisIndex]->Fill(proIndex,daqIndex);
1270  mHitCorrModule[hisIndex]->Fill(proIter->module-1,daqIter->module-1);
1271  } else {
1272  LOG_WARN << " weird tray # " << daqIter->tray << endm;
1273  }
1274  }
1275  }
1276  if( (daqIter->tray==proIter->tray)&&
1277  (daqIter->module==proIter->module) &&
1278  ( ( (proIter->cell==6)&&((proIter->cell==daqIter->cell) ||
1279  (proIter->cell==daqIter->cell+1)) )
1280  || ( (proIter->cell==1)&&((proIter->cell==daqIter->cell) ||
1281  (proIter->cell==daqIter->cell-1)) )
1282  || ( (proIter->cell>=2&&proIter->cell<=6) &&
1283  ( (proIter->cell==daqIter->cell) ||
1284  (proIter->cell==daqIter->cell-1) ||
1285  (proIter->cell==daqIter->cell+1) ) ) ) ) {
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  cellHit.tot = daqIter->tot;
1294  cellHit.index2BTofHit = daqIter->index2BTofHit;
1295  cellHit.theta = proIter->theta;
1296  matchHitCellsVec.push_back(cellHit);
1297  }
1298  }
1299  } //end {sec. C}
1300  if(Debug()) { LOG_INFO << "C: before/after: " << allCellsHitVec.size() << "/" << matchHitCellsVec.size() << endm; }
1301  if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
1302  if (doPrintCpuInfo) {
1303  timer.stop();
1304  LOG_INFO << "CPU time after Step C - matching : "
1305  << timer.elapsedTime() << " sec" << endm;
1306  timer.start();
1307  }
1308 
1309  //.........................................................................
1311  //
1312  Int_t nSingleHitCells(0);
1313  Int_t nMultiHitsCells(0);
1314 
1315  tofCellHitVector singleHitCellsVec;
1316  tofCellHitVector multiHitsCellsVec;
1317 
1318  tofCellHitVector tempVec = matchHitCellsVec;
1319  tofCellHitVector erasedVec = tempVec;
1320  while (tempVec.size() != 0) {
1321  Int_t nTracks = 0;
1322  idVector trackIdVec;
1323 
1324  tofCellHitVectorIter tempIter=tempVec.begin();
1325  tofCellHitVectorIter erasedIter=erasedVec.begin();
1326  while(erasedIter!= erasedVec.end()) {
1327  if(tempIter->tray == erasedIter->tray &&
1328  tempIter->module == erasedIter->module &&
1329  tempIter->cell == erasedIter->cell) {
1330  nTracks++;
1331  trackIdVec.push_back(erasedIter->trackIdVec.back()); // merge
1332  erasedVec.erase(erasedIter);
1333  erasedIter--;
1334  }
1335  erasedIter++;
1336  }
1337 
1338  cellHit.cell = tempIter->cell;
1339  cellHit.module = tempIter->module;
1340  cellHit.tray = tempIter->tray;
1341  cellHit.hitPosition = tempIter->hitPosition;
1342  cellHit.trackIdVec = trackIdVec;
1343  cellHit.zhit = tempIter->zhit;
1344  cellHit.yhit = tempIter->yhit;
1345  cellHit.tot = tempIter->tot;
1346  cellHit.index2BTofHit = tempIter->index2BTofHit;
1347  cellHit.theta = tempIter->theta;
1348 
1349  Float_t ycenter = (tempIter->cell-1-2.5)*mWidthPad;
1350  Float_t dy = tempIter->yhit - ycenter;
1351  Float_t dz = tempIter->zhit;
1352 
1353  if(mHisto) {
1354  mTracksPerCellMatch1->Fill(trackIdVec.size());
1355 // mDaqOccupancyMatch1->Fill((tempIter->module-1)*mNCell+(tempIter->cell-1));
1356  mDeltaHitMatch1->Fill(dy, dz);
1357  }
1358 
1359  if (nTracks==1){
1360  nSingleHitCells++;
1361  singleHitCellsVec.push_back(cellHit);
1362  } else if (nTracks>1){
1363  nMultiHitsCells++;
1364  multiHitsCellsVec.push_back(cellHit);
1365  // for multiple hit cells either discard (yes) or
1366  // find the most likely candidate.
1367  } else {
1368  LOG_WARN << "D: no tracks extrapolate to matched cell ... should not happen!" << endm;
1369  }
1370 
1371  if(Debug()) {
1372  LOG_DEBUG << "D: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
1373  idVectorIter ij=trackIdVec.begin();
1374  while (ij != trackIdVec.end()) { LOG_DEBUG << " " << *ij; ij++; }
1375  LOG_DEBUG << endm;
1376  }
1377 
1378  tempVec = erasedVec;
1379  }
1380  if(Debug()) { LOG_INFO << "D: before/after: " << matchHitCellsVec.size() << "/" << singleHitCellsVec.size() << endm; }
1381  //end of Sect.D
1382 
1383  if(mHisto) {
1384  mCellsPerEventMatch1->Fill(singleHitCellsVec.size()+multiHitsCellsVec.size());
1385  if(singleHitCellsVec.size()) mEventCounterHisto->Fill(9);
1386  }
1387  if (doPrintCpuInfo) {
1388  timer.stop();
1389  LOG_INFO << "CPU time after Step D - erasing : "
1390  << timer.elapsedTime() << " sec" << endm;
1391  timer.start();
1392  }
1393 
1394  //.........................................................................
1396  //
1397  tofCellHitVector FinalMatchedCellsVec;
1398  // FinalMatchedCellsVec.clear();
1399  tempVec = singleHitCellsVec;
1400  if(mHisto) {
1401  mCellsPerEventMatch2->Fill(tempVec.size());
1402  for(unsigned int ii=0;ii<tempVec.size();ii++) {
1403  mTracksPerCellMatch2->Fill(tempVec[ii].trackIdVec.size());
1404 // mDaqOccupancyMatch2->Fill((tempVec[ii].module-1)*mNCell+(tempVec[ii].cell-1));
1405  Float_t ycenter = (tempVec[ii].cell-1-2.5)*mWidthPad;
1406  Float_t dy = tempVec[ii].yhit-ycenter;
1407  Float_t dz = tempVec[ii].zhit;
1408  mDeltaHitMatch2->Fill(dy, dz);
1409  }
1410  }
1411 
1412  erasedVec = tempVec;
1413  while (tempVec.size() != 0) {
1414  StructCellHit cellHit;
1415  Int_t nCells = 0;
1416  idVector vTrackId;
1417  vector<StThreeVectorD> vPosition;
1418  vector<Int_t> vchannel, vtray, vmodule, vcell;
1419  vector<Float_t> vzhit, vyhit;
1420  vector<Double_t> vtot, vtheta;
1421  vector<Int_t> vindex2BTofHit;
1422 
1423  tofCellHitVectorIter tempIter=tempVec.begin();
1424  tofCellHitVectorIter erasedIter=erasedVec.begin();
1425  while(erasedIter!= erasedVec.end()) {
1426  if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
1427  nCells++;
1428  vtray.push_back(erasedIter->tray);
1429  vmodule.push_back(erasedIter->module);
1430  vcell.push_back(erasedIter->cell);
1431  vPosition.push_back(erasedIter->hitPosition);
1432  vTrackId.push_back(erasedIter->trackIdVec.back());
1433  vzhit.push_back(erasedIter->zhit);
1434  vyhit.push_back(erasedIter->yhit);
1435  vtot.push_back(erasedIter->tot);
1436  vindex2BTofHit.push_back(erasedIter->index2BTofHit);
1437  vtheta.push_back(erasedIter->theta);
1438 
1439  erasedVec.erase(erasedIter);
1440  erasedIter--;
1441  }
1442  erasedIter++;
1443  }
1444 
1445  if (nCells==1){
1446  // for singly hit cell, copy data in singleHitCellsVec
1447  cellHit.tray = vtray[0];
1448  cellHit.module = vmodule[0];
1449  cellHit.cell = vcell[0];
1450  cellHit.trackIdVec.push_back(vTrackId[0]);
1451  cellHit.hitPosition = vPosition[0];
1452  cellHit.matchFlag = 1;
1453  cellHit.zhit = vzhit[0];
1454  cellHit.yhit = vyhit[0];
1455  cellHit.tot = vtot[0];
1456  cellHit.index2BTofHit = vindex2BTofHit[0];
1457  cellHit.theta = vtheta[0];
1458 
1459  FinalMatchedCellsVec.push_back(cellHit);
1460 
1461  // debugging output
1462  if(Debug()) {
1463  LOG_DEBUG << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
1464  idVectorIter ij=vTrackId.begin();
1465  while (ij != vTrackId.end()) { LOG_DEBUG << " " << *ij; ij++; }
1466  LOG_DEBUG << endm;
1467  }
1468 
1469  }
1470  else if (nCells>1){ // for multiple hit cells find the most likely candidate.
1471  Int_t thiscandidate(-99);
1472  Int_t thisMatchFlag(0);
1473 
1474  // sort on tot
1475  Float_t tot(0.);
1476  vector<Int_t> ttCandidates;
1477  for (Int_t i=0;i<nCells;i++) {
1478  Double_t tt = vtot[i];
1479  if(tt<40.&&tt>tot) { // open the ToT to 40 ns
1480  tot = tt;
1481  ttCandidates.clear();
1482  ttCandidates.push_back(i);
1483  } else if (tt==tot) {
1484  ttCandidates.push_back(i);
1485  }
1486  }
1487  if (ttCandidates.size()==1) {
1488  thiscandidate = ttCandidates[0];
1489  thisMatchFlag = 2;
1490  } else if (ttCandidates.size()>1) { // sort on hitposition
1491  Float_t ss(99.);
1492  vector<Int_t> ssCandidates;
1493  for(size_t j=0;j<ttCandidates.size();j++) {
1494  Float_t yy = vyhit[ttCandidates[j]];
1495  Float_t ycell = (vcell[ttCandidates[j]]-1-2.5)*mWidthPad;
1496  Float_t ll = fabs(yy-ycell);
1497  if(ll<ss) {
1498  ss = ll;
1499  ssCandidates.clear();
1500  ssCandidates.push_back(ttCandidates[j]);
1501  }else if (ll==ss)
1502  ssCandidates.push_back(ttCandidates[j]);
1503  }
1504  if (ssCandidates.size()==1){
1505  thiscandidate = ssCandidates[0];
1506  thisMatchFlag = 3;
1507  }
1508  }
1509 
1510  if (thiscandidate>=0) {
1511  cellHit.tray = vtray[thiscandidate];
1512  cellHit.module = vmodule[thiscandidate];
1513  cellHit.cell = vcell[thiscandidate];
1514  cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
1515  cellHit.hitPosition = vPosition[thiscandidate];
1516  cellHit.matchFlag = thisMatchFlag;
1517  cellHit.zhit = vzhit[thiscandidate];
1518  cellHit.yhit = vyhit[thiscandidate];
1519  cellHit.tot = vtot[thiscandidate];
1520  cellHit.index2BTofHit = vindex2BTofHit[thiscandidate];
1521  cellHit.theta = vtheta[thiscandidate];
1522 
1523  FinalMatchedCellsVec.push_back(cellHit);
1524 
1525  // debugging output
1526  if(Debug()) { LOG_DEBUG << "E: itray=" << cellHit.tray << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:" << vTrackId[thiscandidate] << endm; }
1527  }
1528 
1529  } else {
1530  LOG_WARN << "E: no cells belong to this track ... should not happen!" << endm;
1531  }
1532 
1533  tempVec = erasedVec;
1534  }
1535 
1536  if(Debug()) { LOG_INFO << "E: before/after: " << singleHitCellsVec.size() << "/" << FinalMatchedCellsVec.size() << endm; }
1537  // end of Sect.E
1538  if (doPrintCpuInfo) {
1539  timer.stop();
1540  LOG_INFO << "CPU time after Step E - sorting : "
1541  << timer.elapsedTime() << " sec" << endm;
1542  timer.start();
1543  }
1544 
1545  //.........................................................................
1547  //
1548  if(mHisto) {
1549  if(FinalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
1550  mCellsPerEventMatch3->Fill(FinalMatchedCellsVec.size());
1551  }
1552 
1553 // StSPtrVecBTofHit& tofHits = theTof->tofHits();
1554  Int_t nValidSingleHitCells(0), nValidSinglePrimHitCells(0);
1555 
1556  for (size_t ii=0; ii < FinalMatchedCellsVec.size(); ii++){
1557  Int_t tray = FinalMatchedCellsVec[ii].tray;
1558  Int_t module = FinalMatchedCellsVec[ii].module;
1559  Int_t cell = FinalMatchedCellsVec[ii].cell;
1560 
1561  Float_t ycenter = (cell-1-2.5)*mWidthPad;
1562  Float_t dy = FinalMatchedCellsVec[ii].yhit - ycenter;
1563  Float_t dz = FinalMatchedCellsVec[ii].zhit;
1564  if (FinalMatchedCellsVec[ii].trackIdVec.size()!=1)
1565  LOG_WARN << "F: WHAT!?! mult.matched cell in single cell list " << tray << " " << module << " " << cell << endm;
1566 
1567  if(mHisto) {
1568  mTracksPerCellMatch3->Fill(FinalMatchedCellsVec[ii].trackIdVec.size());
1569 // mDaqOccupancyMatch3->Fill((module-1)*mNCell+(cell-1));
1570  mDeltaHitMatch3->Fill(dy, dz);
1571  mDeltaHitFinal[tray-1]->Fill(dy,dz);
1572  }
1573 
1574  // get track-id from cell hit vector
1575  int trackNode = FinalMatchedCellsVec[ii].trackIdVec[0];
1576  StMuTrack *globalTrack = mMuDst->globalTracks(trackNode);
1577  if(!globalTrack) {
1578  LOG_WARN << "Wrong global track!" << endm;
1579  continue;
1580  }
1581 
1582  // Fill association in TOF Hit Collection
1583  StMuBTofHit *tofHit = mMuDst->btofHit(FinalMatchedCellsVec[ii].index2BTofHit);
1584  if(tofHit->tray()!=tray || tofHit->module()!=module || tofHit->cell()!=cell) {
1585  LOG_WARN << "Wrong hit in the MuBTofHit!" << endm;
1586  continue;
1587  }
1588  nValidSingleHitCells++;
1589 
1591  tofHit->setAssociatedTrackId(globalTrack->id());
1592  tofHit->setIndex2Global(trackNode);
1593  globalTrack->setIndex2BTofHit(FinalMatchedCellsVec[ii].index2BTofHit);
1594 
1595  int ip = index2Primary[trackNode];
1596  StMuTrack *primaryTrack = 0;
1597  if(ip>=0) {
1598  nValidSinglePrimHitCells++;
1599  tofHit->setIndex2Primary(ip);
1600  primaryTrack = (StMuTrack *)mMuDst->array(muPrimary)->UncheckedAt(ip);
1601  if(primaryTrack) {
1602  primaryTrack->setIndex2BTofHit(FinalMatchedCellsVec[ii].index2BTofHit);
1603  }
1604  }
1605 
1606  // Fill the matched data in StBTofPidTraits
1607  StMuBTofPidTraits pidTof = globalTrack->btofPidTraits();
1608  pidTof.setMatchFlag(FinalMatchedCellsVec[ii].matchFlag);
1609  pidTof.setYLocal(dy);
1610  pidTof.setZLocal(FinalMatchedCellsVec[ii].zhit);
1611  pidTof.setPosition(FinalMatchedCellsVec[ii].hitPosition);
1612  pidTof.setThetaLocal(FinalMatchedCellsVec[ii].theta);
1613  globalTrack->setBTofPidTraits(pidTof);
1614 
1615  if(primaryTrack) {
1616  StMuBTofPidTraits ppidTof = primaryTrack->btofPidTraits();
1617  ppidTof.setMatchFlag(FinalMatchedCellsVec[ii].matchFlag);
1618  ppidTof.setYLocal(dy);
1619  ppidTof.setZLocal(FinalMatchedCellsVec[ii].zhit);
1620  ppidTof.setPosition(FinalMatchedCellsVec[ii].hitPosition);
1621  ppidTof.setThetaLocal(FinalMatchedCellsVec[ii].theta);
1622  primaryTrack->setBTofPidTraits(ppidTof);
1623  }
1624 
1625  } // end final matched cells
1626 
1627  if(mHisto) {
1628  mCellsPrimaryPerEventMatch3->Fill(nValidSinglePrimHitCells);
1629  }
1630 
1631  if(Debug()) { LOG_INFO << "F: before/after" << FinalMatchedCellsVec.size() << "/" <<nValidSinglePrimHitCells << endm; }
1632  // end of Sect.F
1633  if (doPrintCpuInfo) {
1634  timer.stop();
1635  LOG_INFO << "CPU time after Step F - final : "
1636  << timer.elapsedTime() << " sec" << endm;
1637  timer.start();
1638  }
1639 
1640  LOG_INFO << " #(daq hits): " << daqCellsHitVec.size()
1641  << "\t#(proj hits): " << allCellsHitVec.size()
1642  << "\t#(prim proj hits): " << nPrimaryHits
1643  << "\n#(matched hits): " << FinalMatchedCellsVec.size()
1644  << "\n#(single valid hits): " << nValidSingleHitCells
1645  << "\t#(single prim valid hits): " << nValidSinglePrimHitCells
1646  << endm;
1647 
1648  if (doPrintMemoryInfo) {
1649  StMemoryInfo::instance()->snapshot();
1650  StMemoryInfo::instance()->print();
1651  }
1652  if (doPrintCpuInfo) {
1653  timer.stop();
1654  LOG_INFO << "CPU time for StBTofMatchMaker::Make(): "
1655  << timer.elapsedTime() << " sec" << endm;
1656  }
1657 
1658  LOG_INFO << "StBTofMatchMaker -- bye-bye" << endm;
1659 
1660  return;
1661 }
1662 
1663 //---------------------------------------------------------------------------
1664 // Book histograms and create ordered collection for easy manipulation
1665 void StBTofMatchMaker::bookHistograms(void){
1666 
1667  mEventCounterHisto = new TH1D("eventCounter","eventCounter",20,0,20);
1668 
1669  mCellsMultInEvent = new TH1D("cellsPerEvent","cellsPerEvent",1000,0,1000);
1670  mHitsMultInEvent = new TH1D("hitsPerEvent","hitsPerEvent",1000,0,1000);
1671  mHitsPrimaryInEvent = new TH1D("hitsPrimaryPerEvent","hitsPrimaryPerEvent",1000,0,1000);
1672  mHitsMultPerTrack = new TH1D("hitsPerTrack","hitsPerTrack",10,0,10);
1673  mHitsPosition = new TH2D("hitsPosition","hitsPositions",300,-15.,15.,200,-5.,5.);
1674 
1675  // occupancy
1676  for(int i=0;i<mNTray;i++) {
1677  char hisname[100];
1678  sprintf(hisname,"Occupancy_Tray_%d",i+1);
1679  mDaqOccupancy[i] = new TH1D(hisname,"",192,0,192);
1680  sprintf(hisname,"OccupancyProj_Tray_%d",i+1);
1681  mDaqOccupancyProj[i] = new TH1D(hisname,"",192,0,192);
1682  }
1683 
1684  // correlation
1685  for(int i=0;i<mNTray;i++) {
1686  char hisname[100];
1687  sprintf(hisname,"Corr_Tray_%d",i+1);
1688  mHitCorr[i] = new TH2D(hisname,"",192,0,192,192,0,192);
1689  sprintf(hisname,"Corr_Tray_%d_module",i+1);
1690  mHitCorrModule[i] = new TH2D(hisname,"",32,0,32,32,0,32);
1691  }
1692 
1693  // project hit position
1694  for(int i=0;i<mNTray;i++) {
1695  char hisname[100];
1696  sprintf(hisname,"LocalYZ_Tray_%d",i+1);
1697  mDeltaHitFinal[i] = new TH2D(hisname,"",300,-15.,15.,200,-5.,5.);
1698  }
1699 
1700  // TPC tracks
1701  if(mSaveTree) {
1702  mTrackTree = new TTree("track","track");
1703  mTrackTree->Branch("pt",&trackTree.pt,"pt/F");
1704  mTrackTree->Branch("eta",&trackTree.eta,"eta/F");
1705  mTrackTree->Branch("phi",&trackTree.phi,"phi/F");
1706  mTrackTree->Branch("nfitpts",&trackTree.nfitpts,"nfitpts/I");
1707  mTrackTree->Branch("dEdx",&trackTree.dEdx,"dEdx/F");
1708  mTrackTree->Branch("ndEdxpts",&trackTree.ndEdxpts,"ndEdxpts/I");
1709  mTrackTree->Branch("charge",&trackTree.charge,"charge/I");
1710  mTrackTree->Branch("projTrayId",&trackTree.projTrayId,"projTrayId/I");
1711  mTrackTree->Branch("projCellChan",&trackTree.projCellChan,"projCellChan/I");
1712  mTrackTree->Branch("projY",&trackTree.projY,"projY/F");
1713  mTrackTree->Branch("projZ",&trackTree.projZ,"projZ/F");
1714  }
1715 
1716  mTrackPtEta = new TH2D("trackPtEta","",100,0.,5.,60,-1.5,1.5);
1717  //mTrackPtPhi = new TH2D("trackPtPhi","",100,0.,5.,120,0.,2*3.14159);
1718  mTrackPtPhi = new TH2D("trackPtPhi","",100,0.,5.,120,0.,2.*M_PI);
1719  mTrackNFitPts = new TH1D("trackNFitPts","",50,0.,50.);
1720  mTrackdEdxvsp = new TH2D("trackdEdxvsp","",500,0.,5.,1000,0.,10.);
1721  mNSigmaPivsPt = new TH2D("nSigmaPivsPt","",500,0.,5.,1000,-10.,10.);
1722 
1723  // association
1724  mCellsPerEventMatch1 = new TH1D("cellsPerEventMatch1","cellPerEventMatch1",100,0,100);
1725  mHitsPerEventMatch1 = new TH1D("hitsPerEventMatch1","hitsPerEventMatch1",100,0,100);
1726  mCellsPerTrackMatch1 = new TH1D("cellsPerTrackMatch1","cellsPerTrackMatch1",100,0,100);
1727  mTracksPerCellMatch1 = new TH1D("tracksPerCellMatch1","tracksPerCellMatch1",100,0,100);
1728  mDeltaHitMatch1 = new TH2D("deltaHitMatch1","deltaHitMatch1",300,-15,15,200,-5.,5.);
1729 
1730  // kick out multi-hit
1731  mCellsPerEventMatch2 = new TH1D("cellsPerEventMatch2","cellPerEventMatch2",100,0,100);
1732  mHitsPerEventMatch2 = new TH1D("hitsPerEventMatch2","hitsPerEventMatch2",100,0,100);
1733  mCellsPerTrackMatch2 = new TH1D("cellsPerTrackMatch2","cellsPerTrackMatch2",100,0,100);
1734  mTracksPerCellMatch2 = new TH1D("tracksPerCellMatch2","tracksPerCellMatch2",100,0,100);
1735  mDeltaHitMatch2 = new TH2D("deltaHitMatch2","deltaHitMatch2",300,-15,15,200,-5.,5.);
1736 
1737  // sort out multi matched cells
1738  mCellsPerEventMatch3 = new TH1D("cellsPerEventMatch3","cellsPerEventMatch3",100,0,100);
1739  mHitsPerEventMatch3 = new TH1D("hitsPerEventMatch3","hitsPerEventMatch3",100,0,100);
1740  mCellsPerTrackMatch3 = new TH1D("cellsPerTrackMatch3","cellsPerTrackMatch3",100,0,100);
1741  mTracksPerCellMatch3 = new TH1D("tracksPerCellMatch3","tracksPerCellMatch3",100,0,100);
1742  mDeltaHitMatch3 = new TH2D("deltaHitMatch3","deltaHitMatch3",300,-15,15,200,-5.,5.);
1743 
1744  mCellsPrimaryPerEventMatch3 = new TH1D("cellsPrimaryPerEventMatch3","cellsPrimaryPerEventMatch3",100,0,100);
1745 
1746  return;
1747 }
1748 
1749 
1750 //---------------------------------------------------------------------------
1751 // store histograms in a seperate root file
1752 void StBTofMatchMaker::writeHistogramsToFile(){
1753  // Output file
1754  TFile *theHistoFile = new TFile(mHistoFileName.c_str(), "RECREATE");
1755  LOG_INFO << "StBTofMatchMaker::writeHistogramsToFile()"
1756  << " histogram file " << mHistoFileName << endm;
1757 
1758  theHistoFile->cd();
1759 
1760  if(mHisto) {
1761 
1762 
1763  for(int i=0;i<mNTray;i++) {
1764  mDaqOccupancy[i]->Write();
1765  mDaqOccupancyProj[i]->Write();
1766  mHitCorr[i]->Write();
1767  mHitCorrModule[i]->Write();
1768  mDeltaHitFinal[i]->Write();
1769  }
1770 
1771  mEventCounterHisto->Write();
1772  mCellsMultInEvent->Write();
1773  mHitsMultInEvent->Write();
1774  mHitsPrimaryInEvent->Write();
1775  mHitsMultPerTrack->Write();
1776  mHitsPosition->Write();
1777 
1778  mTrackPtEta->Write();
1779  mTrackPtPhi->Write();
1780  mTrackNFitPts->Write();
1781  mTrackdEdxvsp->Write();
1782  mNSigmaPivsPt->Write();
1783 
1784  mCellsPerEventMatch1->Write();
1785  mHitsPerEventMatch1->Write();
1786  mCellsPerTrackMatch1->Write();
1787  mTracksPerCellMatch1->Write();
1788  mDeltaHitMatch1->Write();
1789 
1790  mCellsPerEventMatch2->Write();
1791  mHitsPerEventMatch2->Write();
1792  mCellsPerTrackMatch2->Write();
1793  mTracksPerCellMatch2->Write();
1794  mDeltaHitMatch2->Write();
1795 
1796  mCellsPerEventMatch3->Write();
1797  mHitsPerEventMatch3->Write();
1798  mCellsPerTrackMatch3->Write();
1799  mTracksPerCellMatch3->Write();
1800  mDeltaHitMatch3->Write();
1801 
1802  mCellsPrimaryPerEventMatch3->Write();
1803 
1804  theHistoFile->Write();
1805  theHistoFile->Close();
1806 
1807  }
1808 
1809  if(mSaveTree) {
1810  TFile *theTreeFile = new TFile((mHistoFileName+".tree.root").c_str(), "RECREATE");
1811  theTreeFile->cd();
1812  mTrackTree->Write();
1813  theTreeFile->Write();
1814  theTreeFile->Close();
1815  }
1816 
1817  return;
1818 }
1819 
1820 //---------------------------------------------------------------------------
1821 // determine whether this is a valid TPC track
1822 bool StBTofMatchMaker::validTrack(StTrack *track){
1823  // 1. no track, no go.
1824  if (!track) return false;
1825 
1826  // 2. track quality flag, should be >0
1827  if (track->flag()<=0 || track->flag()>=1000) return false;
1828 
1829  // 3. minimum #hits per track - obsolete
1830  // if (track->topologyMap().numberOfHits(kTpcId) < mMinHitsPerTrack) return false;
1831  // 4. minimum #fit points per track
1832  if (track->fitTraits().numberOfFitPoints(kTpcId) < mMinFitPointsPerTrack) return false;
1833  // 5. minimum #fit points over #maximum points
1834  //fg float ratio = (1.0*track->fitTraits().numberOfFitPoints(kTpcId)) / (1.0*track->numberOfPossiblePoints(kTpcId));
1835  float ratio = (float)track->fitTraits().numberOfFitPoints(kTpcId) / (1.0*track->numberOfPossiblePoints(kTpcId));
1836  if (ratio < mMinFitPointsOverMax) return false;
1837 
1838  return true;
1839 }
1840 
1841 //---------------------------------------------------------------------------
1842 // determine whether this is a valid TPC track
1843 bool StBTofMatchMaker::validTrack(StMuTrack *track){
1844  // 1. no track, no go.
1845  if (!track) return false;
1846 
1847  // 2. track quality flag, should be >0
1848  if (track->flag()<=0 || track->flag()>=1000) return false;
1849 
1850  // 3. minimum #hits per track - obsolete
1851  // if (track->nHits() < mMinHitsPerTrack) return false;
1852  // 4. minimum #fit points per track
1853  if (track->nHitsFit(kTpcId) < mMinFitPointsPerTrack) return false;
1854  // 5. minimum #fit points over #maximum points
1855  //fg float ratio = (1.0*track->fitTraits().numberOfFitPoints(kTpcId)) / (1.0*track->numberOfPossiblePoints(kTpcId));
1856  float ratio = (float)track->nHitsFit(kTpcId) / (1.0*track->nHitsPoss(kTpcId));
1857  if (ratio < mMinFitPointsOverMax) return false;
1858 
1859  return true;
1860 }
1861 
1862 //---------------------------------------------------------------------------
1863 // returns the proper track geometry, based on a global user setting
1864 StTrackGeometry* StBTofMatchMaker::trackGeometry(StTrack* track){
1865  // returns apropriate StTrackGeometry (standard or outerGeometry)
1866  if (!track) return 0;
1867  StTrackGeometry *thisTrackGeometry;
1868  if (mOuterTrackGeometry)
1869  thisTrackGeometry = track->outerGeometry();
1870  else
1871  thisTrackGeometry = track->geometry();
1872  return thisTrackGeometry;
1873 }
StMaker::m_Mode
Int_t m_Mode
counters
Definition: StMaker.h:81
StBTofMatchMaker::setMaxDCA
void setMaxDCA(Float_t)
set maximum distance of closest approach
Definition: StBTofMatchMaker.h:342
StBTofGeometry::HelixCrossCellIds
Bool_t HelixCrossCellIds(const StHelixD &helix, IntVec &idVec, DoubleVec &pathVec, PointVec &crossVec) const
Definition: StBTofGeometry.cxx:1531
StMuDst::btofHit
static StMuBTofHit * btofHit(int i)
returns pointer to the i-th muBTofHit
Definition: StMuDst.h:419
StMuDst::globalTracks
static TObjArray * globalTracks()
returns pointer to the global tracks list
Definition: StMuDst.h:303
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
StMuTrack::index2Global
Int_t index2Global() const
Returns index of associated global track. If not in order can be set with StMuDst::fixTrackIndeces() ...
Definition: StMuTrack.h:231
StBTofMatchMaker::setAlignFileName
void setAlignFileName(const Char_t *infile="")
input file for alignment parameters
Definition: StBTofMatchMaker.h:356
StBTofMatchMaker::Make
Int_t Make()
Main match algorithm.
Definition: StBTofMatchMaker.cxx:297
StMuTrack::id
short id() const
Returns the track id(or key), is unique for a track node, i.e. global and primary tracks have the sam...
Definition: StMuTrack.h:228
StHelix::x
double x(double s) const
coordinates of helix at point s
Definition: StHelix.hh:182
StMuTrack::nHitsDedx
UShort_t nHitsDedx() const
Return number of hits used for dEdx.
Definition: StMuTrack.h:238
StBTofMatchMaker::setMinHitsPerTrack
void setMinHitsPerTrack(Int_t)
set minimum hits per track
Definition: StBTofMatchMaker.h:336
StBTofMatchMaker::setMinFitPointsPerTrack
void setMinFitPointsPerTrack(Int_t)
set minimum fit points per track
Definition: StBTofMatchMaker.h:338
StBTofMatchMaker::Finish
Int_t Finish()
Print run summary, and write QA histograms.
Definition: StBTofMatchMaker.cxx:280
StBTofMatchMaker::Init
Int_t Init()
process start-up options
Definition: StBTofMatchMaker.cxx:190
StBTofMatchMaker::setCreateTreeFlag
void setCreateTreeFlag(Bool_t tree=kTRUE)
enable track-tree filling
Definition: StBTofMatchMaker.h:348
StMuTrack::helix
StPhysicalHelixD helix() const
Returns inner helix (first measured point)
Definition: StMuTrack.cxx:407
StMuTrack::nHitsFit
UShort_t nHitsFit() const
Return total number of hits used in fit.
Definition: StMuTrack.h:239
StBTofMatchMaker::InitRun
Int_t InitRun(Int_t)
initialize DaqMap, Geometry, and INL
Definition: StBTofMatchMaker.cxx:231
StMuTrack::flag
short flag() const
Returns flag, (see StEvent manual for type information)
Definition: StMuTrack.h:230
StMuTrack::charge
Short_t charge() const
Returns charge.
Definition: StMuTrack.h:255
StBTofMatchMaker::setOuterTrackGeometry
void setOuterTrackGeometry()
selection of inner or outer geometry. By default - outerGeometry
Definition: StBTofMatchMaker.h:332
StBTofMatchMaker::setHistoFileName
void setHistoFileName(const Char_t *)
set histogram output file name
Definition: StBTofMatchMaker.h:344
StMuTrack::nSigmaPion
Double_t nSigmaPion() const
Returns Craig&#39;s distance to the calculated dE/dx band for pions in units of sigma.
Definition: StMuTrack.h:245
StBTofMatchMaker::setSaveGeometry
void setSaveGeometry(const Bool_t geomSave=kFALSE)
save geometry if it will be used by following makers in the chain
Definition: StBTofMatchMaker.h:350
StMuTrack::nHitsPoss
UShort_t nHitsPoss() const
Return number of possible hits on track.
Definition: StMuTrack.cxx:242
StMuTrack::momentum
const StThreeVectorF & momentum() const
Returns 3-momentum at dca to primary vertex.
Definition: StMuTrack.h:260
kStOK
Definition: Stypes.h:40
StMuTrack::dEdx
Double_t dEdx() const
Returns measured dE/dx value.
Definition: StMuTrack.h:248
StBTofMatchMaker::setMinFitPointsOverMax
void setMinFitPointsOverMax(Float_t)
set minimum fit-points/max-points ratio
Definition: StBTofMatchMaker.h:340
StBTofMatchMaker::StBTofMatchMaker
StBTofMatchMaker(const Char_t *name="btofMatch")
Default constructor.
Definition: StBTofMatchMaker.cxx:151
StMuDst::array
static TClonesArray * array(int type)
returns pointer to the n-th TClonesArray
Definition: StMuDst.h:262
StBTofMatchMaker::setCreateHistoFlag
void setCreateHistoFlag(Bool_t histos=kTRUE)
enable QA histogram filling
Definition: StBTofMatchMaker.h:346
StMuTrack::outerHelix
StPhysicalHelixD outerHelix() const
Returns outer helix (last measured point)
Definition: StMuTrack.cxx:412
Generated by   doxygen 1.8.5