StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StETofMatchMaker.cxx
1 /***************************************************************************
2  *
3  * $Id: StETofMatchMaker.cxx,v 1.10 2020/01/18 02:37:10 fseck Exp $
4  *
5  * Author: Florian Seck, April 2018
6  ***************************************************************************
7  *
8  * Description: StETofMatchMaker - class to match StETofHits to tracks.
9  * The matching is done in several steps:
10  * - get a list of eTOF hits
11  * - check for each track if its helix has an intersection with
12  * the eTOF active volumes (MRPCs gas gaps)
13  * - resolve matching ambiguities
14  *
15  ***************************************************************************
16  *
17  * $Log: StETofMatchMaker.cxx,v $
18  * Revision 1.10 2020/01/18 02:37:10 fseck
19  * fixing StEvent part of eTOF-only T0 calculation
20  *
21  * Revision 1.9 2020/01/16 03:53:41 fseck
22  * added etof-only and hybrid btof-etof start time calculations for on-the-fly corrections
23  *
24  * Revision 1.8 2019/12/17 03:28:01 fseck
25  * update to histograms for .hist.root files
26  *
27  * Revision 1.7 2019/12/10 16:00:34 fseck
28  * possibility to use step-wise track extrapolation in changing magnetic field via setting a flag
29  *
30  * Revision 1.6 2019/05/09 00:02:46 fseck
31  * match distances as member variables and updated handling for many-tracks-to-one-hit matches
32  *
33  * Revision 1.5 2019/04/24 02:33:48 fseck
34  * start time fix to previous commit
35  *
36  * Revision 1.4 2019/04/24 01:02:11 fseck
37  * fix to start time for simulation and more histograms added to doQA mode
38  *
39  * Revision 1.3 2019/03/25 01:05:48 fseck
40  * added more histograms for offline QA
41  *
42  * Revision 1.2 2019/03/08 19:09:31 fseck
43  * added a few eTOF histograms to the .hist.root files for offline QA
44  *
45  * Revision 1.1 2019/02/19 19:52:28 jeromel
46  * Reviewed code provided by F.Seck
47  *
48  *
49  ***************************************************************************/
50 #include <sstream>
51 #include <cmath>
52 
53 #include "TFile.h"
54 #include "TH1F.h"
55 #include "TH2F.h"
56 
57 #include "StParticleTypes.hh"
58 #include "StParticleDefinition.hh"
59 #include "SystemOfUnits.h"
60 #include "PhysicalConstants.h"
61 
62 #include "StChainOpt.h" // for renaming the histogram file
63 
64 #include "StEvent.h"
65 #include "StTrackNode.h"
66 #include "StTrackGeometry.h"
67 #include "StGlobalTrack.h"
68 #include "StPrimaryTrack.h"
69 #include "StPrimaryVertex.h"
70 #include "StETofCollection.h"
71 #include "StETofHit.h"
72 #include "StTrackPidTraits.h"
73 #include "StETofPidTraits.h"
74 #include "StTpcDedxPidAlgorithm.h"
75 #include "StDedxPidTraits.h"
76 
77 #include "StDetectorDbMaker/St_MagFactorC.h"
78 
79 #include "StBTofCollection.h"
80 #include "StBTofHeader.h"
81 
82 #include "StMuDSTMaker/COMMON/StMuDst.h"
83 #include "StMuDSTMaker/COMMON/StMuArrays.h"
84 #include "StMuDSTMaker/COMMON/StMuTrack.h"
85 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
86 #include "StMuDSTMaker/COMMON/StMuETofHit.h"
87 #include "StMuDSTMaker/COMMON/StMuETofPidTraits.h"
88 #include "StMuDSTMaker/COMMON/StMuETofDigi.h"
89 #include "StMuDSTMaker/COMMON/StMuETofHeader.h"
90 
91 #include "StETofMatchMaker.h"
92 #include "StETofHitMaker/StETofHitMaker.h"
93 #include "StETofUtil/StETofGeometry.h"
94 #include "StETofUtil/StETofConstants.h"
95 
96 #include "tables/St_etofMatchParam_Table.h"
97 
98 
99 // *******************************************************************************************************
100 
101 namespace etofProjection {
102  // safety margins in cm in local x and y direction for track extrapolation to etof modules
103  const double safetyMargins[ 2 ] = { 2., 2. };
104 
105  const float minEtaCut = 0.0;
106 
107  const float minEtaProjCut = -1.0;
108  const float maxEtaProjCut = -1.7;
109 
110 
111  // Legacy hardcoded alignment. Moved to database.
112  const double deltaXoffset[ 108 ] = { 0 };
113 
114  const double deltaYoffset[ 108 ] = { 0 };
115 
116  const double deltaRcut = 4.;
117 }
118 
119 //---------------------------------
120 static const double kaon_minus_mass_c2 = 493.68 * MeV;
121 //---------------------------------
122 
123 namespace etofHybridT0 {
124  const unsigned int minCand = 2;
125  const unsigned int nSwitch = 5;
126  const unsigned int nMin = 10;
127  const unsigned int nUpdate = 5;
128  const float lowCut = 0.10; // 10%
129  const float highCut = 0.85; // 85%
130  const float diffToClear = 2.; // ns
131  const float biasOffset = 0.075; // ns
132 }
133 // *******************************************************************************************************
134 
135 
136 //---------------------------------------------------------------------------
137 StETofMatchMaker::StETofMatchMaker( const char* name )
138 : StMaker( "etofMatch", name ),
139  mEvent( nullptr ),
140  mMuDst( nullptr ),
141  mETofGeom( nullptr ),
142  mFileNameMatchParam( "" ),
143  mFileNameAlignParam( "" ),
144  mIsStEventIn( false ),
145  mIsMuDstIn( false ),
146  mOuterTrackGeometry( true ),
147  mUseHelixSwimmer( false ),
148  mUseOnlyBTofHeaderStartTime( true ),
149  mIsSim( false ),
150  mDoQA( false ),
151  mDebug( false ),
152  mMatchDistX( 5. ),
153  mMatchDistY( 10. ),
154  mMatchDistT( 99999. ),
155  mT0corrVec(),
156  mT0corr( 0. ),
157  mT0switch( 0 ),
158  mNupdatesT0( 0 ),
159  mIndex2Primary(),
160  mMatchRadius( 0. ),
161  mLocalYmax(16.),
162  mClockJumpCand(),
163  mClockJumpDirection(),
164  mHistFileName( "" ),
165  mHistograms(),
166  mHistograms2d(),
167  dx_3sig(2.5),
168  dy_3sig(4.0),
169  dt_3sig(0.22),
170  dy_max(5.0)
171 {
172  mT0corrVec.reserve( 500 );
173  mTrackCuts.push_back( 0. ); // nHitsFit
174  mTrackCuts.push_back( 0. ); // nHitsRatio
175  mTrackCuts.push_back( 0. ); // low pt
176 }
177 
178 
179 
180 //---------------------------------------------------------------------------
181 StETofMatchMaker::~StETofMatchMaker()
182 {
183  /* nope */
184 }
185 
186 
187 //---------------------------------------------------------------------------
188 Int_t
189 StETofMatchMaker::Init()
190 {
191  LOG_INFO << "StETofMatchMaker::Init()" << endm;
192 
193  LOG_INFO << "isSimulation flag was set to: " << mIsSim << endm;
194 
195  bookHistograms();
196 
197  return kStOk;
198 }
199 
200 
201 //---------------------------------------------------------------------------
202 Int_t
203 StETofMatchMaker::InitRun( Int_t runnumber )
204 {
205  LOG_INFO << "StETofMatchMaker::InitRun()" << endm;
206 
207  TDataSet* dbDataSet = nullptr;
208  std::ifstream paramFile;
209 
210  // --------------------------------------------------------------------------------------------
211  // initialize hit building parameters from parameter file (if filename is provided) or database:
212  // -- match param
213  // --------------------------------------------------------------------------------------------
214 
215  // match param
216  if( mFileNameMatchParam.empty() ) {
217  LOG_INFO << "etofMatchParam: no filename provided --> load database table" << endm;
218 
219  dbDataSet = GetDataBase( "Calibrations/etof/etofMatchParam" );
220 
221  St_etofMatchParam* etofMatchParam = static_cast< St_etofMatchParam* > ( dbDataSet->Find( "etofMatchParam" ) );
222  if( !etofMatchParam ) {
223  LOG_ERROR << "unable to get the match params from the database" << endm;
224  return kStFatal;
225  }
226 
227  etofMatchParam_st* matchParamTable = etofMatchParam->GetTable();
228 
229  mMatchRadius = matchParamTable->matchRadius;
230 
231  mTrackCuts.at( 0 ) = matchParamTable->trackCutNHitsFit;
232  mTrackCuts.at( 1 ) = matchParamTable->trackCutNHitsRatio;
233  mTrackCuts.at( 2 ) = matchParamTable->trackCutLowPt;
234 
235  }
236  else {
237  LOG_INFO << "etofMatchParam: filename provided --> use parameter file: " << mFileNameMatchParam.c_str() << endm;
238 
239  paramFile.open( mFileNameMatchParam.c_str() );
240 
241  if( !paramFile.is_open() ) {
242  LOG_ERROR << "unable to get the 'etofMatchParam' parameters from file --> file does not exist" << endm;
243  return kStFatal;
244  }
245 
246  std::vector< float > param;
247  float temp;
248  while( paramFile >> temp ) {
249  param.push_back( temp );
250  }
251 
252  paramFile.close();
253 
254  if( param.size() != 4 ) {
255  LOG_ERROR << "parameter file for 'etofMatchParam' has not the right amount of entries: ";
256  LOG_ERROR << param.size() << " instead of 4 !!!!" << endm;
257  return kStFatal;
258  }
259 
260  mMatchRadius = param.at( 0 );
261 
262  mTrackCuts.at( 0 ) = param.at( 1 );
263  mTrackCuts.at( 1 ) = param.at( 2 );
264  mTrackCuts.at( 2 ) = param.at( 3 );
265 
266  }
267 
268  LOG_INFO << " match radius: " << mMatchRadius << endm;
269  LOG_INFO << " track cut (nHitsFit): " << mTrackCuts.at( 0 ) << endm;
270  LOG_INFO << " track cut (nHitsRatio): " << mTrackCuts.at( 1 ) << endm;
271  LOG_INFO << " track cut (low pt): " << mTrackCuts.at( 2 ) << endm;
272 
273  // --------------------------------------------------------------------------------------------
274 
275 
276  // --------------------------------------------------------------------------------------------
277  // initializie etof geometry
278  // --------------------------------------------------------------------------------------------
279  if( !mETofGeom ) {
280  LOG_INFO << " creating a new eTOF geometry . . . " << endm;
281  mETofGeom = new StETofGeometry( "etofGeometry", "etofGeometry in MatchMaker" );
282  }
283 
284  if( mETofGeom && !mETofGeom->isInitDone() ) {
285  LOG_INFO << " eTOF geometry initialization ... " << endm;
286 
287  if( !gGeoManager ) GetDataBase( "VmcGeometry" );
288 
289  if( !gGeoManager ) {
290  LOG_ERROR << "Cannot get GeoManager" << endm;
291  return kStFatal;
292  }
293 
294  LOG_DEBUG << " gGeoManager: " << gGeoManager << endm;
295 
296  if (mFileNameAlignParam != ""){
297  LOG_DEBUG << " gGeoManager: Setting alignment file: " << mFileNameAlignParam << endm;
298  mETofGeom->setFileNameAlignParam(mFileNameAlignParam);
299  }
300 
301  mETofGeom->init( gGeoManager, etofProjection::safetyMargins, mUseHelixSwimmer ); //provide backline to initiating maker to load DB tables
302  }
303 
304  if ( mETofGeom && !mETofGeom->isInitDone() ) { //if initialization attempt above failed.
305  LOG_ERROR << "Initialization of StEtofGeometry failed" << endm;
306  return kStFatal;
307  }
308 
309 
310  // --------------------------------------------------------------------------------------------
311  // pointer to eTOF hit maker
312  // --------------------------------------------------------------------------------------------
313 
314  mETofHitMaker = ( StETofHitMaker* ) GetMaker( "etofHit" );
315  LOG_DEBUG << "StETofMatchMaker::InitRun() -- pointer to eTOF hit maker: " << mETofHitMaker << endm;
316 
317  if( mDoQA ) {
318  // for geometry debugging
319  for( unsigned int i=0; i<mETofGeom->nValidModules(); i++ ) {
320 
321  StThreeVectorF pos = mETofGeom->module( i )->centerPos();
322 
323  int sector = mETofGeom->module( i )->sector();
324  int plane = mETofGeom->module( i )->plane();
325 
326  for( int j=0; j<sector; j++ ) {
327  mHistograms.at( "eTofSectors" )->Fill( pos.x(), pos.y() );
328  }
329  for( int j=0; j<plane; j++ ) {
330  mHistograms.at( "eTofModules" )->Fill( pos.x(), pos.y() );
331  }
332  }
333 
334 
335  // for magnetic field debugging
336  int nBinsX = mHistograms2d.at( "bfield_z" )->GetNbinsX();
337  int nBinsY = mHistograms2d.at( "bfield_z" )->GetNbinsY();
338  for( int i=0; i<nBinsX; i++ ) {
339  for( int j=0; j<nBinsY; j++ ) {
340  double z = mHistograms2d.at( "bfield_z" )->GetXaxis()->GetBinCenter( i );
341  double y = mHistograms2d.at( "bfield_z" )->GetYaxis()->GetBinCenter( j );
342 
343  mHistograms2d.at( "bfield_z" )->Fill( z, y, mETofGeom->getFieldZ( 0., y, z ) );
344  }
345  }
346  }
347 
348  return kStOk;
349 }
350 
351 
352 //---------------------------------------------------------------------------
353 Int_t
354 StETofMatchMaker::FinishRun( Int_t runnumber )
355 {
356  LOG_DEBUG << "StETofMatchMaker::FinishRun() -- cleaning up the geometry" << endm;
357 
358  if( mETofGeom ) {
359  mETofGeom->reset();
360  }
361 
362  return kStOk;
363 }
364 
365 
366 //---------------------------------------------------------------------------
367 Int_t
368 StETofMatchMaker::Finish()
369 {
370  LOG_DEBUG << "StETofMatchMaker::Finish()" << endm;
371 
372  if( mDoQA ) {
373  LOG_INFO << "Finish() - writing *.etofMatch.root ..." << endm;
374  setHistFileName();
375  writeHistograms();
376  }
377 
378  if( mETofGeom ) {
379  delete mETofGeom;
380  mETofGeom = nullptr;
381  }
382 
383  return kStOk;
384 }
385 
386 
387 //---------------------------------------------------------------------------
388 Int_t
389 StETofMatchMaker::Make()
390 {
391  LOG_DEBUG << "StETofMatchMaker::Make(): starting ..." << endm;
392 
393  mIsStEventIn = false;
394  mIsMuDstIn = false;
395 
396  mEvent = ( StEvent* ) GetInputDS( "StEvent" );
397  // mEvent = NULL; //don't check for StEvent for genDst.C testing. PW
398 
399  if ( mEvent ) {
400  LOG_DEBUG << "Make() - running on StEvent" << endm;
401 
402  StETofCollection* etofCollection = mEvent->etofCollection();
403 
404  if( !etofCollection ) { //additional check for empty StEvents structures produced by other Makers. Needed for genDst.C
405  LOG_WARN << "Make() - Found StEvent data structure, but no eTOF collection. Try MuDst processing instead" << endm;
406  mMuDst = ( StMuDst* ) GetInputDS( "MuDst" );
407 
408  if( mMuDst ) {
409  LOG_DEBUG << "Make() - running on MuDsts" << endm;
410 
411  mIsMuDstIn = true;
412 
413  fillIndexToPrimaryMap();
414 
415  cleanUpTraits();
416  }
417  }else{
418  mIsStEventIn = true;
419 
420  cleanUpTraits();
421  }
422  }
423  else {
424  LOG_DEBUG << "Make(): no StEvent found" << endm;
425 
426  mMuDst = ( StMuDst* ) GetInputDS( "MuDst" );
427 
428  if( mMuDst ) {
429  LOG_DEBUG << "Make() - running on MuDsts" << endm;
430 
431  mIsMuDstIn = true;
432 
433  fillIndexToPrimaryMap();
434 
435  cleanUpTraits();
436  }
437  else {
438  LOG_DEBUG << "Make() - no StMuDst found" << endm;
439  return kStOk;
440  }
441  }
442 
443  if( !mIsStEventIn && !mIsMuDstIn ) {
444  LOG_WARN << "Make() - neither StEvent nor MuDst available ... bye-bye" << endm;
445  return kStOk;
446  }
447 
448  if( mDoQA ) {
449  mHistograms.at( "eventCounter" )->Fill( 1 );
450  }
451 
452  //.........................................................................
453  // A. read data from StETofHit & build vector of candidate hits
454  //
455  eTofHitVec detectorHitVec;
456 
457  readETofDetectorHits( detectorHitVec );
458 
459  if( detectorHitVec.size() == 0 ) {
460  //LOG_INFO << "Make() -- event done ... bye-bye" << endm;
461 
462  return kStOk;
463  }
464 
465  //.........................................................................
466  // B. loop over global tracks & determine all track intersections with active eTof volumes
467  //
468  eTofHitVec intersectionVec;
469  int nPrimaryWithIntersection = 0;
470 
471  float bFieldFromGeom = mETofGeom->getFieldZ( 100., 100., 0. );
472  float bField = 0;
473  if( mIsMuDstIn ) {
474  bField = mMuDst->event()->runInfo().magneticField();
475  }
476  else {
477  bField = mEvent->runInfo()->magneticField();
478  }
479 
480  if( mUseHelixSwimmer && fabs( bFieldFromGeom - bField ) > 0.2 ) {
481  LOG_WARN << "Wrong magnetc field bField = " << bField << " bFieldFromGeom = " << bFieldFromGeom << " --> check the magF input!" << endm;
482  }
483 
484  findTrackIntersections( intersectionVec, nPrimaryWithIntersection );
485 
486  if( intersectionVec.size() == 0 ) {
487  //LOG_INFO << "Make() -- event done ... bye-bye" << endm;
488 
489  return kStOk;
490  }
491 
492  mHistograms.at( "intersectionMult_etofMult" )->Fill( detectorHitVec.size(), intersectionVec.size() );
493 
494  //.........................................................................
495  // C. match detector hits to track intersections
496  //
497  eTofHitVec matchCandVec;
498 
499  matchETofHits( detectorHitVec, intersectionVec, matchCandVec );
500 
501  if( matchCandVec.size() == 0 ) {
502  //LOG_INFO << "Make() -- event done ... bye-bye" << endm;
503 
504  return kStOk;
505  }
506 
507  //Single Sided Hit Matching and clustering
508  eTofHitVec finalMatchVec;
509  sortandcluster(matchCandVec , detectorHitVec , intersectionVec , finalMatchVec);
510 
511 
512  //.........................................................................
513  // D. sort matchCand vector and deal with (discard) hits matched by multiple tracks
514  //
515  // define the vectors that store matchCandidates with a single / multiple tracks pointing to an eTof hit
516  eTofHitVec singleTrackMatchVec;
517  vector< eTofHitVec > multiTrackMatchVec;
518 
519  //sortSingleMultipleHits( matchCandVec, singleTrackMatchVec, multiTrackMatchVec ); // old matching procedure
520 
521  //if( singleTrackMatchVec.size() == 0 ) {
522  //LOG_INFO << "Make() -- event done ... bye-bye" << endm
523  // return kStOk;
524  // }
525 
526  //.........................................................................
527  // E. sort singleTrackMatchVector for multiple hits associated to single tracks and determine the best match
528  //
529 
530  //finalizeMatching( singleTrackMatchVec, finalMatchVec ); // old matching procedure
531 
532  if( finalMatchVec.size() == 0 ) {
533  //LOG_INFO << "Make() -- event done ... bye-bye" << endm;
534 
535  return kStOk;
536  }
537  else{
538  //LOG_INFO << "Make() -- number of found matches of eTOF hits with tracks: " << finalMatchVec.size() << endm;
539  }
540 
541  //.........................................................................
542  // F. fill ETofPidTraits for global and primary tracks and assign associated track to hits
543  //
544  fillPidTraits( finalMatchVec );
545 
546  //.........................................................................
547  // G. calculate pid variables for primary tracks and update PidTraits
548  //
549  int nPrimaryWithPid = 0;
550 
551  calculatePidVariables( finalMatchVec, nPrimaryWithPid );
552 
553  mHistograms.at( "primaryIntersect_validMatch" )->Fill( nPrimaryWithIntersection, nPrimaryWithPid );
554 
555  //.........................................................................
556  // H. fill QA histograms
557  //
558  fillQaHistograms( finalMatchVec );
559 
560  fillSlewHistograms( finalMatchVec );
561 
562  checkClockJumps();
563 
564  //LOG_INFO << "Make() -- event done ... bye-bye" << endm;
565 
566  return kStOk;
567 }
568 
569 
570 void StETofMatchMaker::fillIndexToPrimaryMap()
571 {
572  // clear and fill index2Primary map
573  mIndex2Primary.clear();
574 
575  for( int i = 0; i < mMuDst->array( muPrimary )->GetEntries(); i++ ) {
576  StMuTrack* pTrack = ( StMuTrack* ) mMuDst->array( muPrimary )->UncheckedAt( i );
577  if( !pTrack ) {
578  continue;
579  }
580  int index2Global = pTrack->index2Global();
581  if( index2Global < 0 ) {
582  continue;
583  }
584  mIndex2Primary[ index2Global ] = i;
585  }
586 }
587 
588 
589 void StETofMatchMaker::cleanUpTraits()
590 {
591  // StEvent processing
592  if( mIsStEventIn ) {
593  StSPtrVecTrackNode& nodes = mEvent->trackNodes();
594  size_t nNodes = nodes.size();
595  for( size_t iNode = 0; iNode < nNodes; iNode++ ) {
596  StGlobalTrack* gTrack = dynamic_cast< StGlobalTrack* > ( nodes[ iNode ]->track( global ) );
597  if( !gTrack ) continue;
598 
599  //clean up any association done before
600  StSPtrVecTrackPidTraits& traits = gTrack->pidTraits();
601  for( auto it = traits.begin(); it != traits.end(); it++ ) {
602  if( ( *it )->detector() == kETofId ) {
603 
604  if( mDebug ) {
605  StETofPidTraits* etofTraits = ( StETofPidTraits* ) *it;
606 
607  LOG_INFO << "cleanUpTraits() - etof traits on global track " << iNode << " already exist" << endm;
608  LOG_INFO << "matchFlag: " << etofTraits->matchFlag() << " localX: " << etofTraits->localX() << " localY: " << etofTraits->localY();
609  LOG_INFO << " tof: " << etofTraits->timeOfFlight() << " pathlength: " << etofTraits->pathLength() << " beta: " << etofTraits->beta() << endm;
610 
611  if( etofTraits->etofHit() ) {
612  LOG_INFO << "time: " << etofTraits->etofHit()->time() << endm;
613  }
614  }
615 
616  traits.erase( it );
617 
618  if( mDebug ) {
619  LOG_INFO << " ... erased" << endm;
620  }
621 
622  break;
623  }
624  }
625 
626  StPrimaryTrack* pTrack = dynamic_cast< StPrimaryTrack* > ( gTrack->node()->track( primary ) );
627  if( pTrack ) {
628  //clean up any association done before
629  StSPtrVecTrackPidTraits& ptraits = pTrack->pidTraits();
630  for( auto it = ptraits.begin(); it != ptraits.end(); it++ ) {
631  if( ( *it )->detector() == kETofId ) {
632 
633  if( mDebug ) {
634  StETofPidTraits* etofTraits = ( StETofPidTraits* ) *it;
635 
636  LOG_INFO << "cleanUpTraits() - etof traits on primary track corresponding to node " << iNode << " already exist" << endm;
637  LOG_INFO << "matchFlag: " << etofTraits->matchFlag() << " localX: " << etofTraits->localX() << " localY: " << etofTraits->localY();
638  LOG_INFO << " tof: " << etofTraits->timeOfFlight() << " pathlength: " << etofTraits->pathLength() << " beta: " << etofTraits->beta() << endm;
639 
640  if( etofTraits->etofHit() ) {
641  LOG_INFO << "time: " << etofTraits->etofHit()->time() << endm;
642  }
643  }
644 
645  ptraits.erase( it );
646 
647  if( mDebug ) {
648  LOG_INFO << " ... erased" << endm;
649  }
650 
651  break;
652  }
653  }
654  }
655 
656  }
657 
658  size_t nHits = mEvent->etofCollection()->etofHits().size();
659  for( size_t i=0; i<nHits; i++ ) {
660  StETofHit* aHit = mEvent->etofCollection()->etofHits().at( i );
661  if( !aHit ) continue;
662  aHit->setAssociatedTrack( nullptr );
663  }
664 
665  }
666  else { // MuDst processing
667  size_t nNodes = mMuDst->numberOfGlobalTracks();
668  for( size_t iNode=0; iNode<nNodes; iNode++ ) {
669  StMuTrack* track = mMuDst->globalTracks( iNode );
670  if( !track ) continue;
671  if( track->index2ETofHit() < 0 ) continue;
672 
673  if( mDebug ) {
674  StMuETofPidTraits etofTraits = track->etofPidTraits();
675 
676  LOG_INFO << "cleanUpTraits() - etof traits on global track " << iNode << " already exist" << endm;
677  LOG_INFO << "matchFlag: " << etofTraits.matchFlag() << " localX: " << etofTraits.localX() << " localY: " << etofTraits.localY();
678  LOG_INFO << " tof: " << etofTraits.timeOfFlight() << " pathlength: " << etofTraits.pathLength() << " beta: " << etofTraits.beta() << endm;
679 
680  if( mMuDst->etofHit( track->index2ETofHit() ) ) {
681  LOG_INFO << "time: " << mMuDst->etofHit( track->index2ETofHit() )->time() << endm;
682  }
683  }
684 
685  //clean up any association done before
686  StMuETofPidTraits pidETof;
687  track->setETofPidTraits( pidETof );
688  track->setIndex2ETofHit( -1 );
689 
690  if( mDebug ) {
691  LOG_INFO << " ... erased" << endm;
692  }
693 
694  int pIndex = -1;
695  auto it = mIndex2Primary.find( iNode );
696  if( it != mIndex2Primary.end() ) {
697  pIndex = it->second;
698  }
699  if( pIndex >= 0 ) {
700  StMuTrack* pTrack = ( StMuTrack* ) mMuDst->array( muPrimary )->UncheckedAt( pIndex );
701  if( pTrack ) {
702 
703  if( mDebug ) {
704  StMuETofPidTraits etofTraits = pTrack->etofPidTraits();
705 
706  LOG_INFO << "cleanUpTraits() - etof traits on primary track " << pIndex << " already exist" << endm;
707  LOG_INFO << "matchFlag: " << etofTraits.matchFlag() << " localX: " << etofTraits.localX() << " localY: " << etofTraits.localY();
708  LOG_INFO << " tof: " << etofTraits.timeOfFlight() << " pathlength: " << etofTraits.pathLength() << " beta: " << etofTraits.beta() << endm;
709 
710  if( mMuDst->etofHit( pTrack->index2ETofHit() ) ) {
711  LOG_INFO << "time: " << mMuDst->etofHit( pTrack->index2ETofHit() )->time() << endm;
712  }
713  }
714 
715  //clean up any association done before
716  pTrack->setETofPidTraits( pidETof );
717  pTrack->setIndex2ETofHit( -1 );
718 
719  if( mDebug ) {
720  LOG_INFO << " ... erased" << endm;
721  }
722  }
723  }
724  }
725 
726  size_t nHits = mMuDst->numberOfETofHit();
727  for( size_t i=0; i<nHits; i++ ) {
728  StMuETofHit* aHit = mMuDst->etofHit( i );
729  if( !aHit ) continue;
730  aHit->setIndex2Primary( -1 );
731  aHit->setIndex2Global( -1 );
732  aHit->setAssociatedTrackId( -1 );
733  }
734  }
735 }
736 
737 
738 
739 //.........................................................................
740 // A. read data from StETofHit & build vector of candidate hits
741 //
742 //---------------------------------------------------------------------------
743 void
744 StETofMatchMaker::readETofDetectorHits( eTofHitVec& detectorHitVec )
745 {
746  size_t nHits = 0;
747 
748  // event selection ... only continue with events that have eTOF hits
749  if( mIsStEventIn ) {
750  if( !mEvent->etofCollection() ) {
751  LOG_WARN << "readETofDetectorHits() - no etof collection --> nothing to do ... bye-bye" << endm;
752  return;
753  }
754  if( !mEvent->etofCollection()->hitsPresent() ) {
755  LOG_WARN << "readETofDetectorHits() - no etof hits present --> nothing to do ... bye-bye" << endm;
756  return;
757  }
758 
759  nHits = mEvent->etofCollection()->etofHits().size();
760  //LOG_INFO << " number of ETOF hits: " << nHits << endm;
761  }
762  else if( mIsMuDstIn ) {
763  if( !mMuDst->etofArray( muETofHit ) ) {
764  LOG_WARN << "readETofDetectorHits() - no digi array --> nothing to do ... bye-bye" << endm;
765  return;
766  }
767 
768  if( !mMuDst->numberOfETofHit() ) {
769  LOG_WARN << "readETofDetectorHits() - no hits present --> nothing to do ... bye-bye" << endm;
770  return;
771  }
772 
773  nHits = mMuDst->numberOfETofHit();
774  //LOG_INFO << " number of ETOF hits: " << nHits << endm;
775  }
776 
777  if( mDoQA ) {
778  mHistograms.at( "eventCounter" )->Fill( 2 );
779  }
780 
781  // fill hits into the detectorHitVec structure
782  // (depending on if StEvent or MuDst is used as input)
783  if( mIsStEventIn ) {
784  for( size_t i=0; i<nHits; i++ ) {
785  StETofHit* aHit = mEvent->etofCollection()->etofHits().at( i );
786 
787  if( !aHit ) {
788  continue;
789  }
790 
791  if( fabs(aHit->localY()) > mLocalYmax ) {//reject unphysical hits outside of detector surface from mismatches PW
792  continue;
793  }
794 
795  StructETofHit detectorHit;
796 
797  detectorHit.sector = aHit->sector();
798  /*
799  // ---------------------------------------------------
800  // rotate hits by 60 degree to avaluate random matches
801  int sec = rotateHit( aHit->sector(), 2 );
802  detectorHit.sector = sec;
803  // ---------------------------------------------------
804  */
805  detectorHit.plane = aHit->zPlane();
806  detectorHit.counter = aHit->counter();
807  detectorHit.hitTime = aHit->time();
808  detectorHit.localX = aHit->localX();
809  detectorHit.localY = aHit->localY();
810  detectorHit.tot = aHit->totalTot();
811  detectorHit.clusterSize = aHit->clusterSize();
812  detectorHit.index2ETofHit = i;
813 
814  detectorHitVec.push_back( detectorHit );
815  }
816  }
817  else {
818  for( size_t i=0; i<nHits; i++ ) {
819  StMuETofHit* aHit = mMuDst->etofHit( i );
820 
821  if( !aHit ) {
822  continue;
823  }
824 
825  if( fabs(aHit->localY()) > mLocalYmax ) {//reject unphysical hits outside of detector surface from mismatches PW
826  continue;
827  }
828 
829  StructETofHit detectorHit;
830 
831  detectorHit.sector = aHit->sector();
832  /*
833  // ---------------------------------------------------
834  // rotate hits by 60 degree to avaluate random matches
835  int sec = rotateHit( aHit->sector(), 2 );
836  detectorHit.sector = sec;
837  // ---------------------------------------------------
838  */
839  detectorHit.plane = aHit->zPlane();
840  detectorHit.counter = aHit->counter();
841  detectorHit.hitTime = aHit->time();
842  detectorHit.localX = aHit->localX();
843  detectorHit.localY = aHit->localY();
844  detectorHit.tot = aHit->totalTot();
845  detectorHit.clusterSize = aHit->clusterSize();
846  detectorHit.index2ETofHit = i;
847 
848  detectorHitVec.push_back( detectorHit );
849 
850 
851 
852 
853  }
854  }
855 
856 
857  // fill the hits in the structre with more information e.g. global position
858  // (independent from used input file format )
859  for( auto& hit : detectorHitVec ) {
860  double xl[ 3 ] = { hit.localX, hit.localY, 0 };
861 
862  int moduleId = mETofGeom->calcModuleIndex( hit.sector, hit.plane );
863  int counterId = hit.counter - 1;
864  double xg[ 3 ];
865 
866  mETofGeom->hitLocal2Master( moduleId, counterId, xl, xg );
867 
868  StThreeVectorF globalPos( xg[ 0 ], xg[ 1 ], xg[ 2 ] );
869 
870  hit.globalPos = globalPos;
871 
872  hit.strip = ( ( StETofGeomCounter* ) mETofGeom->findETofNode( moduleId, counterId ) )->findStrip( xl );
873 
874  if( mDebug ) {
875  LOG_DEBUG << "hit( " << hit.index2ETofHit << " ) at sector: " << hit.sector;
876  LOG_DEBUG << " zPlane: " << hit.plane << " counter: " << hit.counter;
877  LOG_DEBUG << " with local (X, Y): (" << xl[ 0 ] << ", " << xl[ 1 ] << ")" << endm;
878  LOG_DEBUG << "global (X, Y, Z): " << xg[ 0 ] << ", " << xg[ 1 ] << ", " << xg[ 2 ] << endm;
879  LOG_DEBUG << " strip: " << hit.strip << endm;
880  }
881 
882  // some more histogramming for QA
883  float hit_eta = globalPos.pseudoRapidity();
884  float hit_phi = globalPos.phi();
885 
886  if ( hit_phi < 0. ) hit_phi += 2. * M_PI;
887 
888  LOG_DEBUG << "global (eta, phi): " << hit_eta << ", " << hit_phi << endm;
889 
890  if( fabs( hit.localY ) < eTofConst::stripLength / 2. * 1.5 ) {
891  mHistograms.at( "eTofHits_globalXY" )->Fill( globalPos.x(), globalPos.y() );
892  }
893 
894  if( mDoQA ) {
895  if( fabs( hit.localY ) < eTofConst::stripLength / 2. * 1.5 ) {
896  mHistograms.at( "eTofHits_phi_eta" )->Fill( hit_phi, hit_eta );
897  }
898 
899  if( hit.sector == 18 || hit.sector == 24 ) {
900  mHistograms.at( "eTofHits_globalYZ" )->Fill( globalPos.y(), globalPos.z() );
901  }
902 
903  std::string histName_hit_localXY = "eTofHits_localXY_s" + std::to_string( hit.sector ) + "m" + std::to_string( hit.plane ) + "c" + std::to_string( hit.counter );
904  std::string histName_hit_globalXY = "eTofHits_globalXY_s" + std::to_string( hit.sector ) + "m" + std::to_string( hit.plane ) + "c" + std::to_string( hit.counter );
905  std::string histName_hit_eta_phi = "eTofHits_phi_eta_s" + std::to_string( hit.sector ) + "m" + std::to_string( hit.plane ) + "c" + std::to_string( hit.counter );
906 
907  mHistograms.at( histName_hit_localXY )->Fill( hit.localX, hit.localY );
908  mHistograms.at( histName_hit_globalXY )->Fill( hit.globalPos.x(), hit.globalPos.y() );
909  mHistograms.at( histName_hit_eta_phi )->Fill( hit_phi, hit_eta );
910  }
911  }
912 
913  if( mDoQA ) {
914  mHistograms.at( "detectorHitMult" )->Fill( detectorHitVec.size() );
915  if( detectorHitVec.size() > 0 ) {
916  mHistograms.at( "eventCounter" )->Fill( 3 );
917  }
918  }
919 }
920 
921 
922 //.........................................................................
923 // B. loop over global tracks & determine all track intersections with active eTof volumes
924 //
925 //---------------------------------------------------------------------------
926 void
927 StETofMatchMaker::findTrackIntersections( eTofHitVec& intersectionVec, int& nPrimaryWithIntersection )
928 {
929  size_t nNodes = 0;
930  size_t nPrimaryCrossings = 0;
931 
932  // StEvent processing
933  if( mIsStEventIn ) {
934  StSPtrVecTrackNode& nodes = mEvent->trackNodes();
935 
936  nNodes = nodes.size();
937 
938  if( mDebug ) {
939  LOG_INFO << "# tracks in the event: " << nNodes << endm;
940  }
941 
942  for( size_t iNode = 0; iNode < nNodes; iNode++ ) {
943  if( mDebug ) {
944  LOG_DEBUG << "track : " << iNode << endm;
945  }
946  StGlobalTrack* theTrack = dynamic_cast< StGlobalTrack* > ( nodes[ iNode ]->track( global ) );
947 
948  if( !validTrack( theTrack ) ) continue;
949 
950  bool isPrimary = false;
951  StPrimaryTrack* pTrack = dynamic_cast< StPrimaryTrack* > ( theTrack->node()->track( primary ) );
952  if( pTrack ) {
953  isPrimary = true;
954  if( mDebug ) {
955  LOG_DEBUG << "track : " << iNode << " is a primary track" << endm;
956  }
957  }
958 
959  StPhysicalHelixD theHelix = mOuterTrackGeometry ? theTrack->outerGeometry()->helix() : theTrack->geometry()->helix();
960 
961  int nCrossings = 0;
962 
963  extrapolateTrackToETof( intersectionVec, theHelix, iNode, nCrossings, isPrimary );
964 
965 
966  if( isPrimary ) {
967  nPrimaryCrossings += nCrossings;
968  if( nCrossings > 0 ) {
969  nPrimaryWithIntersection++;
970 
971  if( mDoQA ) {
972  int charge = pTrack->geometry()->charge();
973  float pMom = pTrack->geometry()->momentum().mag();
974 
975  mHistograms.at( "intersection_primaryTrack_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
976  if( charge > 0 ) mHistograms.at( "intersection_primaryTrackpos_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
977  else mHistograms.at( "intersection_primaryTrackneg_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
978 
979  if( pMom > 1 ) {
980  mHistograms.at( "intersection_primaryTrackMom0_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
981  if( charge > 0 ) mHistograms.at( "intersection_primaryTrackMom0pos_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
982  else mHistograms.at( "intersection_primaryTrackMom0neg_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
983  }
984  else if( pMom > 0.5 ) {
985  mHistograms.at( "intersection_primaryTrackMom1_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
986  if( charge > 0 ) mHistograms.at( "intersection_primaryTrackMom1pos_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
987  else mHistograms.at( "intersection_primaryTrackMom1neg_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
988  }
989  else {
990  mHistograms.at( "intersection_primaryTrackMom2_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
991  if( charge > 0 ) mHistograms.at( "intersection_primaryTrackMom2pos_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
992  else mHistograms.at( "intersection_primaryTrackMom2neg_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
993  }
994  }
995  }
996  }
997 
998  if( mDoQA && nCrossings > 0 ) {
999  ETofTrack eTrack( theTrack );
1000 
1001  mHistograms.at( "intersection_track_pt_eta" )->Fill( eTrack.pt, eTrack.eta );
1002  mHistograms.at( "intersection_track_pt_phi" )->Fill( eTrack.pt, eTrack.phi );
1003  mHistograms.at( "intersection_track_phi_eta" )->Fill( eTrack.phi, eTrack.eta );
1004 
1005  mHistograms.at( "intersection_track_nHitsTpc" )->Fill( eTrack.nFtPts );
1006 
1007  if( eTrack.dEdx > -999. ) mHistograms.at( "intersection_track_mom_dEdx" )->Fill( eTrack.pt * cosh( eTrack.eta ), eTrack.dEdx );
1008  if( eTrack.nSigmaPion > -999. ) mHistograms.at( "intersection_track_mom_nsigmaPi" )->Fill( eTrack.pt * cosh( eTrack.eta ), eTrack.nSigmaPion );
1009  }
1010  }
1011  } // end of StEvent processing
1012  else { // MuDst processing
1013  nNodes = mMuDst->numberOfGlobalTracks();
1014 
1015  if( mDebug ) {
1016  LOG_INFO << "# tracks in the event: " << nNodes << endm;
1017  }
1018 
1019  for( size_t iNode = 0; iNode < nNodes; iNode++ ) {
1020  if( mDebug ) {
1021  LOG_DEBUG << "track : " << iNode << endm;
1022  }
1023 
1024  StMuTrack* theTrack = mMuDst->globalTracks( iNode );
1025 
1026  if( !validTrack( theTrack ) ) continue;
1027 
1028  bool isPrimary= false;
1029 
1030  int pIndex = -1;
1031  auto it = mIndex2Primary.find( iNode );
1032  if( it != mIndex2Primary.end() ) {
1033  pIndex = it->second;
1034  }
1035  if( pIndex >= 0 ) {
1036  isPrimary = true;
1037  if( mDebug ) {
1038  LOG_DEBUG << "track : " << iNode << " is a primary track" << endm;
1039  }
1040  }
1041 
1042  StPhysicalHelixD theHelix = mOuterTrackGeometry ? theTrack->outerHelix() : theTrack->helix();
1043 
1044  int nCrossings = 0;
1045 
1046  extrapolateTrackToETof( intersectionVec, theHelix, iNode, nCrossings, isPrimary );
1047 
1048  if( isPrimary ) {
1049  nPrimaryCrossings += nCrossings;
1050  if( nCrossings > 0 ) {
1051  nPrimaryWithIntersection++;
1052 
1053  if( mDoQA ) {
1054  int charge = theTrack->primaryTrack()->charge();
1055  float pMom = theTrack->primaryTrack()->momentum().mag();
1056 
1057  mHistograms.at( "intersection_primaryTrack_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1058  if( charge > 0 ) mHistograms.at( "intersection_primaryTrackpos_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1059  else mHistograms.at( "intersection_primaryTrackneg_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1060 
1061  if( pMom > 1 ) {
1062  mHistograms.at( "intersection_primaryTrackMom0_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1063  if( charge > 0 ) mHistograms.at( "intersection_primaryTrackMom0pos_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1064  else mHistograms.at( "intersection_primaryTrackMom0neg_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1065  }
1066  else if( pMom > 0.5 ) {
1067  mHistograms.at( "intersection_primaryTrackMom1_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1068  if( charge > 0 ) mHistograms.at( "intersection_primaryTrackMom1pos_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1069  else mHistograms.at( "intersection_primaryTrackMom1neg_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1070  }
1071  else {
1072  mHistograms.at( "intersection_primaryTrackMom2_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1073  if( charge > 0 ) mHistograms.at( "intersection_primaryTrackMom2pos_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1074  else mHistograms.at( "intersection_primaryTrackMom2neg_globalXY" )->Fill( intersectionVec.back().globalPos.x() , intersectionVec.back().globalPos.y() );
1075  }
1076  }
1077  }
1078  }
1079 
1080  if( mDoQA && nCrossings > 0 ) {
1081  ETofTrack eTrack( theTrack );
1082 
1083  mHistograms.at( "intersection_track_pt_eta" )->Fill( eTrack.pt, eTrack.eta );
1084  mHistograms.at( "intersection_track_pt_phi" )->Fill( eTrack.pt, eTrack.phi );
1085  mHistograms.at( "intersection_track_phi_eta" )->Fill( eTrack.phi, eTrack.eta );
1086 
1087  mHistograms.at( "intersection_track_nHitsTpc" )->Fill( eTrack.nFtPts );
1088 
1089  if( eTrack.dEdx > 0. ) mHistograms.at( "intersection_track_mom_dEdx" )->Fill( eTrack.pt * cosh( eTrack.eta ), eTrack.dEdx );
1090  if( eTrack.nSigmaPion > -999. ) mHistograms.at( "intersection_track_mom_nsigmaPi" )->Fill( eTrack.pt * cosh( eTrack.eta ), eTrack.nSigmaPion );
1091  }
1092  }
1093  } // end of MuDst processing
1094 
1095  //LOG_INFO << "# tracks in the event: " << nNodes << " ... out of which " << intersectionVec.size() << " intersect with eTOF" << endm;
1096 
1097  if( mDoQA ) {
1098  mHistograms.at( "intersectionMult" )->Fill( intersectionVec.size() );
1099  mHistograms.at( "intersectionMult_primary" )->Fill( nPrimaryCrossings );
1100 
1101  if( intersectionVec.size() > 0 ) {
1102  mHistograms.at( "eventCounter" )->Fill( 4 );
1103  }
1104  }
1105 }
1106 
1107 
1108 //---
1110 bool
1111 StETofMatchMaker::validTrack( const StTrack* track )
1112 {
1113  if( track ) {
1114  return validTrack( ETofTrack( track ) );
1115  }
1116  else {
1117  return false;
1118  }
1119 }
1120 
1121 //---
1123 bool
1124 StETofMatchMaker::validTrack( const StMuTrack* track )
1125 {
1126  if( track ) {
1127  return validTrack( ETofTrack( track ) );
1128  }
1129  else {
1130  return false;
1131  }
1132 }
1133 
1134 //---
1136 bool
1137 StETofMatchMaker::validTrack( const ETofTrack& track )
1138 {
1139  double ratioFitToPoss = 1. * track.nFtPts / track.nHitsPoss;
1140 
1141  if( mDoQA ) {
1142  mHistograms.at( "track_phi_eta" ) ->Fill( track.phi, track.eta );
1143  mHistograms.at( "track_phi_pt" ) ->Fill( track.phi, track.pt );
1144  mHistograms.at( "nHits" ) ->Fill( track.nFtPts );
1145  mHistograms.at( "track_pt_nHits" )->Fill( track.pt, track.nFtPts );
1146  }
1147 
1148  // kick out tracks that will not point to the eTOF
1149  // TODO: more carefull eta acceptance cut in the future (performance vs. efficientcy)
1150  if( track.eta > etofProjection::minEtaCut ) return false;
1151 
1152  if( mDoQA && track.eta > -1.65 && track.eta < -1.05 ) {
1153  mHistograms.at( "nHits_etofregion" )->Fill( track.nFtPts );
1154  }
1155 
1156  // kick out low quality tracks
1157  //if( track.flag <= flagMinCut ) return false;
1158  //if( track.flag >= flagMaxCut ) return false;
1159  if( track.nFtPts < mTrackCuts.at( 0 ) ) return false;
1160  if( ratioFitToPoss < mTrackCuts.at( 1 ) ) return false;
1161  if( track.pt < mTrackCuts.at( 2 ) ) return false;
1162  //if( track.mom < mTrackCuts.at( 2 ) ) return false;
1163 
1164  if( mDebug ) {
1165  LOG_DEBUG << "valid track for extrapolation to eTOF with nHitsFit: " << track.nFtPts;
1166  LOG_DEBUG << " mom: " << track.mom << " pt: " << track.pt << " phi: " << track.phi << " eta: " << track.eta << endm;
1167  }
1168 
1169  return true;
1170 }
1171 
1172 //---
1173 void
1174 StETofMatchMaker::extrapolateTrackToETof( eTofHitVec& intersectionVec, const StPhysicalHelixD& theHelix, const int& iNode, int& nCrossings, bool isPrimary )
1175 {
1176  // first project helix to the middle eTof plane to get the sector(s) of possible intersections
1177  StThreeVectorD projection = mETofGeom->helixCrossPlane( theHelix, eTofConst::zplanes[ 1 ] );
1178 
1179  // get rid of tracks that do not project to the rough eta region of the eTof
1180  float projEta = projection.pseudoRapidity();
1181 
1182  if( projEta < etofProjection::maxEtaProjCut ) return;
1183  if( projEta > etofProjection::minEtaProjCut ) return;
1184 
1185  float projPhi = projection.phi();
1186  if ( projPhi < 0. ) projPhi += 2. * M_PI;
1187 
1188  if( mDoQA ) {
1189  // histogramming for QA
1190  mHistograms.at( "trackProj_globalXY" )->Fill( projection.x(), projection.y() );
1191  mHistograms.at( "trackProj_phi_eta" )->Fill( projPhi, projEta );
1192  }
1193 
1194  vector< int > idVec;
1195  vector< StThreeVectorD > globalVec;
1196  vector< StThreeVectorD > localVec;
1197  vector< double > thetaVec;
1198  vector< double > pathLenVec;
1199 
1200  // look for intersections of the extrapolated helix with the active eTof volumes
1201  mETofGeom->helixCrossCounter( theHelix, idVec, globalVec, localVec, thetaVec, pathLenVec );
1202 
1203  nCrossings = idVec.size();
1204 
1205  // loop backwards through the vectors, so that one can remove the
1206  // current entry if it doesn't pass e.g. a local Y cut in the future
1207  for( int i = nCrossings-1; i >= 0; i-- ) {
1208  if( mDebug ) {
1209  LOG_INFO << " ------> crossing with volume index: " << idVec.at( i ) << endm;
1210  LOG_INFO << "track intersection: " << globalVec.at( i ).x() << ", " << globalVec.at( i ).y() << ", " << globalVec.at( i ).z() << endm;
1211  LOG_INFO << "local coordinates: " << localVec.at( i ).x() << ", " << localVec.at( i ).y() << ", " << localVec.at( i ).z() << endm;
1212  }
1213 
1214  StructETofHit intersect;
1215 
1216  mETofGeom->decodeVolumeIndex( idVec.at( i ), intersect.sector, intersect.plane, intersect.counter, intersect.strip );
1217 
1218  intersect.localX = localVec.at( i ).x();
1219  intersect.localY = localVec.at( i ).y();
1220  intersect.globalPos = globalVec.at( i );
1221  intersect.trackId = iNode;
1222  intersect.theta = thetaVec.at( i );
1223  intersect.pathLength = pathLenVec.at( i );
1224  intersect.isPrimary = isPrimary;
1225 
1226  // fill intersection into storage vector
1227  intersectionVec.push_back( intersect );
1228 
1229 
1230  if( mDoQA ) {
1231  float intersectPhi = intersect.globalPos.phi();
1232  if( intersectPhi < 0. ) intersectPhi += 2. * M_PI;
1233 
1234  mHistograms.at( "intersection_globalXY" )->Fill( intersect.globalPos.x(), intersect.globalPos.y() );
1235  mHistograms.at( "intersection_phi_eta" )->Fill( intersectPhi, intersect.globalPos.pseudoRapidity() );
1236  }
1237  }
1238 
1239  if( mDoQA ) {
1240  mHistograms.at( "intersection_perTrack" )->Fill( idVec.size() );
1241  }
1242 }
1243 
1244 
1245 //.........................................................................
1246 // C. match detector hits to track intersections
1247 //
1248 //---------------------------------------------------------------------------
1249 void
1250 StETofMatchMaker::matchETofHits( eTofHitVec& detectorHitVec, eTofHitVec& intersectionVec, eTofHitVec& matchCandVec )
1251 {
1252  for( auto detHitIter = detectorHitVec.begin(); detHitIter != detectorHitVec.end(); detHitIter++ ) {
1253  for( auto interIter = intersectionVec.begin(); interIter != intersectionVec.end(); interIter++ ) {
1254 
1255  //fill correlation histograms
1256  int sector = detHitIter->sector;
1257  int plane = detHitIter->plane;
1258 
1259  int moduleId = ( sector - eTofConst::sectorStart ) * eTofConst::nPlanes + plane - eTofConst::zPlaneStart;
1260 
1261  if( mDoQA ) {
1262  int detHitIndex = ( detHitIter->counter - eTofConst::counterStart ) * eTofConst::nStrips + detHitIter->strip - eTofConst::stripStart;
1263  int intersecIndex = ( interIter->counter - eTofConst::counterStart ) * eTofConst::nStrips + interIter->strip - eTofConst::stripStart;
1264 
1265  mHistograms.at( "detHitvsInter_strip_s" + std::to_string( sector ) + "m" + std::to_string( plane ) )->Fill( detHitIndex, intersecIndex );
1266 
1267  mHistograms.at( "detHitvsInter_X" )->Fill( detHitIter->globalPos.x(), interIter->globalPos.x() );
1268  mHistograms.at( "detHitvsInter_Y" )->Fill( detHitIter->globalPos.y(), interIter->globalPos.y() );
1269 
1270  mHistograms.at( "moduleIndex_deltaX" )->Fill( moduleId, detHitIter->localX - interIter->localX );
1271  mHistograms.at( "moduleIndex_deltaY" )->Fill( moduleId, detHitIter->localY - interIter->localY );
1272 
1273  mHistograms.at( "detHitvsInter_localX" )->Fill(detHitIter->localX, interIter->localX );
1274  mHistograms.at( "detHitvsInter_localY" )->Fill(detHitIter->localY, interIter->localY );
1275  }
1276 
1277 
1278  // store match candidates
1279  bool isMatch = false;
1280 
1281  // deltaX, deltaY (subtract offset until alignment is done properly)
1282  float deltaX = detHitIter->localX - interIter->localX;
1283  float deltaY = detHitIter->localY - interIter->localY;
1284  double tstart = startTimeBTof(); //no eToF start time available here!
1285  double deltaT = detHitIter->hitTime - tstart; //basic cut to reject hits far of in time
1286 
1287  int counterIndex = ( detHitIter->sector - eTofConst::sectorStart ) * eTofConst::nPlanes * eTofConst::nCounters
1288  + ( detHitIter->plane - eTofConst::zPlaneStart ) * eTofConst::nCounters
1289  + ( detHitIter->counter - eTofConst::counterStart );
1290 
1291  deltaX -= etofProjection::deltaXoffset[ counterIndex ];
1292  deltaY -= etofProjection::deltaYoffset[ counterIndex ];
1293 
1294  bool corrTime=false; // for single sided hit time corr
1295 
1296  if( detHitIter->sector == interIter->sector ) {
1297  if( detHitIter->plane == interIter->plane ) {
1298  if( detHitIter->counter == interIter->counter ) {
1299 
1300  if(detHitIter->clusterSize < 999){
1301 
1302  // if( fabs( deltaX ) < mMatchDistX ) {
1303  // if( fabs( deltaY ) < mMatchDistY ) {
1304 
1305  if( ( ( (deltaY*deltaY) / (mMatchDistY*mMatchDistY) ) + ( (deltaX*deltaX) / (mMatchDistX*mMatchDistX) ) ) < 2. ) {
1306  if( fabs( deltaT ) < mMatchDistT ) {
1307  isMatch = true;
1308  }
1309  }
1310  }else{
1311 
1312  float mMatchDistYSingleSided = 15;
1313 
1314 
1315 
1316  if( fabs( deltaX ) < mMatchDistX ) {
1317  if( fabs( deltaY ) < mMatchDistYSingleSided ) {
1318  if( fabs( deltaT ) < mMatchDistT ) {
1319 
1320 
1321  isMatch = true;
1322  deltaY = 27; // keep SHs out of NHs way while sorting
1323  corrTime = true;
1324 
1325  }
1326  }
1327  }
1328  }
1329  }
1330  }
1331  }
1332 
1333  if( isMatch ) {
1334  StructETofHit matchCand;
1335 
1336  matchCand.sector = detHitIter->sector;
1337  matchCand.plane = detHitIter->plane;
1338  matchCand.counter = detHitIter->counter;
1339  matchCand.strip = detHitIter->strip;
1340  matchCand.localX = detHitIter->localX;
1341  matchCand.localY = detHitIter->localY;
1342  matchCand.hitTime = detHitIter->hitTime;
1343  matchCand.tot = detHitIter->tot;
1344  matchCand.clusterSize = detHitIter->clusterSize;
1345  matchCand.index2ETofHit = detHitIter->index2ETofHit;
1346  matchCand.IdTruthHit = detHitIter->IdTruth;
1347 
1348  matchCand.globalPos = interIter->globalPos;
1349  matchCand.trackId = interIter->trackId;
1350  matchCand.theta = interIter->theta;
1351  matchCand.pathLength = interIter->pathLength;
1352  matchCand.isPrimary = interIter->isPrimary;
1353  matchCand.IdTruth = interIter->IdTruth;
1354 
1355  matchCand.matchFlag = 0;
1356  matchCand.deltaX = deltaX;
1357  matchCand.deltaY = deltaY;
1358 
1359  matchCand.tof = -999.;
1360  matchCand.beta = -999.;
1361 
1362  // correct single sided matches
1363  if(corrTime){
1364  matchCand.localY = interIter->localY;
1365 
1366  // if side A
1367  double corr ;
1368  float tcorr = 0;
1369  if(sector == 15 || sector == 17 || sector == 21 || sector == 22 ){
1370  tcorr = 16.49;
1371  }else{
1372  tcorr = 18.23;
1373  }
1374  if(detHitIter->localY < 0){
1375  matchCand.hitTime = detHitIter->hitTime - (((13.5 + interIter->localY ) / tcorr )) + (13.5/tcorr);
1376  // matchCand.totDiff = 1;
1377  corr = (((13.5 - interIter->localY ) / tcorr ));
1378  // if side B
1379  }else{
1380  matchCand.hitTime = detHitIter->hitTime - (((13.5 - interIter->localY ) / tcorr )) + (13.5/tcorr);
1381  // matchCand.totDiff = -1;
1382  corr = (((13.5 + interIter->localY ) / tcorr ));
1383  }
1384 
1385  matchCand.totDiff = matchCand.totDiff * corr;
1386 
1387  // cout << "interIter->localY " << interIter->localY<< endl;
1388  // cout << "corr " << corr << endl;
1389  // cin.get();
1390 
1391  }
1392 
1393  matchCandVec.push_back( matchCand );
1394 
1395  if( mDebug ) {
1396  LOG_INFO << " * * FOUND MATCH CAND : " << matchCand.sector << " " << matchCand.plane << " " << matchCand.counter;
1397  LOG_INFO << " with (deltaX, deltaY) = (" << deltaX << ", " << deltaY << ")" << endm;
1398  }
1399 
1400  mHistograms.at( "matchCand_globalXY" )->Fill( matchCand.globalPos.x(), matchCand.globalPos.y() );
1401 
1402  if( mDoQA ) {
1403  float matchCandPhi = matchCand.globalPos.phi();
1404  if ( matchCandPhi < 0. ) matchCandPhi += 2. * M_PI;
1405 
1406  mHistograms.at( "matchCand_phi_eta" )->Fill( matchCandPhi, matchCand.globalPos.pseudoRapidity() );
1407 
1408  mHistograms.at( "matchCand_deltaX" )->Fill( deltaX );
1409  mHistograms.at( "matchCand_deltaY" )->Fill( deltaY );
1410 
1411  std::string histName_deltaX = "matchCand_deltaX_s" + std::to_string( matchCand.sector ) + "m" + std::to_string( matchCand.plane ) + "c" + std::to_string( matchCand.counter );
1412  std::string histName_deltaY = "matchCand_deltaY_s" + std::to_string( matchCand.sector ) + "m" + std::to_string( matchCand.plane ) + "c" + std::to_string( matchCand.counter );
1413 
1414  mHistograms.at( histName_deltaX )->Fill( deltaX );
1415  mHistograms.at( histName_deltaY )->Fill( deltaY );
1416 
1417 
1418  // get track corresponding to the trackId stored for matchCand
1419  if( mIsStEventIn ) {
1420  StSPtrVecTrackNode& nodes = mEvent->trackNodes();
1421  StGlobalTrack *matchCandTrack = dynamic_cast< StGlobalTrack* > ( nodes[ matchCand.trackId ]->track( global ) );
1422  if( matchCandTrack ) {
1423  int nFitPts = matchCandTrack->fitTraits().numberOfFitPoints( kTpcId );
1424 
1425  mHistograms.at( "matchCand_deltaX_nHitsTpc" )->Fill( nFitPts, deltaX );
1426  mHistograms.at( "matchCand_deltaY_nHitsTpc" )->Fill( nFitPts, deltaY );
1427  }
1428  }
1429  else {
1430  StMuTrack* matchCandTrack = mMuDst->globalTracks( matchCand.trackId );
1431  if( matchCandTrack ) {
1432  int nFitPts = matchCandTrack->nHitsFit( kTpcId );
1433 
1434  mHistograms.at( "matchCand_deltaX_nHitsTpc" )->Fill( nFitPts, deltaX );
1435  mHistograms.at( "matchCand_deltaY_nHitsTpc" )->Fill( nFitPts, deltaY );
1436  }
1437  }
1438  }
1439  }
1440 
1441  }
1442  }
1443 
1444  if( mDoQA ) {
1445  mHistograms.at( "matchCandMult" )->Fill( matchCandVec.size() );
1446 
1447  if( matchCandVec.size() > 0 ) {
1448  mHistograms.at( "eventCounter" )->Fill( 5 );
1449  }
1450  }
1451 }
1452 
1453 
1454 //.........................................................................
1455 // D. sort matchCand vector and deal with (discard) hits matched by multiple tracks
1456 //
1457 //---------------------------------------------------------------------------
1458 void
1459 StETofMatchMaker::sortSingleMultipleHits( eTofHitVec& matchCandVec, eTofHitVec& singleTrackMatchVec, std::vector< eTofHitVec >& multiTrackMatchVec )
1460 {
1461  int nSingleTrackMatch = 0;
1462  int nMultiTrackMatch = 0;
1463 
1464 
1465  // define temporary vectors for iterating through matchCandVec
1466  eTofHitVec tempVec = matchCandVec;
1467  eTofHitVec erasedVec = tempVec;
1468 
1469  while( tempVec.size() != 0 ) {
1470  int nTracks = 0;
1471 
1472  // define temporary storage vectors
1473  eTofHitVec candVec;
1474 
1475  eTofHitVecIter tempIter = tempVec.begin();
1476  eTofHitVecIter erasedIter = erasedVec.begin();
1477 
1478  while( erasedIter != erasedVec.end() ) {
1479 
1480  if( tempIter->index2ETofHit == erasedIter->index2ETofHit ) {
1481  nTracks++;
1482  candVec.push_back( *erasedIter );
1483 
1484  erasedVec.erase( erasedIter );
1485  erasedIter--;
1486  }
1487  erasedIter++;
1488  }
1489  // reset tempVec for next iteration in the while-loop
1490  tempVec = erasedVec;
1491 
1492  StructETofHit cand = candVec.front();
1493 
1494  if( mDebug ) {
1495  LOG_INFO << "matchCand at sector " << cand.sector << " plane " << cand.plane << " counter " << cand.counter << " with local (X,Y) = (" << cand.localX << "," << cand.localY << ")";
1496  LOG_INFO << " has " << nTracks << " track(s) pointing to it:" << endm;
1497  }
1498 
1499  if( nTracks == 1 ) {
1500  nSingleTrackMatch++;
1501  singleTrackMatchVec.push_back( cand );
1502 
1503  if( mDebug ) {
1504  LOG_INFO << "track id: " << cand.trackId << " and matching distance " << cand.deltaX << " " << cand.deltaY << endm;
1505  }
1506  }
1507  else if( nTracks > 1 ) {
1508  // for multiple tracks pointing to the same detector hit: either discard or take "most likely" / "best" match candidate
1509  // for now: take match with smallest deltaR
1510  nMultiTrackMatch++;
1511 
1512  multiTrackMatchVec.push_back( candVec );
1513 
1514  float bestResidual = 999.;
1515  int bestIndex = -1;
1516 
1517  for( size_t i=0; i<candVec.size(); i++ ) {
1518  if( mDebug ) {
1519  LOG_INFO << "track id: " << candVec.at( i ).trackId << " and matching distance " << candVec.at( i ).deltaX << " " << candVec.at( i ).deltaY << endm;
1520  }
1521  float residual = pow( candVec.at( i ).deltaX, 2 ) + pow( candVec.at( i ).deltaY, 2 );
1522 
1523  if( residual < bestResidual ) {
1524  bestResidual = residual;
1525  bestIndex = i;
1526  }
1527  }
1528 
1529  if( bestIndex > -1 ) {
1530  singleTrackMatchVec.push_back( candVec.at( bestIndex ) );
1531  if( mDebug ) {
1532  LOG_INFO << "best candidate has track id: " << candVec.at( bestIndex ).trackId << endm;
1533  }
1534  }
1535 
1536  if( mDebug ) {
1537  for( const auto& c: candVec ) {
1538  LOG_INFO << "track id: " << c.trackId << " and matching distance " << c.deltaX << " " << c.deltaY << endm;
1539  }
1540  }
1541  }
1542 
1543  if( mDoQA ) {
1544  mHistograms.at( "trackMatchMultPerDetectorHit" )->Fill( nTracks );
1545  }
1546  }
1547 
1548  //LOG_INFO << "nSingleTrackMatch: " << nSingleTrackMatch << " ... nMultiTrackMatch: " << nMultiTrackMatch << endm;
1549 
1550  if( mDoQA ) {
1551  mHistograms.at( "singleTrackMatchMult" )->Fill( singleTrackMatchVec.size() );
1552 
1553  if( singleTrackMatchVec.size() > 0 ) {
1554  mHistograms.at( "eventCounter" )->Fill( 6 );
1555  }
1556  }
1557 }
1558 
1559 
1560 //.........................................................................
1561 // E. sort singleTrackMatchVector for multiple hits associated to single tracks and determine the best match
1562 // also set the match flag ( TODO: set it more sophisticated when dealing with multi-hits )
1563 //
1564 //---------------------------------------------------------------------------
1565 void
1566 StETofMatchMaker::finalizeMatching( eTofHitVec& singleTrackMatchVec, eTofHitVec& finalMatchVec )
1567 {
1568  // setup temporary vectors for iterating trough singleTrackMatchVec
1569  eTofHitVec tempVec = singleTrackMatchVec;
1570  eTofHitVec erasedVec = tempVec;
1571 
1572  while( tempVec.size() != 0 ) {
1573  int nHits = 0;
1574 
1575  // define temporary storage vectors
1576  eTofHitVec candVec;
1577 
1578  eTofHitVecIter tempIter = tempVec.begin();
1579  eTofHitVecIter erasedIter = erasedVec.begin();
1580 
1581  while( erasedIter != erasedVec.end() ) {
1582 
1583  if( tempIter->trackId == erasedIter->trackId ) {
1584  nHits++;
1585  candVec.push_back( *erasedIter );
1586 
1587  erasedVec.erase( erasedIter );
1588  erasedIter--;
1589  }
1590  erasedIter++;
1591  }
1592  // reset tempVec for next iteration in the while-loop
1593  tempVec = erasedVec;
1594 
1595  // for a track matched to only one eTof hit -> copy match candidate to finalMatchVec
1596  if( nHits == 1 ) {
1597  candVec.front().matchFlag = 1;
1598  finalMatchVec.push_back( candVec.front() );
1599 
1600  if( mDebug ) {
1601  LOG_INFO << "one-to-one match (track id: " << candVec.front().trackId << ") -> pushed into final match vector" << endm;
1602  }
1603  }
1604  // for a track matched to many eTof hits -> only take the most likely / best match
1605  else if ( nHits > 1 ) {
1606  if( mDebug ) {
1607  LOG_INFO << "one-to-many match -> needs further treatment" << endm;
1608  }
1609 
1610  // for the moment sort on distance between track intersection and eTof detector hit:
1611  double bestMatchDist = pow( candVec.front().deltaX, 2 ) + pow( candVec.front().deltaY, 2 );
1612  StructETofHit bestCand = candVec.front();
1613 
1614  bool isOverlapHit = false;
1615 
1616  for( const auto& c: candVec ) {
1617  double candMatchDist = pow( c.deltaX, 2 ) + pow( c.deltaY, 2 );
1618 
1619  if( mDebug ) {
1620  LOG_INFO << "candidate match distance: " << sqrt( candMatchDist ) << endm;
1621  }
1622 
1623  if( candMatchDist < bestMatchDist ) {
1624  bestMatchDist = candMatchDist;
1625  bestCand = c;
1626 
1627  if( mDebug ) {
1628  LOG_INFO << " --> new best match candidate" << endm;
1629  }
1630  }
1631 
1632  if( ( bestCand.sector * 100 + bestCand.plane * 10 + bestCand.counter ) != ( c.sector * 100 + c.plane * 10 + c.counter ) ) {
1633  isOverlapHit = true;
1634  }
1635  }
1636 
1637  if( isOverlapHit ) {
1638  bestCand.matchFlag = 3;
1639 
1640  if( mDoQA ) {
1641  mHistograms.at( "overlapHit_globalXY" )->Fill( bestCand.globalPos.x(), bestCand.globalPos.y() );
1642  }
1643  }
1644  else {
1645  bestCand.matchFlag = 2;
1646  }
1647 
1648  finalMatchVec.push_back( bestCand );
1649 
1650  if( mDebug ) {
1651  LOG_INFO << "one-to-many match resolved (track id: " << bestCand.trackId << ", cand deltaX: " << bestCand.deltaX << ") -> pushed into final match vector" << endm;
1652  }
1653  }
1654 
1655  if( mDoQA ) {
1656  mHistograms.at( "hitMultPerTrack" )->Fill( nHits );
1657  }
1658  }
1659 
1660  if( mDoQA ) {
1661  if( mIsStEventIn ) {
1662  StSPtrVecTrackNode& nodes = mEvent->trackNodes();
1663  for( const auto& v: finalMatchVec ) {
1664  StGlobalTrack* track = dynamic_cast< StGlobalTrack* > ( nodes[ v.trackId ]->track( global ) );
1665 
1666  mHistograms.at( "finalMatch_pt" )->Fill( track->geometry()->momentum().perp() );
1667  }
1668  }
1669  else {
1670  for( const auto& v: finalMatchVec ) {
1671  StMuTrack* track = mMuDst->globalTracks( v.trackId );
1672 
1673  mHistograms.at( "finalMatch_pt" )->Fill( track->momentum().perp() );
1674  }
1675  }
1676 
1677  mHistograms.at( "finalMatchMult" )->Fill( finalMatchVec.size() );
1678 
1679  if( finalMatchVec.size() > 0 ) {
1680  mHistograms.at( "eventCounter" )->Fill( 7 );
1681  }
1682  }
1683 }
1684 
1685 
1686 //.........................................................................
1687 // F. fill ETofPidTraits for global and primary tracks and assign associated track to hits
1688 //
1689 //---------------------------------------------------------------------------
1690 void
1691 StETofMatchMaker::fillPidTraits( eTofHitVec& finalMatchVec )
1692 {
1693  size_t nFinalMatchesGlobal = 0;
1694  size_t nFinalMatchesPrimary = 0;
1695 
1696  if( mIsStEventIn ) {
1697  //get track & hit collection
1698  StSPtrVecETofHit& etofHits = mEvent->etofCollection()->etofHits();
1699  StSPtrVecTrackNode& nodes = mEvent->trackNodes();
1700 
1701  for( const auto& match : finalMatchVec ) {
1702  // get global track
1703  StGlobalTrack* globalTrack = dynamic_cast< StGlobalTrack* > ( nodes[ match.trackId ]->track( global ) );
1704  if( !globalTrack ) {
1705  LOG_WARN << "fillPidTraits() - global track does not exist" << endm;
1706  continue;
1707  }
1708 
1709  // fill association in ETofCollection
1710  StETofHit* hit = etofHits[ match.index2ETofHit ];
1711  if( !hit ) {
1712  LOG_WARN << "fillPidTraits() - eTof hit does not exist" << endm;
1713  continue;
1714  }
1715  hit->setAssociatedTrack( globalTrack );
1716 
1717  nFinalMatchesGlobal++;
1718 
1719  // fill the matching data into StETofPidTraits
1720  StETofPidTraits* pidTraits = new StETofPidTraits();
1721 
1722  pidTraits->setETofHit( hit );
1723  pidTraits->setMatchFlag( match.matchFlag );
1724  pidTraits->setLocalX( match.localX );
1725  pidTraits->setLocalY( match.localY );
1726  pidTraits->setThetaLocal( match.theta );
1727  pidTraits->setDeltaX( match.deltaX );
1728  pidTraits->setDeltaY( match.deltaY );
1729  pidTraits->setPosition( match.globalPos );
1730 
1731  pidTraits->setPathLength( match.pathLength );
1732 
1733  pidTraits->setTimeOfFlight( -999. );
1734  pidTraits->setBeta( -999. );
1735 
1736  globalTrack->addPidTraits( pidTraits );
1737 
1738  // get primary track if exists and fill matching data into StETofPidTraits
1739  StPrimaryTrack *pTrack = dynamic_cast< StPrimaryTrack* > ( nodes[ match.trackId ]->track( primary ) );
1740  if( pTrack ) {
1741  nFinalMatchesPrimary++;
1742 
1743  if( mDoQA ) {
1744  mHistograms.at( "finalMatchPrimary_globalXY" )->Fill( match.globalPos.x(), match.globalPos.y() );
1745 
1746  float mom = pTrack->geometry()->momentum().mag();
1747  if( mom > 1 ) {
1748  mHistograms.at( "finalMatchPrimaryMom0_globalXY" )->Fill( match.globalPos.x() , match.globalPos.y() );
1749  }
1750  else if( mom > 0.5 ) {
1751  mHistograms.at( "finalMatchPrimaryMom1_globalXY" )->Fill( match.globalPos.x() , match.globalPos.y() );
1752  }
1753  else {
1754  mHistograms.at( "finalMatchPrimaryMom2_globalXY" )->Fill( match.globalPos.x() , match.globalPos.y() );
1755  }
1756  }
1757 
1758  StETofPidTraits* ppidTraits = new StETofPidTraits();
1759 
1760  ppidTraits->setETofHit( hit );
1761  ppidTraits->setMatchFlag( match.matchFlag );
1762  ppidTraits->setLocalX( match.localX );
1763  ppidTraits->setLocalY( match.localY );
1764  ppidTraits->setThetaLocal( match.theta );
1765  ppidTraits->setDeltaX( match.deltaX );
1766  ppidTraits->setDeltaY( match.deltaY );
1767  ppidTraits->setPosition( match.globalPos );
1768 
1769  ppidTraits->setPathLength( match.pathLength );
1770 
1771  ppidTraits->setTimeOfFlight( -999. );
1772  ppidTraits->setBeta( -999. );
1773 
1774  pTrack->addPidTraits( ppidTraits );
1775  }
1776  }
1777  }
1778  else {
1779  for( const auto& match : finalMatchVec ) {
1780  // get global track
1781  StMuTrack* gTrack = mMuDst->globalTracks( match.trackId );
1782  if( !gTrack ) {
1783  LOG_WARN << "fillPidTraits() - global track does not exist" << endm;
1784  continue;
1785  }
1786 
1787  // fill association to hit
1788  StMuETofHit* hit = mMuDst->etofHit( match.index2ETofHit );
1789  if( !hit ) {
1790  LOG_WARN << "fillPidTraits() - eTof hit does not exist" << endm;
1791  continue;
1792  }
1793  hit->setAssociatedTrackId( gTrack->id() );
1794  hit->setIndex2Global( match.trackId );
1795  gTrack->setIndex2ETofHit( match.index2ETofHit );
1796 
1797  nFinalMatchesGlobal++;
1798 
1799  // fill the matching data into StETofPidTraits
1800  StMuETofPidTraits pidTraits = gTrack->etofPidTraits();
1801  pidTraits.setMatchFlag( match.matchFlag );
1802  pidTraits.setLocalX( match.localX );
1803  pidTraits.setLocalY( match.localY );
1804  pidTraits.setThetaLocal( match.theta );
1805  pidTraits.setDeltaX( match.deltaX );
1806  pidTraits.setDeltaY( match.deltaY );
1807  pidTraits.setPosition( match.globalPos );
1808 
1809  pidTraits.setPathLength( match.pathLength );
1810 
1811  pidTraits.setTimeOfFlight( -999. );
1812  pidTraits.setBeta( -999. );
1813 
1814  gTrack->setETofPidTraits ( pidTraits );
1815 
1816  int pIndex = -1;
1817  auto it = mIndex2Primary.find( match.trackId );
1818  if( it != mIndex2Primary.end() ) {
1819  pIndex = it->second;
1820  }
1821  if( pIndex < 0 ) continue;
1822 
1823  StMuTrack* pTrack = ( StMuTrack* ) mMuDst->array( muPrimary )->UncheckedAt( pIndex );
1824  if( pTrack ) {
1825  nFinalMatchesPrimary++;
1826 
1827  if( mDoQA ) {
1828  mHistograms.at( "finalMatchPrimary_globalXY" )->Fill( match.globalPos.x(), match.globalPos.y() );
1829 
1830  float mom = pTrack->momentum().mag();
1831  if( mom > 1 ) {
1832  mHistograms.at( "finalMatchPrimaryMom0_globalXY" )->Fill( match.globalPos.x() , match.globalPos.y() );
1833  }
1834  else if( mom > 0.5 ) {
1835  mHistograms.at( "finalMatchPrimaryMom1_globalXY" )->Fill( match.globalPos.x() , match.globalPos.y() );
1836  }
1837  else {
1838  mHistograms.at( "finalMatchPrimaryMom2_globalXY" )->Fill( match.globalPos.x() , match.globalPos.y() );
1839  }
1840  }
1841 
1842 
1843  hit->setIndex2Primary( pIndex );
1844  pTrack->setIndex2ETofHit( match.index2ETofHit );
1845 
1846  StMuETofPidTraits ppidTraits = pTrack->etofPidTraits();
1847  ppidTraits.setMatchFlag( match.matchFlag );
1848  ppidTraits.setLocalX( match.localX );
1849  ppidTraits.setLocalY( match.localY );
1850  ppidTraits.setThetaLocal( match.theta );
1851  ppidTraits.setDeltaX( match.deltaX );
1852  ppidTraits.setDeltaY( match.deltaY );
1853  ppidTraits.setPosition( match.globalPos );
1854 
1855  ppidTraits.setPathLength( match.pathLength );
1856 
1857  ppidTraits.setTimeOfFlight( -999. );
1858  ppidTraits.setBeta( -999. );
1859 
1860  pTrack->setETofPidTraits ( ppidTraits );
1861  }
1862  }
1863  }
1864 
1865  if( mDoQA ) {
1866  mHistograms.at( "finalMatchMultGlobal" )->Fill( nFinalMatchesGlobal );
1867  mHistograms.at( "finalMatchMultPrimary" )->Fill( nFinalMatchesPrimary );
1868  }
1869 }
1870 
1871 
1872 //.........................................................................
1873 // G. calculate pid variables for primary tracks and update PidTraits
1874 //
1875 //---------------------------------------------------------------------------
1876 void
1877 StETofMatchMaker::calculatePidVariables( eTofHitVec& finalMatchVec, int& nPrimaryWithPid )
1878 {
1879  // get start time
1880  double tstart = startTime( finalMatchVec );
1881 
1882  if( fabs( tstart + 9999. ) < 1.e-5 ) {
1883  LOG_WARN << "calculatePidVariables() -- no valid start time available ... skip filling pidTraits with more information" << endm;
1884  nPrimaryWithPid = -1;
1885  return;
1886  }
1887 
1888  if( mIsStEventIn ) { // StEvent processing ...
1889  // assign pathlength, time-of-flight, beta ... to the match candidates
1890  for( auto& matchCand : finalMatchVec ) {
1891 
1892  StETofHit* aHit = dynamic_cast< StETofHit* > ( mEvent->etofCollection()->etofHits().at( matchCand.index2ETofHit ) );
1893  if( !aHit ) {
1894  LOG_ERROR << "calculatePidVariables() - no hit" << endm;
1895  continue;
1896  }
1897 
1898  // global track
1899  StTrack* gTrack = aHit->associatedTrack();
1900  if( !gTrack ) {
1901  LOG_ERROR << "calculatePidVariables() - no associated track" << endm;
1902  continue;
1903  }
1904 
1905 
1906  StSPtrVecTrackPidTraits& traits = gTrack->pidTraits();
1907  StETofPidTraits* pidTraits = 0;
1908  for( size_t i=0; i<traits.size(); i++ ) {
1909  if( traits[ i ]->detector() == kETofId ) {
1910  pidTraits = dynamic_cast< StETofPidTraits* > ( traits[ i ] );
1911  break;
1912  }
1913  }
1914  if( !pidTraits ) continue;
1915 
1916 
1917  double tof = timeOfFlight( tstart, aHit->time() );
1918 
1919  // set time-of-flight
1920  matchCand.tof = tof;
1921 
1922  pidTraits->setTimeOfFlight( tof );
1923 
1924 
1925  // primary track
1926  StTrack* pTrack = gTrack->node()->track( primary );
1927  StETofPidTraits* ppidTraits = 0;
1928  if( pTrack ) {
1929  StSPtrVecTrackPidTraits& ptraits = pTrack->pidTraits();
1930  for( size_t i=0; i<ptraits.size(); i++ ) {
1931  if( ptraits[ i ]->detector() == kETofId ) {
1932  ppidTraits = dynamic_cast< StETofPidTraits* > ( ptraits[ i ] );
1933  break;
1934  }
1935  }
1936  if( ppidTraits ) {
1937  ppidTraits->setTimeOfFlight( tof );
1938  }
1939  }
1940 
1941  if( mDebug ) {
1942  LOG_INFO << "calculatePidVariables() - time-of-flight assigned to pid traits: " << tof << " ..." << endm;
1943  }
1944 
1945  // pid calculation if the track is a primary track
1946  // switch indicating to calculate PID or not
1947  bool doPID = false;
1948 
1949  double pathLength = matchCand.pathLength;
1950 
1951 
1952  if( !pTrack ) {
1953  if( mDebug ) {
1954  LOG_INFO << "calculatePidVariables() - the associated track is not a primary one. Skip PID calculation! " << endm;
1955  }
1956  }
1957  else {
1958  const StVertex* pVtx = pTrack->vertex();
1959  if( !pVtx ) {
1960  if( mDebug ) {
1961  LOG_INFO << "calculatePidVariables() - the associated track is not coming from any vertex. Skip PID calculation! " << endm;
1962  }
1963  }
1964  else {
1965  StThreeVectorD vtxPos = pVtx->position();
1966  StPhysicalHelixD theHelix = mOuterTrackGeometry ? gTrack->outerGeometry()->helix() : gTrack->geometry()->helix();
1967 
1968  pathLength -= theHelix.pathLength( vtxPos );
1969 
1970  doPID = true;
1971  }
1972  }
1973 
1974  if( !doPID ) continue;
1975 
1976  // set pathlength
1977  matchCand.pathLength = pathLength;
1978 
1979 
1980  double beta = pathLength / ( tof * nanosecond * c_light );
1981 
1982  // set beta
1983  matchCand.beta = beta;
1984 
1985  if( mDebug ) {
1986  LOG_INFO << "calculatePidVariables() - pathlength: " << pathLength << " time-of-flight: " << tof << " and beta: " << beta << " are set" << endm;
1987  }
1988 
1989  if( mDoQA ) {
1990  mHistograms.at( "matchCand_timeOfFlight" )->Fill( tof );
1991  mHistograms.at( "matchCand_timeOfFlight_pathLength" )->Fill( tof, pathLength );
1992  }
1993  mHistograms.at( "matchCand_timeOfFlight_pathLength_zoom" )->Fill( tof, pathLength );
1994 
1995  pidTraits->setPathLength( pathLength );
1996  pidTraits->setBeta( beta );
1997 
1998  if( ppidTraits ) {
1999  ppidTraits->setPathLength( pathLength );
2000  ppidTraits->setBeta( beta );
2001 
2002  nPrimaryWithPid++;
2003  }
2004 
2005  if( mDebug ) {
2006  LOG_INFO << "calculatePidVariables() - pathlength and beta assigned to pid traits ..." << endm;
2007  }
2008  }
2009  }
2010  else { // MuDst processing ...
2011  // assign pathlength, time-of-flight, beta ... to the match candidates
2012  for( auto& matchCand : finalMatchVec ) {
2013 
2014  StMuETofHit* aHit = mMuDst->etofHit( matchCand.index2ETofHit );
2015  if( !aHit ) {
2016  LOG_ERROR << "calculatePidVariables() - no hit" << endm;
2017  continue;
2018  }
2019 
2020  // global track
2021  StMuTrack* gTrack = aHit->globalTrack();
2022  if( !gTrack ) {
2023  LOG_ERROR << "calculatePidVariables() - no associated track" << endm;
2024  continue;
2025  }
2026 
2027  StMuETofPidTraits pidTraits = gTrack->etofPidTraits();
2028 
2029 
2030  //double tof = timeOfFlight( tstart, aHit->time() );
2031  double tof = timeOfFlight( tstart, matchCand.hitTime );
2032 
2033 
2034  // set time-of-flight
2035  matchCand.tof = tof;
2036 
2037 
2038  pidTraits.setTimeOfFlight( tof );
2039  gTrack->setETofPidTraits( pidTraits );
2040 
2041  // primary track
2042  StMuTrack* pTrack = aHit->primaryTrack();
2043  StMuETofPidTraits ppidTraits;
2044  if( pTrack ) {
2045  ppidTraits = pTrack->etofPidTraits();
2046 
2047  ppidTraits.setTimeOfFlight( tof );
2048  pTrack->setETofPidTraits( ppidTraits );
2049  }
2050 
2051  if( mDebug ) {
2052  LOG_INFO << "calculatePidVariables() - time-of-flight assigned to pid traits: " << tof << " ..." << endm;
2053  }
2054 
2055  // pid calculation if the track is a primary track
2056  // switch indicating to calculate PID or not
2057  bool doPID = false;
2058 
2059  double pathLength = matchCand.pathLength;
2060 
2061  if( !pTrack ) {
2062  if( mDebug ) {
2063  LOG_INFO << "calculatePidVariables() - the associated track is not a primary one. Skip PID calculation!" << endm;
2064  }
2065  }
2066  else {
2067  StMuPrimaryVertex* pVtx = mMuDst->primaryVertex( pTrack->vertexIndex() );
2068  if( !pVtx ) {
2069  if( mDebug ) {
2070  LOG_INFO << "calculatePidVariables() - the associated track is not coming from any vertex. Skip PID calculation!" << endm;
2071  }
2072  }
2073  else {
2074  StThreeVectorD vtxPos = pVtx->position();
2075  StPhysicalHelixD theHelix = mOuterTrackGeometry ? gTrack->outerHelix() : gTrack->helix();
2076 
2077  pathLength -= theHelix.pathLength( vtxPos );
2078 
2079  doPID = true;
2080  }
2081  }
2082 
2083  if( !doPID ) continue;
2084 
2085  // set pathlength
2086  matchCand.pathLength = pathLength;
2087 
2088 
2089  double beta = pathLength / ( tof * nanosecond * c_light );
2090 
2091  // set beta
2092  matchCand.beta = beta;
2093 
2094 
2095  if( matchCand.clusterSize > 999 ){
2096 
2097  mHistograms.at( "AAA_beta_mom_SD")->Fill( pTrack->momentum().mag() , 1/beta );
2098 
2099  }
2100 
2101 
2102 
2103  if( mDebug ) {
2104  LOG_INFO << "calculatePidVariables() - pathlength: " << pathLength << " time-of-flight: " << tof << " and beta: " << beta << " are set" << endm;
2105  }
2106 
2107  if( mDoQA ) {
2108  mHistograms.at( "matchCand_timeOfFlight" )->Fill( tof );
2109  mHistograms.at( "matchCand_timeOfFlight_pathLength" )->Fill( tof, pathLength );
2110  }
2111  mHistograms.at( "matchCand_timeOfFlight_pathLength_zoom" )->Fill( tof, pathLength );
2112 
2113  pidTraits.setPathLength( pathLength );
2114  pidTraits.setBeta( beta );
2115  gTrack->setETofPidTraits( pidTraits );
2116 
2117  if( pTrack ) {
2118  ppidTraits.setPathLength( pathLength );
2119  ppidTraits.setBeta( beta );
2120  pTrack->setETofPidTraits( ppidTraits );
2121 
2122  nPrimaryWithPid++;
2123  }
2124 
2125  if( mDebug ) {
2126  LOG_INFO << "calculatePidVariables() - pathlength and beta assigned to pid traits ..." << endm;
2127  }
2128  }
2129  }
2130 }
2131 
2132 
2133 //---------------------------------------------------------------------------
2134 // get the start time -- from bTOF header (default)
2135 double
2136 StETofMatchMaker::startTimeBTof()
2137 {
2138  if( mDebug ) {
2139  LOG_INFO << "startTimeBTof(): -- loading start time from bTOF header" << endm;
2140  }
2141 
2142  StBTofHeader* btofHeader = nullptr;
2143 
2144  if( mIsStEventIn ) {
2145  StBTofCollection* btofCollection = ( StBTofCollection* ) mEvent->btofCollection();
2146 
2147  if ( btofCollection ) {
2148  btofHeader = btofCollection->tofHeader();
2149  }
2150  else {
2151  LOG_DEBUG << "startTimeBTof(): -- no StBTofCollection found by getTstart" << endm;
2152  return -9999.;
2153  }
2154  }
2155  else if( mIsMuDstIn ) {
2156  btofHeader = mMuDst->btofHeader();
2157  }
2158 
2159  if( !btofHeader ) {
2160  LOG_DEBUG << "startTimeBTof(): -- no bTOF header --> no start time avaiable" << endm;
2161  return -9999.;
2162  }
2163 
2164  double tstart = btofHeader->tStart();
2165 
2166  if( !isfinite( tstart ) ) {
2167  LOG_DEBUG << "startTimeBTof(): -- from bTOF header is NaN" << endm;
2168  return -9999.;
2169  }
2170 
2171  if( tstart != -9999. ) {
2172  tstart = fmod( tstart, eTofConst::bTofClockCycle );
2173  if( tstart < 0 ) tstart += eTofConst::bTofClockCycle;
2174  }
2175 
2176  if( mDebug ) {
2177  LOG_DEBUG << "startTimeBTof(): -- start time: " << tstart << endm;
2178  }
2179 
2180  return tstart;
2181 }
2182 
2183 //---------------------------------------------------------------------------
2184 // get the start time -- only from eTOF matches
2185 double
2186 StETofMatchMaker::startTimeETof( const eTofHitVec& finalMatchVec, unsigned int& nCand_etofT0 )
2187 {
2188  if( mDebug ) {
2189  LOG_INFO << "startTimeETof(): -- calculating start time from eTOF matches" << endm;
2190  }
2191  std::vector< double > t0_cand;
2192  double t0_sum = 0.;
2193 
2194  nCand_etofT0 = 0;
2195 
2196  for( auto& match : finalMatchVec ) {
2197  if( !match.isPrimary ) continue;
2198 
2199  if( sqrt( pow( match.deltaX, 2 ) + pow( match.deltaY, 2 ) ) > etofProjection::deltaRcut ) continue;
2200 
2201  double pathLength = match.pathLength;
2202 
2203  double momentum;
2204  int charge;
2205 
2206  double nsigmaPi = -999.;
2207  double nsigmaK = -999.;
2208  double nsigmaP = -999.;
2209 
2210  if( mIsStEventIn ) { // StEvent processing ...
2211  StSPtrVecTrackNode& nodes = mEvent->trackNodes();
2212 
2213  // global track
2214  StGlobalTrack* gTrack = dynamic_cast< StGlobalTrack* > ( nodes[ match.trackId ]->track( global ) );
2215  if( !gTrack ) {
2216  continue;
2217  }
2218 
2219  // primary track
2220  StPrimaryTrack* pTrack = dynamic_cast< StPrimaryTrack* > ( gTrack->node()->track( primary ) );
2221  if( !pTrack ) {
2222  continue;
2223  }
2224 
2225 
2226  momentum = pTrack->geometry()->momentum().mag();
2227  charge = pTrack->geometry()->charge();
2228 
2229  static StTpcDedxPidAlgorithm PidAlgorithm;
2230  static StPionPlus* Pion = StPionPlus::instance();
2231  static StKaonPlus* Kaon = StKaonPlus::instance();
2232  static StProton* Proton = StProton::instance();
2233  const StParticleDefinition* pd = pTrack->pidTraits( PidAlgorithm );
2234 
2235  if( pd && PidAlgorithm.traits() ) {
2236  nsigmaPi = PidAlgorithm.numberOfSigma( Pion );
2237  nsigmaK = PidAlgorithm.numberOfSigma( Kaon );
2238  nsigmaP = PidAlgorithm.numberOfSigma( Proton );
2239  }
2240 
2241  StPhysicalHelixD theHelix = mOuterTrackGeometry ? gTrack->outerGeometry()->helix() : gTrack->geometry()->helix();
2242  pathLength -= theHelix.pathLength( pTrack->vertex()->position() );
2243  }
2244  else { // MuDst processing ...
2245  StMuETofHit* aHit = mMuDst->etofHit( match.index2ETofHit );
2246  if( !aHit ) continue;
2247 
2248  // global track
2249  StMuTrack* gTrack = aHit->globalTrack();
2250  if( !gTrack ) continue;
2251 
2252  // primary track
2253  StMuTrack* pTrack = aHit->primaryTrack();
2254  if( !pTrack ) continue;
2255 
2256 
2257  momentum = pTrack->momentum().mag();
2258  charge = pTrack->charge();
2259 
2260  nsigmaPi = pTrack->nSigmaPion();
2261  nsigmaK = pTrack->nSigmaKaon();
2262  nsigmaP = pTrack->nSigmaProton();
2263 
2264  StPhysicalHelixD theHelix = mOuterTrackGeometry ? gTrack->outerHelix() : gTrack->helix();
2265  pathLength -= theHelix.pathLength( mMuDst->primaryVertex( pTrack->vertexIndex() )->position() );
2266  }
2267 
2268  if( momentum < 0.2 ) continue;
2269 
2270  double tofExpect;
2271 
2272  // identyfy paricles via dE/dx
2273  if( momentum < 0.6 && fabs( nsigmaPi ) < 2. ) {
2274  tofExpect = expectedTimeOfFlight( pathLength, momentum, pion_minus_mass_c2 );
2275  if( mDoQA ) {
2276  LOG_INFO << "startTimeETof(): -- using a pion as start time candidate" << endm;
2277  }
2278  }
2279  else if( momentum < 0.5 && fabs( nsigmaK ) < 1. ) {
2280  tofExpect = expectedTimeOfFlight( pathLength, momentum, kaon_minus_mass_c2 );
2281  if( mDoQA ) {
2282  LOG_INFO << "startTimeETof(): -- using a kaon as start time candidate" << endm;
2283  }
2284  }
2285  else if( momentum < 0.9 && charge == 1 && fabs( nsigmaP ) < 2. ) {
2286  tofExpect = expectedTimeOfFlight( pathLength, momentum, proton_mass_c2 );
2287  if( mDoQA ) {
2288  LOG_INFO << "startTimeETof(): -- using a proton as start time candidate" << endm;
2289  }
2290  }
2291  else{
2292  continue;
2293  }
2294 
2295  t0_cand.push_back( match.hitTime - tofExpect );
2296  t0_sum += t0_cand.back();
2297 
2298  if( mDebug ) {
2299  LOG_INFO << match.hitTime << " " << tofExpect << " " << t0_cand.back() << endm;
2300  }
2301  }
2302 
2303  nCand_etofT0 = t0_cand.size();
2304 
2305  if( nCand_etofT0 == 0 ) {
2306  return -9999.;
2307  }
2308  else if( nCand_etofT0 > 1 ) { // remove hits too far from others
2309  for( unsigned int i=0; i < nCand_etofT0; i++ ) {
2310 
2311  double t0_diff = t0_cand.at( i ) - ( t0_sum - t0_cand.at( i ) ) / ( nCand_etofT0 - 1 );
2312 
2313  if( fabs( t0_diff ) > 5.0 ) {
2314  t0_sum -= t0_cand.at( i );
2315  nCand_etofT0--;
2316  }
2317  }
2318  }
2319 
2320  if( nCand_etofT0 < 2 ) {
2321  return -9999.;
2322  }
2323 
2324 
2325  double tStart = fmod( t0_sum / nCand_etofT0, eTofConst::bTofClockCycle );
2326  if( tStart < 0 ) tStart += eTofConst::bTofClockCycle;
2327 
2328  return tStart;
2329 }
2330 
2331 
2332 //---------------------------------------------------------------------------
2333 // distance on a circle
2334 double
2335 StETofMatchMaker::moduloDist( const double& dist, const double& mod ) {
2336  return dist - mod * round( dist / mod );
2337 }
2338 
2339 //---------------------------------------------------------------------------
2340 // calculate hybrid start time with eTOF and bTOF information
2341 double
2342 StETofMatchMaker::startTime( const eTofHitVec& finalMatchVec ) {
2343  // for simulation tstart == 0
2344  if( mIsSim ) {
2345  return 0.;
2346  }
2347 
2348  double tstartBTof = startTimeBTof();
2349 
2350  if( mUseOnlyBTofHeaderStartTime ) {
2351  return tstartBTof;
2352  }
2353 
2354 
2355  unsigned int nCand_etofT0 = 0;
2356  double tstartETof = startTimeETof( finalMatchVec, nCand_etofT0 );
2357 
2358  double t0Diff = moduloDist( tstartETof - tstartBTof, eTofConst::bTofClockCycle );
2359 
2360  //LOG_INFO << "startTime(): -- start time comparison: bTOF " << tstartBTof << " eTOF " << tstartETof;
2361  //LOG_INFO << " nCand_etofT0: " << nCand_etofT0 << " difference: " << t0Diff << " mT0corr: " << mT0corr << endm;
2362 
2363  if( mDoQA && tstartBTof != -9999. && tstartETof != -9999. ) {
2364  mHistograms.at( "startTimeDiff" )->Fill( t0Diff );
2365  mHistograms.at( "startTimeDiff_nCand" )->Fill( t0Diff, nCand_etofT0 );
2366  }
2367 
2368  //---------------------------------
2369  // calculate T0 corr for hybrid start time
2370  if( tstartBTof != -9999. && tstartETof != -9999. && nCand_etofT0 > etofHybridT0::minCand ) {
2371  if( mT0corr != 0 && fabs( moduloDist( t0Diff - mT0corr, eTofConst::bTofClockCycle ) ) > etofHybridT0::diffToClear ) {
2372  mT0switch++;
2373  if( mDebug ) {
2374  LOG_INFO << "mT0switch: " << mT0switch << endm;
2375  }
2376 
2377  if( mT0switch > etofHybridT0::nSwitch ) {
2378  mT0corrVec.clear();
2379  mT0corrVec.push_back( t0Diff );
2380  mT0switch = 0;
2381  mT0corr = 0.;
2382  if( mDebug ) {
2383  LOG_INFO << "startTime(): -- cleared T0 correction vector since the start time (eTof <-> bTOF) seems to have changed" << endm;
2384  }
2385  }
2386  }
2387  else {
2388  mT0corrVec.push_back( t0Diff );
2389  mT0switch = 0;
2390  }
2391 
2392  if( ( mT0corr != 0. || mT0corrVec.size() >= etofHybridT0::nMin ) && mT0corrVec.size() % etofHybridT0::nUpdate == 0 ) {
2393  std::sort( mT0corrVec.begin(), mT0corrVec.end() );
2394 
2395  if( mDebug ) {
2396  for( const auto& v : mT0corrVec ) {
2397  LOG_DEBUG << "startTime(): -- T0corrVec: " << v << endm;
2398  }
2399  }
2400 
2401 
2402  // preselection if mT0corr is not yet filled with information:
2403  // reject outliers far away from the median
2404  if( mT0corr == 0. ) {
2405  mT0corr = mT0corrVec.at( mT0corrVec.size() / 2 );
2406  for( auto it = mT0corrVec.begin(); it != mT0corrVec.end(); it++ ) {
2407  if( fabs( *it - mT0corr ) > 0.5 ) {
2408  mT0corrVec.erase( it-- );
2409  }
2410  }
2411  mT0corr = mT0corrVec.at( mT0corrVec.size() / 2 );
2412 
2413  //reject outliers
2414  for( auto it = mT0corrVec.begin(); it != mT0corrVec.end(); it++ ) {
2415  if( fabs( *it - mT0corr ) > 0.2 ) {
2416  mT0corrVec.erase( it-- );
2417  }
2418  }
2419  mT0corr = mT0corrVec.at( mT0corrVec.size() / 2 );
2420  }
2421 
2422 
2423  //reject outliers
2424  for( auto it = mT0corrVec.begin(); it != mT0corrVec.end(); it++ ) {
2425  if( fabs( *it - mT0corr ) > 0.5 ) {
2426  mT0corrVec.erase( it-- );
2427  }
2428  }
2429 
2430  if( mDebug ) {
2431  for( const auto& v : mT0corrVec ) {
2432  LOG_DEBUG << "startTime(): -- cleaned T0corrVec: " << v << endm;
2433  }
2434  }
2435 
2436 
2437  // reject some fraction of the early and late matches
2438  int low = floor( mT0corrVec.size() * etofHybridT0::lowCut );
2439  int high = floor( mT0corrVec.size() * etofHybridT0::highCut );
2440 
2441  double temp = 0.;
2442  int nAccept = 0;
2443  for( int i = low; i < high; i++ ) {
2444  if( mDebug ) {
2445  LOG_DEBUG << "startTime(): -- T0corrVec: " << mT0corrVec.at( i ) << endm;
2446  }
2447 
2448  temp += mT0corrVec.at( i );
2449  nAccept++;
2450  }
2451 
2452  if( nAccept >= 5 ) {
2453  mT0corr = ( temp / nAccept ) - etofHybridT0::biasOffset;
2454  mNupdatesT0++;
2455 
2456  if( mDoQA ) {
2457  LOG_INFO << "startTime(): -- hybrid T0 correction between eTOF & bTOF (update " << mNupdatesT0 << "): " << mT0corr << " (" << nAccept << ")" << endm;
2458 
2459  mHistograms.at( "startTime_T0corr" )->SetBinContent( mNupdatesT0 , mT0corr );
2460  }
2461  }
2462  }
2463  }
2464  //---------------------------------
2465 
2466  double tstart;
2467 
2468  if( mT0corr != 0. && tstartBTof != -9999. ) {
2469  tstart = fmod( moduloDist( tstartBTof + mT0corr, eTofConst::bTofClockCycle ), eTofConst::bTofClockCycle );
2470  if( tstart < 0. ) tstart += eTofConst::bTofClockCycle;
2471 
2472  if( mDoQA ) {
2473  LOG_INFO << "startTime(): -- returning hybrid start time (BTof/VPD T0: " << tstartBTof << " T0 corr: " << mT0corr << "): " << tstart << endm;
2474  }
2475  }
2476  else if( tstartETof != -9999. && nCand_etofT0 > etofHybridT0::minCand ) {
2477  tstart = tstartETof;
2478  if( mDoQA ) {
2479  LOG_INFO << "startTime(): -- returning eTOF-only start time: " << tstart << endm;
2480  }
2481  }
2482  else {
2483  tstart = tstartBTof;
2484 
2485  if( mDoQA ) {
2486  LOG_INFO << "startTime(): -- returning bTOF/VPD start time: " << tstart << endm;
2487  }
2488  }
2489 
2490  return tstart;
2491 }
2492 
2493 
2494 //---------------------------------------------------------------------------
2495 // calculate measured time of flight
2496 double
2497 StETofMatchMaker::timeOfFlight( const double& startTime, const double& stopTime )
2498 {
2499  return moduloDist( stopTime - startTime, eTofConst::bTofClockCycle );
2500 }
2501 
2502 
2503 //---------------------------------------------------------------------------
2504 // calculate expected time of flight of particle hypothesis
2505 double
2506 StETofMatchMaker::expectedTimeOfFlight( const double& pathLength, const double& momentum, const double& mass )
2507 {
2508  double inverseBeta = sqrt( 1 + mass * mass / ( momentum * momentum ) );
2509 
2510  return pathLength * centimeter * ( 1. / c_light ) * inverseBeta / nanosecond;
2511 }
2512 
2513 
2514 //.........................................................................
2515 // H. fill QA histograms
2516 //
2517 //---------------------------------------------------------------------------
2518 void
2519 StETofMatchMaker::fillQaHistograms( eTofHitVec& finalMatchVec )
2520 {
2521 
2522  vector< int > nPidMatches( 36 );
2523 
2524  for( auto& matchCand : finalMatchVec ) {
2525 
2526  int charge;
2527  float mom;
2528 
2529  // int sector = 0; //
2530  // int plane = 0; //
2531  // int counter = 0; //
2532 
2533  float dEdx = -999.;
2534  float nSigmaPion = -999;
2535 
2536  if( mIsStEventIn ) {
2537 
2538  StSPtrVecTrackNode& nodes = mEvent->trackNodes();
2539 
2540  StGlobalTrack* gTrack = dynamic_cast< StGlobalTrack* > ( nodes[ matchCand.trackId ]->track( global ) );
2541  if( !gTrack ) continue;
2542 
2543  StPrimaryTrack* pTrack = dynamic_cast< StPrimaryTrack* > ( gTrack->node()->track( primary ) );
2544  if( !pTrack ) continue;
2545 
2546  charge = pTrack->geometry()->charge();
2547  mom = pTrack->geometry()->momentum().mag();
2548 
2549  static StTpcDedxPidAlgorithm PidAlgorithm;
2550  static StPionPlus* Pion = StPionPlus::instance();
2551  const StParticleDefinition* pd = pTrack->pidTraits( PidAlgorithm );
2552 
2553  if( pd && PidAlgorithm.traits() ) {
2554  dEdx = PidAlgorithm.traits()->mean() * 1.e6;
2555  nSigmaPion = PidAlgorithm.numberOfSigma( Pion );
2556  }
2557  }
2558  else {
2559 
2560  StMuETofHit* aHit = mMuDst->etofHit( matchCand.index2ETofHit );
2561  if( !aHit ) continue;
2562 
2563  StMuTrack* pTrack = aHit->primaryTrack();
2564  if( !pTrack ) continue;
2565 
2566  //sector = aHit->sector();
2567  //plane = aHit->zPlane();
2568  //counter = aHit->counter();
2569 
2570  charge = pTrack->charge();
2571  mom = pTrack->momentum().mag();
2572 
2573  dEdx = pTrack->dEdx() * 1.e6;
2574  nSigmaPion = pTrack->nSigmaPion();
2575  }
2576 
2577  int sign = ( charge < 0 ) ? -1 : ( charge > 0 );
2578  float beta = matchCand.beta;
2579 
2580  // skip events with no valid start time: beta of matchCand structure != -999.
2581  if( fabs( beta + 999. ) < 0.001 ) {
2582  if( mDoQA ) {
2583  LOG_WARN << "beta not set --> no start time available???" << endm;
2584  }
2585  continue;
2586  }
2587 
2588  // expected time of flight for mass hypothesis pion
2589  float tofpi = expectedTimeOfFlight( matchCand.pathLength, mom, pion_plus_mass_c2 );
2590 
2591  mHistograms.at( "matchCand_beta_signmom" )->Fill( sign * mom, 1. / beta );
2592  if( mom > 0.6 && mom < 1.5 ) {
2593  mHistograms.at( "matchCand_t0corr_1d" )->Fill( matchCand.tof - tofpi );
2594  }
2595 
2596  if( mDoQA ) {
2597  float tof = matchCand.tof;
2598  float pathlength = matchCand.pathLength;
2599  float m2 = mom * mom * ( -1 + 1 / ( beta * beta ) );
2600 
2601  //LOG_INFO << "momentum: " << mom << " ... beta: " << beta << " ... m^2: " << m2 << " ... dEdx: " << dEdx << endm;
2602  //LOG_DEBUG << "pathlength: " << pathlength << " ... time-of-flight: " << matchCand.tof << " ... mom: " << mom << " ... beta: " << beta << " ... charge: " << charge << " ... m^2: " << m2 << " ... dEdx: " << dEdx << endm;
2603 
2604  if( fabs( tof+999. ) < 1.e-5 || fabs( pathlength+999. ) < 1.e-5 ) {
2605  LOG_INFO << "tof==0 or pathlength==0 ..." << endl;
2606  continue;
2607  }
2608 
2609  mHistograms.at( "matchCand_beta_mom" )->Fill( mom, 1. / beta );
2610  mHistograms.at( "matchCand_m2_mom" )->Fill( mom, m2 );
2611  mHistograms.at( "matchCand_m2_signmom" )->Fill( sign * mom, m2 );
2612 
2613 
2614 
2615  // plots per counter
2616  std::string histName_beta_mom = "matchCand_beta_mom_s" + std::to_string( matchCand.sector ) + "m" + std::to_string( matchCand.plane ) + "c" + std::to_string( matchCand.counter );
2617  mHistograms.at( histName_beta_mom )->Fill( mom, 1. / beta );
2618 
2619  std::string histName_t0corr_mom = "matchCand_t0corr_mom_s" + std::to_string( matchCand.sector ) + "m" + std::to_string( matchCand.plane ) + "c" + std::to_string( matchCand.counter );
2620  mHistograms.at( histName_t0corr_mom )->Fill( mom, tof - tofpi );
2621 
2622  std::string histName_t0corr_mom_zoom = "matchCand_t0corr_mom_zoom_s" + std::to_string( matchCand.sector ) + "m" + std::to_string( matchCand.plane ) + "c" + std::to_string( matchCand.counter );
2623  mHistograms.at( histName_t0corr_mom_zoom )->Fill( mom, tof - tofpi );
2624 
2625  if( nSigmaPion < 1.5 ) {
2626  std::string histName_t0corr_strip = "matchCand_t0corr_strip_s" + std::to_string( matchCand.sector ) + "m" + std::to_string( matchCand.plane ) + "c" + std::to_string( matchCand.counter );
2627  mHistograms.at( histName_t0corr_strip )->Fill( matchCand.localX, tof - tofpi );
2628  }
2629 
2630  if( nSigmaPion < 2. ) {
2631  if( matchCand.clusterSize >= 100 ) {
2632  // hits with clock jump based on local Y position
2633  std::string histName_t0corr_jump = "matchCand_t0corr_jump_s" + std::to_string( matchCand.sector ) + "m" + std::to_string( matchCand.plane ) + "c" + std::to_string( matchCand.counter );
2634  mHistograms.at( histName_t0corr_jump )->Fill( matchCand.localX, tof - tofpi ); //cutdownYS
2635 
2636  int get4Index = matchCand.sector * 1000 + matchCand.plane * 100 + matchCand.counter * 10 + ( matchCand.localX + 16 ) / 4 + 1;
2637  mClockJumpCand[ get4Index ]++;
2638  }
2639  }
2640 
2641 
2642  if( sqrt( pow( matchCand.deltaX, 2 ) + pow( matchCand.deltaY, 2 ) ) < etofProjection::deltaRcut ) {
2643 
2644  mHistograms.at( "matchCand_beta_mom_matchDistCut" )->Fill( mom, 1. / beta );
2645  mHistograms.at( "matchCand_m2_mom_matchDistCut" )->Fill( mom, m2 );
2646 
2647  std::string histName_t0corr_mom_zoom_cut = "matchCand_t0corr_mom_zoom_cut_s" + std::to_string( matchCand.sector ) + "m" + std::to_string( matchCand.plane ) + "c" + std::to_string( matchCand.counter );
2648 
2649  mHistograms.at( histName_t0corr_mom_zoom_cut )->Fill( mom, tof - tofpi );
2650 
2651  nPidMatches.at( ( matchCand.sector - 13 ) * 3 + matchCand.plane - 1 ) += 1;
2652  }
2653 
2654  if( fabs( mom - 1 ) < 0.1 && dEdx > 0 ) mHistograms.at( "matchCand_dEdx_beta_mom1gev" )->Fill( 1. / beta, dEdx );
2655  if( fabs( mom - 2 ) < 0.1 && dEdx > 0 ) mHistograms.at( "matchCand_dEdx_beta_mom2gev" )->Fill( 1. / beta, dEdx );
2656  }
2657  }
2658 
2659  if( mDoQA ) {
2660  for( size_t i=0; i<36; i++ ) {
2661  mHistograms.at( "matchCandMultPerSector_matchDistCut" )->Fill( i, nPidMatches.at( i ) );
2662  }
2663  }
2664 }
2665 
2666 void
2667 StETofMatchMaker::fillSlewHistograms( eTofHitVec& finalMatchVec )
2668 {
2669  if( !mDoQA || !mIsMuDstIn ) return;
2670 
2671  map< int, float > hitIndexToDeltaTmap;
2672 
2673  for( auto& matchCand : finalMatchVec ) {
2674  StMuETofHit* aHit = mMuDst->etofHit( matchCand.index2ETofHit );
2675  if( !aHit ) continue;
2676 
2677  StMuTrack* pTrack = aHit->primaryTrack();
2678  if( !pTrack ) continue;
2679 
2680  float mom = pTrack->momentum().mag();
2681  float nSigmaPion = pTrack->nSigmaPion();
2682 
2683  if( fabs( nSigmaPion > 1.5 ) ) continue;
2684 
2685 
2686  // skip events with no valid start time: beta of matchCand structure != -999.
2687  if( fabs( matchCand.beta + 999. ) < 0.001 ) {
2688  LOG_INFO << "beta not set --> no start time available???" << endm;
2689  continue;
2690  }
2691 
2692  float tof = matchCand.tof;
2693  float pathlength = matchCand.pathLength;
2694 
2695  if( fabs( tof+999. ) < 1.e-5 || fabs( pathlength+999. ) < 1.e-5 ) {
2696  LOG_INFO << "tof==0 or pathlength==0 ..." << endl;
2697  continue;
2698  }
2699 
2700  LOG_DEBUG << "for slewing: momentum: " << mom << " ... nsigmaPi: " << nSigmaPion << " tof:" << tof << endm;
2701 
2702  // expected time of flight for mass hypothesis pion
2703  float tofpi = expectedTimeOfFlight( pathlength , mom, pion_plus_mass_c2 );
2704  hitIndexToDeltaTmap[ matchCand.index2ETofHit ] = tof - tofpi;
2705  }
2706 
2707  //find digis that belong to the matched hits
2708  size_t nDigis = mMuDst->numberOfETofDigi();
2709  for( size_t i=0; i<nDigis; i++ ) {
2710  StMuETofDigi* aDigi = mMuDst->etofDigi( i );
2711  if( !aDigi ) continue;
2712 
2713  if( hitIndexToDeltaTmap.count( aDigi->associatedHitId() ) ) {
2714  std::string histName_slewing = "matchCand_slewing_digi_s" + std::to_string( aDigi->sector() ) + "m" + std::to_string( aDigi->zPlane() ) + "c" + std::to_string( aDigi->counter() );
2715  mHistograms.at( histName_slewing )->Fill( aDigi->calibTot(), hitIndexToDeltaTmap.at( aDigi->associatedHitId() ) );
2716  }
2717  }
2718 }
2719 
2720 
2721 //---------------------------------------------------------------------------
2722 // returns the proper track geometry, based on a global user setting
2723 StTrackGeometry*
2724 StETofMatchMaker::trackGeometry( StTrack* track ) const
2725 {
2726  // returns apropriate StTrackGeometry (standard or outerGeometry)
2727  if ( !track ) return nullptr;
2728 
2729  StTrackGeometry* thisTrackGeometry;
2730 
2731  if ( mOuterTrackGeometry ) {
2732  thisTrackGeometry = track->outerGeometry();
2733  }
2734  else {
2735  thisTrackGeometry = track->geometry();
2736  }
2737 
2738  return thisTrackGeometry;
2739 }
2740 
2741 
2742 //---------------------------------------------------------------------------
2743 // setHistFileName uses the string argument from the chain being run to set
2744 // the name of the output histogram file.
2745 void
2746 StETofMatchMaker::setHistFileName()
2747 {
2748  string extension = ".etofMatch.root";
2749 
2750  if( GetChainOpt()->GetFileOut() != nullptr ) {
2751  TString outFile = GetChainOpt()->GetFileOut();
2752 
2753  mHistFileName = ( string ) outFile;
2754 
2755  // get rid of .root
2756  size_t lastindex = mHistFileName.find_last_of( "." );
2757  mHistFileName = mHistFileName.substr( 0, lastindex );
2758 
2759  // get rid of .MuDst or .event if necessary
2760  lastindex = mHistFileName.find_last_of( "." );
2761  mHistFileName = mHistFileName.substr( 0, lastindex );
2762 
2763  // get rid of directories
2764  lastindex = mHistFileName.find_last_of( "/" );
2765  mHistFileName = mHistFileName.substr( lastindex + 1, mHistFileName.length() );
2766 
2767  mHistFileName = mHistFileName + extension;
2768  } else {
2769  LOG_ERROR << "Cannot set the output filename for histograms" << endm;
2770  }
2771 }
2772 
2773 
2774 //---------------------------------------------------------------------------
2775 void
2776 StETofMatchMaker::bookHistograms()
2777 {
2778  LOG_INFO << "bookHistograms" << endm;
2779 
2780  mHistograms[ "eTofHits_globalXY" ] = new TH2F( "A_eTofHits_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2781  mHistograms[ "intersectionMult_etofMult" ] = new TH2F( "B_intersectionMult_etofMult", "multiplicity correlation;# eTOF hits;#track intersections;# events", 300, 0., 300., 200, 0., 200. );
2782  mHistograms[ "matchCand_globalXY" ] = new TH2F( "C_matchCand_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2783  mHistograms[ "matchCand_beta_signmom" ] = new TH2F( "G_matchCand_beta_signmom" , "match candidate 1/beta vs. momentum;q/|q| * p (GeV/c);1/#beta", 400, -4., 6., 1000, 0.8, 2. );
2784  mHistograms[ "matchCand_timeOfFlight_pathLength_zoom" ] = new TH2F( "G_matchCand_timeOfFlight_pathLength", "match candidate pathlength vs. time of flight;ToF (ns);pathlength (cm)", 800, -25., 75., 800, 200., 600. );
2785  mHistograms[ "primaryIntersect_validMatch" ] = new TH2F( "G_primary_Intersection_validMatch", "primary tracks at eTOF;# tracks with intersection;# tracks with valid PID", 200, 0., 200., 100, 0., 100. );
2786  mHistograms[ "matchCand_t0corr_1d" ] = new TH1F( "H_matchCand_t0corr_1d", "measured tof - tof_{#pi};#Delta time (ns);counts", 3000, -15., 15. );
2787 
2788 
2789  AddHist( mHistograms.at( "eTofHits_globalXY" ) );
2790  AddHist( mHistograms.at( "intersectionMult_etofMult" ) );
2791  AddHist( mHistograms.at( "matchCand_globalXY" ) );
2792  AddHist( mHistograms.at( "matchCand_beta_signmom" ) );
2793  AddHist( mHistograms.at( "matchCand_timeOfFlight_pathLength_zoom" ) );
2794  AddHist( mHistograms.at( "primaryIntersect_validMatch" ) );
2795  AddHist( mHistograms.at( "matchCand_t0corr_1d" ) );
2796 
2797  if( mDoQA ) {
2798  // histograms to test sector & module numbering
2799  mHistograms[ "eTofSectors" ] = new TH2F( "QA_eTofSectors", "center of modules in global XY;x (cm);y (cm)", 100, -300, 300, 100, -300, 300 );
2800  mHistograms[ "eTofModules" ] = new TH2F( "QA_eTofModules", "center of modules in global XY;x (cm);y (cm)", 100, -300, 300, 100, -300, 300 );
2801 
2802  // histograms testing the magnetic field
2803  mHistograms2d[ "bfield_z" ] = new TH2F( "QA_bField_z", "magnetic field (z);z (cm);y (cm)", 350, -350., 0., 300, 0., 300. );
2804 
2805 
2806  // event-level QA
2807  mHistograms[ "eventCounter" ] = new TH1F( "QA_eventCounter","eventCounter;;# events", 10, 0.5, 10.5 );
2808 
2809 
2810  // ----------
2811  // step - A -
2812  // ----------
2813  mHistograms[ "eTofHits_phi_eta" ] = new TH2F( "A_eTofHits_phi_eta", "eta vs. phi; #phi; #eta", 400, 0., 2 * M_PI, 400, -1.7, -0.9 );
2814 
2815  mHistograms[ "eTofHits_globalYZ" ] = new TH2F( "A_eTofHits_globalYZ", "global YZ for sector 18 & 24;y (cm);z (cm)", 400, -300., 300., 100, -310., -270. );
2816 
2817  for( int sector = eTofConst::sectorStart; sector <= eTofConst::sectorStop; sector++ ) {
2818  for( int plane = eTofConst::zPlaneStart; plane <= eTofConst::zPlaneStop; plane++ ) {
2819  for( int counter = eTofConst::counterStart; counter <= eTofConst::counterStop; counter++ ) {
2820 
2821  //single sided matching qa
2822  std::string histName_t0corr_mom_zoom = "matched_t0corr_mom_zoom_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2823 
2824  mHistograms2d[ histName_t0corr_mom_zoom ] = new TH2F( Form( "T_matched_t0corr_mom_zoom_s%dm%dc%d", sector, plane, counter ), Form( "measured tof - tof_{#pi} vs. momentum in sector %d module %d counter %d;mom (GeV/c);#Delta time (ns)", sector, plane, counter ), 200, 0., 3., 500, -5., 5. );
2825 
2826 
2827  std::string histName_t0corr_mom_zoom_SD = "matched_t0corr_mom_zoom_SD_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2828 
2829  mHistograms2d[ histName_t0corr_mom_zoom_SD ] = new TH2F( Form( "T_matched_t0corr_mom_zoom_SD_s%dm%dc%d", sector, plane, counter ), Form( "measured tof - tof_{#pi} vs. momentum in sector %d module %d counter %d;mom (GeV/c);#Delta time (ns)", sector, plane, counter ), 200, 0., 3., 500, -5., 5. );
2830 
2831 
2832 
2833 
2834 
2835 
2836  std::string histName_hit_localXY = "eTofHits_localXY_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2837  std::string histName_hit_globalXY = "eTofHits_globalXY_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2838  std::string histName_hit_eta_phi = "eTofHits_phi_eta_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2839 
2840  // local XY per counter
2841  mHistograms[ histName_hit_localXY ] = new TH2F( Form( "A_eTofHits_localXY_s%dm%dc%d", sector, plane, counter ), Form( "local XY sector %d module %d counter %d;x (cm);y (cm)", sector, plane, counter ), 40, -20., 20., 50, -25., 25. );
2842 
2843  // global XY per counter
2844  mHistograms[ histName_hit_globalXY ] = new TH2F( Form( "A_eTofHits_globalXY_s%dm%dc%d", sector, plane, counter ), Form( "global XY sector %d module %d counter %d;x (cm);y (cm)", sector, plane, counter ), 200, -300., 300., 200, -300., 300. );
2845 
2846  // eta vs. phi per counter
2847  mHistograms[ histName_hit_eta_phi ] = new TH2F( Form( "A_eTofHits_phi_eta_s%dm%dc%d", sector, plane, counter ), Form( "eta vs. phi sector %d module %d counter %d; #phi; #eta", sector, plane, counter ), 200, 0., 2 * M_PI, 200, -1.7, -0.9 );
2848 
2849 
2850  std::string histName_t0corr_jump = "matchCand_t0corr_jump_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2851  mHistograms[ histName_t0corr_jump ] = new TH2F( Form( "H_matchCand_t0corr_jump_s%dm%dc%d", sector, plane, counter ), Form( "measured tof - tof_{#pi} vs. momentum in sector %d module %d counter %d;localX (cm) grouped by Get4;#Delta time (ns)", sector, plane, counter ), 8, -16., 16., 80, -20., 20. );
2852 
2853 
2854 
2855  }
2856  }
2857  }
2858 
2859  mHistograms[ "detectorHitMult" ] = new TH1F( "A_detectorHitMult", "detectorHitMult;multiplicity;# events", 200, 0., 200. );
2860 
2861 
2862  // --------------
2863  // track-level QA
2864  // --------------
2865  mHistograms[ "track_phi_eta" ] = new TH2F( "QA_track_phi_eta", "eta vs. phi; #phi; #eta", 400, 0., 2 * M_PI, 800, -1.9, 1.9 );
2866  mHistograms[ "track_phi_pt" ] = new TH2F( "QA_track_phi_pt", "pt vs. phi; #phi; p_{T}", 400, 0., 2 * M_PI, 400, 0., 5. );
2867 
2868  mHistograms[ "nHits" ] = new TH1F( "QA_nHits", "nHitsTpc;nHitsFit;# tracks", 75, 0., 75. );
2869  mHistograms[ "nHits_etofregion" ] = new TH1F( "QA_nHits_etofregion", "nHitsTpc in etof region;nHitsFit;# tracks", 75, 0., 75. );
2870 
2871  mHistograms[ "track_pt_nHits" ] = new TH2F( "QA_track_pt_nHits", "track nHitsTpc vs. p_{T};p_{T} (GeV/c);nHitsFit", 400, 0., 2., 75, 0., 75. );
2872 
2873  mHistograms[ "trackProj_globalXY" ] = new TH2F( "QA_trackProj_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2874  mHistograms[ "trackProj_phi_eta" ] = new TH2F( "QA_trackProj_phi_eta", "eta vs. phi; #phi; #eta", 400, 0., 2 * M_PI, 400, -1.7, -1.0 );
2875 
2876 
2877  // ----------
2878  // step - B -
2879  // ----------
2880 
2881  mHistograms[ "intersection_globalXY" ] = new TH2F( "B_intersection_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2882  mHistograms[ "intersection_phi_eta" ] = new TH2F( "B_intersection_phi_eta", "eta vs. phi; #phi; #eta", 400, 0., 2 * M_PI, 400, -1.7, -1.0 );
2883 
2884  mHistograms[ "intersectionMult" ] = new TH1F( "B_intersectionMult", "intersectionMult;multiplicity;# events", 200, 0., 200. );
2885  mHistograms[ "intersectionMult_primary" ] = new TH1F( "B_intersectionMult_primary", "intersectionMult for primary tracks;multiplicity;# events", 200, 0., 200. );
2886 
2887  mHistograms[ "intersection_perTrack" ] = new TH1F( "B_intersection_perTrack", "intersections per track;# intersections;# tracks", 50, 0., 50. );
2888 
2889 
2890  // track-level QA plots for tracks that have an intersection with eTof volume(s)
2891  mHistograms[ "intersection_primaryTrack_globalXY" ] = new TH2F( "B_intersection_primaryTrack_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2892  mHistograms[ "intersection_primaryTrackMom0_globalXY" ] = new TH2F( "B_intersection_primaryTrackMom0_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2893  mHistograms[ "intersection_primaryTrackMom1_globalXY" ] = new TH2F( "B_intersection_primaryTrackMom1_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2894  mHistograms[ "intersection_primaryTrackMom2_globalXY" ] = new TH2F( "B_intersection_primaryTrackMom2_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2895 
2896  mHistograms[ "intersection_primaryTrackpos_globalXY" ] = new TH2F( "B_intersection_primaryTrackpos_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2897  mHistograms[ "intersection_primaryTrackMom0pos_globalXY" ] = new TH2F( "B_intersection_primaryTrackMom0pos_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2898  mHistograms[ "intersection_primaryTrackMom1pos_globalXY" ] = new TH2F( "B_intersection_primaryTrackMom1pos_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2899  mHistograms[ "intersection_primaryTrackMom2pos_globalXY" ] = new TH2F( "B_intersection_primaryTrackMom2pos_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2900 
2901  mHistograms[ "intersection_primaryTrackneg_globalXY" ] = new TH2F( "B_intersection_primaryTrackneg_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2902  mHistograms[ "intersection_primaryTrackMom0neg_globalXY" ] = new TH2F( "B_intersection_primaryTrackMom0neg_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2903  mHistograms[ "intersection_primaryTrackMom1neg_globalXY" ] = new TH2F( "B_intersection_primaryTrackMom1neg_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2904  mHistograms[ "intersection_primaryTrackMom2neg_globalXY" ] = new TH2F( "B_intersection_primaryTrackMom2neg_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2905 
2906 
2907 
2908  mHistograms[ "intersection_track_pt_eta" ] = new TH2F( "B_intersection_track_pt_eta", "eta vs. pt;p_{T} (GeV/c);#eta", 400, 0., 5., 400, -2., -0.7 );
2909  mHistograms[ "intersection_track_pt_phi" ] = new TH2F( "B_intersection_track_pt_phi", "phi vs. pt;p_{T} (GeV/c);#phi", 400, 0., 5., 400, 0., 2 * M_PI );
2910  mHistograms[ "intersection_track_phi_eta" ] = new TH2F( "B_intersection_track_phi_eta", "eta vs. phi;#phi;#eta", 400, 0., 2 * M_PI, 400, -2., -0.9 );
2911 
2912  mHistograms[ "intersection_track_nHitsTpc" ] = new TH1F( "B_intersection_track_nHitsTpc", "nHitsTpc;nHitsFit;# tracks", 75, 0., 75. );
2913 
2914  mHistograms[ "intersection_track_mom_dEdx" ] = new TH2F( "B_intersection_track_mom_dEdx", "dE/dx vs. mom;mom (GeV/c);dE/dx (keV/cm)", 100, 0., 5., 100, 0., 10. );
2915  mHistograms[ "intersection_track_mom_nsigmaPi" ] = new TH2F( "B_intersection_track_mom_nsigmaPi", "n#sigma_{#pi} vs. mom; mom (GeV/c);n#sigma_{#pi}", 100, 0., 5., 100, -10., 10. );
2916 
2917 
2918  // ----------
2919  // step - C -
2920  // ----------
2921  mHistograms[ "detHitvsInter_X" ] = new TH2F( "C_detHitvsInter_X" , "detectorHit vs. intersection X;detectorHit X (cm);intersection X (cm)", 400, -300., 300., 400, -300., 300.);
2922  mHistograms[ "detHitvsInter_Y" ] = new TH2F( "C_detHitvsInter_Y" , "detectorHit vs. intersection Y;detectorHit Y (cm);intersection Y (cm)", 400, -300., 300., 400, -300., 300.);
2923 
2924  mHistograms[ "detHitvsInter_localX" ] = new TH2F( "C_detHitvsInter_localX" , "detectorHit vs. intersection X;detectorHit local X (cm);intersection local X (cm)", 40, -20., 20., 80, -20., 20.);
2925  mHistograms[ "detHitvsInter_localY" ] = new TH2F( "C_detHitvsInter_localY" , "detectorHit vs. intersection Y;detectorHit local Y (cm);intersection local Y (cm)", 200, -50., 50., 80, -20., 20.);
2926 
2927 
2928  for( int sector = eTofConst::sectorStart; sector <= eTofConst::sectorStop; sector++ ) {
2929  for( int plane = eTofConst::zPlaneStart; plane <= eTofConst::zPlaneStop; plane++ ) {
2930  std::string histName_detHitvsInter_strip = "detHitvsInter_strip_s" + std::to_string( sector ) + "m" + std::to_string( plane );
2931 
2932  mHistograms[ histName_detHitvsInter_strip ] = new TH2F( Form( "C_detHitvsInter_strip_s%dm%d", sector, plane ), Form( "detectorHit vs. intersection on sector %d module %d;detectorHit strip;intersection strip", sector, plane ), 96, -0.5, 95.5, 96, -0.5, 95.5 );
2933  }
2934  }
2935 
2936  mHistograms[ "moduleIndex_deltaX" ] = new TH2F( "C_moduleIndex_deltaX", "module index vs. local #Delta X;module index;#DeltaX (cm)", 36, -0.5, 35.5, 100, -50., 50. );
2937  mHistograms[ "moduleIndex_deltaY" ] = new TH2F( "C_moduleIndex_deltaY", "module index vs. local #Delta Y;module index;#DeltaY (cm)", 36, -0.5, 35.5, 100, -50., 50. );
2938 
2939  mHistograms[ "matchCand_phi_eta" ] = new TH2F( "C_matchCand_phi_eta", "eta vs. phi; #phi; #eta", 400, 0., 2 * M_PI, 400, -1.7, -0.9 );
2940 
2941  mHistograms[ "matchCand_deltaX" ] = new TH1F( "C_matchCand_deltaX" , "match candidate delta X;match candidate #DeltaX (cm); #match cand", 400, -15., 15. );
2942  mHistograms[ "matchCand_deltaY" ] = new TH1F( "C_matchCand_deltaY" , "match candidate delta Y;match candidate #DeltaY (cm); #match cand", 400, -15., 15. );
2943 
2944 
2945  for( int sector = eTofConst::sectorStart; sector <= eTofConst::sectorStop; sector++ ) {
2946  for( int plane = eTofConst::zPlaneStart; plane <= eTofConst::zPlaneStop; plane++ ) {
2947  for( int counter = eTofConst::counterStart; counter <= eTofConst::counterStop; counter++ ) {
2948  std::string histName_deltaX = "matchCand_deltaX_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2949  std::string histName_deltaY = "matchCand_deltaY_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2950 
2951  mHistograms[ histName_deltaX ] = new TH1F( Form( "C_matchCand_deltaX_s%dm%dc%d", sector, plane, counter ) , Form( "match candidate delta X in sector %d module %d counter %d;match candidate #DeltaX (cm); #match cand", sector, plane, counter ), 400, -15., 15. );
2952  mHistograms[ histName_deltaY ] = new TH1F( Form( "C_matchCand_deltaY_s%dm%dc%d", sector, plane, counter ) , Form( "match candidate delta Y in sector %d module %d counter %d;match candidate #DeltaY (cm); #match cand", sector, plane, counter ), 400, -15., 15. );
2953  }
2954  }
2955  }
2956 
2957  mHistograms[ "matchCand_deltaX_nHitsTpc" ] = new TH2F( "C_matchCand_deltaX_nHitsTpc" , "match candidate delta X vs. nHitsFit in TPC;nHitsFit in TPC;match candidate #DeltaX (cm)", 75, 0., 75., 400, -15., 15. );
2958  mHistograms[ "matchCand_deltaY_nHitsTpc" ] = new TH2F( "C_matchCand_deltaY_nHitsTpc" , "match candidate delta Y vs. nHitsFit in TPC;nHitsFit in TPC;match candidate #DeltaY (cm)", 75, 0., 75., 400, -15., 15. );
2959 
2960  mHistograms[ "matchCandMult" ] = new TH1F( "C_matchCandMult", "matchCandMult;multiplicity;# events", 200, 0., 200. );
2961 
2962 
2963  // ----------
2964  // step - D -
2965  // ----------
2966  mHistograms[ "trackMatchMultPerDetectorHit" ] = new TH1F( "D_trackMatchMultPerDetectorHit", "multiplicity of tracks pointing to the same detector hit;#tracks;#detector hits", 15, 0., 15. );
2967 
2968  mHistograms[ "singleTrackMatchMult" ] = new TH1F( "D_singleTrackMatchMult", "singleTrackMatchMult;multiplicity;# events", 200, 0., 200. );
2969 
2970 
2971  // ----------
2972  // step - E -
2973  // ----------
2974  mHistograms[ "hitMultPerTrack" ] = new TH1F( "E_hitMultPerTrack", "multiplicity of hit matched to the same track;#hits;#tracks", 15, 0., 15. );
2975 
2976  mHistograms[ "finalMatch_pt" ] = new TH1F( "E_finalMatch_pt", "p_{T} distribution of matched tracks", 200, 0., 2. );
2977 
2978  mHistograms[ "finalMatchMult" ] = new TH1F( "E_finalMatchMult", "finalMatchMult;multiplicity;# events", 200, 0., 200. );
2979 
2980  mHistograms[ "overlapHit_globalXY" ] = new TH2F( "E_overlapHit_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2981 
2982 
2983  // ----------
2984  // step - F -
2985  // ----------
2986  mHistograms[ "finalMatchMultGlobal" ] = new TH1F( "F_finalMatchMultGlobal", "finalMatchMultGlobal;multiplicity;# events", 200, 0., 200. );
2987  mHistograms[ "finalMatchMultPrimary" ] = new TH1F( "F_finalMatchMultPrimary", "finalMatchMultPrimary;multiplicity;# events", 200, 0., 200. );
2988 
2989  mHistograms[ "finalMatchPrimary_globalXY" ] = new TH2F( "F_finalMatchPrimary_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2990  mHistograms[ "finalMatchPrimaryMom0_globalXY" ] = new TH2F( "F_finalMatchPrimaryMom0_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2991  mHistograms[ "finalMatchPrimaryMom1_globalXY" ] = new TH2F( "F_finalMatchPrimaryMom1_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2992  mHistograms[ "finalMatchPrimaryMom2_globalXY" ] = new TH2F( "F_finalMatchPrimaryMom2_globalXY", "global XY;x (cm);y (cm)", 400, -300., 300., 400, -300., 300. );
2993 
2994 
2995  // ----------
2996  // step - G -
2997  // ----------
2998  mHistograms[ "startTimeDiff" ] = new TH1F( "G_startTimeDiff", "bTOF - eTOF start time;#DeltaT;#events", 1000, -20., 20. );
2999  mHistograms[ "startTimeDiff_nCand" ] = new TH2F( "G_startTimeDiff_nCand", "bTOF - eTOF start time;#DeltaT;#nCand tracks", 1000, -20., 20., 50, 0, 50 );
3000  mHistograms[ "startTime_T0corr" ] = new TH1F( "G_startTime_T0corr", "T0corr evolution;#updates;T0 corr. (ns)", 1000, 0., 1000. );
3001 
3002  mHistograms[ "matchCand_timeOfFlight" ] = new TH1F( "G_matchCand_timeOfFlight", "match candidate time of flight;ToF (ns);# match candidates", 2000, -400., 600. );
3003  mHistograms[ "matchCand_timeOfFlight_pathLength" ] = new TH2F( "G_matchCand_timeOfFlight_pathLength", "match candidate pathlength vs. time of flight;ToF (ns);pathlength (cm)", 1000, -400., 600., 800, 200., 600. );
3004 
3005 
3006  // ----------
3007  // step - H -
3008  // ----------
3009  mHistograms[ "matchCand_beta_mom" ] = new TH2F( "H_matchCand_beta_mom" , "match candidate 1/beta vs. momentum;p (GeV/c);1/#beta", 400, 0., 10., 1000, 0.0, 2. );
3010 
3011  mHistograms[ "matchCand_beta_mom_matchDistCut" ] = new TH2F( "H_matchCand_beta_mom_matchDistCut" , "match candidate 1/beta vs. momentum;p (GeV/c);1/#beta", 400, 0., 10., 1000, 0.8, 2. );
3012 
3013  mHistograms[ "matchCandMultPerSector_matchDistCut" ] = new TH2F( "H_matchCandMultPerSector_matchDistCut" , "matchCandMultPerSector_matchDistCut;module;# matches per event", 36, 0, 36, 25, 0, 25 );
3014 
3015 
3016  mHistograms[ "matchCand_m2_signmom" ] = new TH2F( "H_matchCand_m2_signmom" , "match candidate m^{2} vs. momentum;q/|q| * p (GeV/c);m^{2} (GeV^{2}/c^{4})", 400, -10., 10., 1000, -0.2, 1.3 );
3017  mHistograms[ "matchCand_m2_mom" ] = new TH2F( "H_matchCand_m2_mom" , "match candidate m^{2} vs. momentum;p (GeV/c);m^{2} (GeV^{2}/c^{4})", 400, 0., 10., 1000, -0.2, 1.3 );
3018 
3019  mHistograms[ "matchCand_m2_mom_matchDistCut" ] = new TH2F( "H_matchCand_m2_mom_matchDistCut" , "match candidate m^{2} vs. momentum;p (GeV/c);m^{2} (GeV^{2}/c^{4})", 400, 0., 10., 1000, -0.2, 1.3 );
3020 
3021 
3022  for( int sector = eTofConst::sectorStart; sector <= eTofConst::sectorStop; sector++ ) {
3023  for( int plane = eTofConst::zPlaneStart; plane <= eTofConst::zPlaneStop; plane++ ) {
3024  for( int counter = eTofConst::counterStart; counter <= eTofConst::counterStop; counter++ ) {
3025 
3026  std::string histName_beta_mom = "matchCand_beta_mom_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
3027  mHistograms[ histName_beta_mom ] = new TH2F( Form( "H_matchCand_beta_mom_s%dm%dc%d", sector, plane, counter ), Form( "match candidate 1/beta vs. momentum in sector %d module %d counter %d;p (GeV/c);1/#beta", sector, plane, counter), 200, 0., 10., 500, 0.0, 2. );
3028 
3029  std::string histName_t0corr_mom = "matchCand_t0corr_mom_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
3030  mHistograms[ histName_t0corr_mom ] = new TH2F( Form( "H_matchCand_t0corr_mom_s%dm%dc%d", sector, plane, counter ), Form( "measured tof - tof_{#pi} vs. momentum in sector %d module %d counter %d;mom (GeV/c);#Delta time (ns)", sector, plane, counter ), 400, 0., 10., 1000, -500., 500. );
3031 
3032  std::string histName_t0corr_mom_zoom = "matchCand_t0corr_mom_zoom_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
3033  std::string histName_t0corr_mom_zoom_cut = "matchCand_t0corr_mom_zoom_cut_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
3034 
3035  mHistograms[ histName_t0corr_mom_zoom ] = new TH2F( Form( "H_matchCand_t0corr_mom_zoom_s%dm%dc%d", sector, plane, counter ), Form( "measured tof - tof_{#pi} vs. momentum in sector %d module %d counter %d;mom (GeV/c);#Delta time (ns)", sector, plane, counter ), 200, 0., 3., 2000, -10., 10. );
3036  mHistograms[ histName_t0corr_mom_zoom_cut ] = new TH2F( Form( "H_matchCand_t0corr_mom_zoom_cut_s%dm%dc%d", sector, plane, counter ), Form( "measured tof - tof_{#pi} vs. momentum in sector %d module %d counter %d;mom (GeV/c);#Delta time (ns)", sector, plane, counter ), 200, 0., 3., 2000, -10., 10. );
3037 
3038  std::string histName_t0corr_strip = "matchCand_t0corr_strip_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
3039  mHistograms[ histName_t0corr_strip ] = new TH2F( Form( "H_matchCand_t0corr_strip_s%dm%dc%d", sector, plane, counter ), Form( "measured tof - tof_{#pi} vs. momentum in sector %d module %d counter %d;local X (cm);#Delta time (ns)", sector, plane, counter ), 32, -16., 16., 2000, -10., 10. );
3040 
3041  std::string histName_slewing_digi = "matchCand_slewing_digi_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
3042  mHistograms[ histName_slewing_digi ] = new TH2F( Form( "H_matchCand_slewing_digi_s%dm%dc%d", sector, plane, counter ), Form( "measured tof - tof_{#pi} vs. digi ToT in sector %d module %d counter %d;digi Tot (arb. units);#Delta time (ns)", sector, plane, counter ), 400, 0., 20., 1000, -3., 3. );
3043  }
3044  }
3045  }
3046 
3047  mHistograms[ "matchCand_dEdx_beta_mom1gev" ] = new TH2F( "H_matchCand_dEdx_beta_mom1gev", "dE/dx vs. 1/#beta at p=1 GeV;1/#beta;dE/dx (keV/cm)", 800, 0.8, 1.8, 800, 0., 15. );
3048  mHistograms[ "matchCand_dEdx_beta_mom2gev" ] = new TH2F( "H_matchCand_dEdx_beta_mom2gev", "dE/dx vs. 1/#beta at p=2 GeV;1/#beta;dE/dx (keV/cm)", 800, 0.8, 1.8, 800, 0., 15. );
3049  }
3050 
3051 
3052  for( auto& kv : mHistograms ) {
3053  kv.second->SetDirectory( 0 );
3054  }
3055  for( auto& kv : mHistograms2d ) {
3056  kv.second->SetDirectory( 0 );
3057  }
3058 }
3059 
3060 
3061 //---------------------------------------------------------------------------
3062 void
3063 StETofMatchMaker::writeHistograms()
3064 {
3065  if( mHistFileName != "" ) {
3066  LOG_INFO << "writing histograms to: " << mHistFileName.c_str() << endm;
3067 
3068  TFile histFile( mHistFileName.c_str(), "RECREATE", "etofMatch" );
3069  histFile.cd();
3070 
3071  for ( const auto& kv : mHistograms ) {
3072  if( kv.second->GetEntries() > 0 ) kv.second->Write();
3073  }
3074  for ( const auto& kv : mHistograms2d ) {
3075  if( kv.second->GetEntries() > 0 ) kv.second->Write();
3076  }
3077  histFile.Close();
3078  }
3079  else {
3080  LOG_INFO << "histogram file name is empty string --> cannot write histograms" << endm;
3081  }
3082 }
3083 
3084 
3085 //---------------------------------------------------------------------------
3086 void
3087 StETofMatchMaker::setMatchDistXYT( double x, double y, double t = 99999. )
3088 {
3089  mMatchDistX = fabs( x );
3090  mMatchDistY = fabs( y );
3091  mMatchDistT = fabs( t );
3092 }
3093 
3094 //---------------------------------------------------------------------------
3095 int
3096 StETofMatchMaker::rotateHit( const int& sector, const int& rot )
3097 {
3098  int sec = 13 + ( sector - 13 + rot ) % 12;
3099  LOG_DEBUG << "hit rotated from: " << sector << " to " << sec << endm;
3100  return sec;
3101 }
3102 
3103 
3104 //---------------------------------------------------------------------------
3105 ETofTrack::ETofTrack( const StTrack* sttrack )
3106 {
3107  mom = -999.;
3108  pt = -999.;
3109  eta = -999.;
3110  phi = -999.;
3111  nFtPts = 0;
3112  nDedxPts = 0;
3113  flag = 0;
3114  nHitsPoss = 999;
3115  dEdx = -999.;
3116  nSigmaPion = -999.;
3117 
3118  if( sttrack ) {
3119  mom = sttrack->geometry()->momentum().mag();
3120  pt = sttrack->geometry()->momentum().perp();
3121  eta = sttrack->geometry()->momentum().pseudoRapidity();
3122  phi = sttrack->geometry()->momentum().phi();
3123  nFtPts = sttrack->fitTraits().numberOfFitPoints( kTpcId );
3124 
3125  static StTpcDedxPidAlgorithm PidAlgorithm;
3126  static StPionPlus* Pion = StPionPlus::instance();
3127  const StParticleDefinition* pd = sttrack->pidTraits( PidAlgorithm );
3128  if( pd ){
3129  nSigmaPion = PidAlgorithm.numberOfSigma( Pion );
3130  if(PidAlgorithm.traits()){
3131  dEdx = PidAlgorithm.traits()->mean() * 1.e6;
3132  nDedxPts = PidAlgorithm.traits()->numberOfPoints();
3133  }
3134  }
3135  flag = sttrack->flag();
3136  nHitsPoss = sttrack->numberOfPossiblePoints( kTpcId );
3137 
3138  if ( phi < 0. ) phi += 2. * M_PI;
3139  }
3140 }
3141 
3142 //---------------------------------------------------------------------------
3143 ETofTrack::ETofTrack( const StMuTrack* mutrack )
3144 {
3145  mom = -999.;
3146  pt = -999.;
3147  eta = -999.;
3148  phi = -999.;
3149  nFtPts = 0;
3150  nDedxPts = 0;
3151  flag = 0;
3152  nHitsPoss = 999;
3153  dEdx = -999.;
3154  nSigmaPion = -999.;
3155 
3156  if( mutrack ) {
3157  mom = mutrack->momentum().mag();
3158  pt = mutrack->momentum().perp();
3159  eta = mutrack->momentum().pseudoRapidity();
3160  phi = mutrack->momentum().phi();
3161  nFtPts = mutrack->nHitsFit( kTpcId );
3162  nDedxPts = mutrack->nHitsDedx();
3163  flag = mutrack->flag();
3164  nHitsPoss = mutrack->nHitsPoss( kTpcId );
3165  dEdx = mutrack->dEdx() * 1.e6;
3166  nSigmaPion = mutrack->nSigmaPion();
3167 
3168  if ( phi < 0. ) phi += 2. * M_PI;
3169  }
3170 }
3171 
3172 //---------------------------------------------------------------------------
3173 void StETofMatchMaker::checkClockJumps()
3174 {
3175  if (mClockJumpCand.size() == 0) return;
3176 
3177  // histogram filled with all hits including clock jump candidates.
3178  int binmax = mHistograms.at("matchCand_t0corr_1d")->GetMaximumBin();
3179  float mainPeakT0 = mHistograms.at("matchCand_t0corr_1d")->GetXaxis()->GetBinCenter(binmax);
3180 
3181  LOG_DEBUG << "mainPeakT0: " << mainPeakT0 << endm;
3182 
3183  bool needsUpdate = false;
3184 
3185  for (const auto& kv : mClockJumpCand) {
3186  LOG_DEBUG << "clock jump candidate: " << kv.first << " " << kv.second << endm;
3187 
3188  int sector = kv.first / 1000;
3189  int plane = (kv.first % 1000) / 100;
3190  int counter = (kv.first % 100) / 10;
3191  int binX = kv.first % 10;
3192 
3193  std::string histName_t0corr_jump = "matchCand_t0corr_jump_s" + std::to_string(sector) + "m" + std::to_string(plane) + "c" + std::to_string(counter);
3194  TH2D* h = (TH2D*) mHistograms.at(histName_t0corr_jump);
3195 
3196  int binYmain_neg = h->GetYaxis()->FindBin(mainPeakT0 - 1.5);
3197  int binYmain_pos = h->GetYaxis()->FindBin(mainPeakT0 + 3.0);
3198  int binYearly_neg = h->GetYaxis()->FindBin(mainPeakT0 - 1.5 - eTofConst::coarseClockCycle);
3199  int binYearly_pos = h->GetYaxis()->FindBin(mainPeakT0 + 3.0 - eTofConst::coarseClockCycle);
3200  int binYlate_neg = h->GetYaxis()->FindBin(mainPeakT0 - 1.5 + eTofConst::coarseClockCycle); // tight cut to reduce interference with slow particles of the main band
3201  int binYlate_pos = h->GetYaxis()->FindBin(mainPeakT0 + 3.0 + eTofConst::coarseClockCycle);
3202 
3203  LOG_DEBUG << "binYmain_neg " << binYmain_neg << " " << binYmain_pos << " binYmain_pos " << binYearly_neg << " binYearly_neg " << binYearly_pos << " binYearly_pos " << endm;
3204 
3205  int nMain = h->Integral(binX, binX, binYmain_neg, binYmain_pos);
3206  int nEarly = h->Integral(binX, binX, binYearly_neg, binYearly_pos);
3207  int nLate = h->Integral(binX, binX, binYlate_neg, binYlate_pos);
3208 
3209  LOG_DEBUG << "nMain " << nMain << " " << nEarly << " nEarly " << endm;
3210 
3211  if (nEarly > nMain && nEarly >= 2) { // first two clock jumped hits in one Get4 IN EACH FILE are not are not detected. Could be changed by changing default direction of jumps, but then wrongly corrected clock jumps are LATE! Cut on late hits being above 1.5 GeV or Pion dEdX could help separate late hits from slow particles.
3212  LOG_DEBUG << "clock jump detected --> give it to hit maker" << endm;
3213 
3214  mClockJumpDirection[ kv.first ] = -1.;
3215  needsUpdate = true;
3216  }
3217 
3218  if (nLate > nMain && nLate >= 2) { //allows to change default to correct backwards in time, as it is electronically most likely. Unlikely to find real proton at wrong position, but correct time.
3219  LOG_DEBUG << "clock jump detected --> give it to hit maker" << endm;
3220 
3221  mClockJumpDirection[ kv.first ] = 1.;
3222  needsUpdate = true;
3223  }
3224  }
3225 
3226  mClockJumpCand.clear();
3227 
3228  // if there was a new entry to the map --> push it to the hit maker (if available)
3229  if( needsUpdate && mETofHitMaker ) {
3230  mETofHitMaker->updateClockJumpMap( mClockJumpDirection );
3231  }
3232 }
3233 
3234 //---------------------------------------------------------------------------
3235 void
3236 StETofMatchMaker::sortMatchCases( eTofHitVec inputVec , std::map< Int_t, eTofHitVec >& outputMap )
3237 {
3238 
3239 
3240  // sort & flag Match candidates
3241 
3242  // define temporary vectors for iterating through matchCandVec
3243  eTofHitVec tempVec = inputVec;
3244  eTofHitVec erasedVec = tempVec;
3245  eTofHitVec tempMMVec;
3246  tempMMVec.clear();
3247  std::map< Int_t, eTofHitVec > MMMap;
3248  MMMap.clear();
3249 
3250  eTofHitVec ssVec;
3251 
3252  // get multi Hit sets
3253  // int deltaSize = 0;
3254 
3255  eTofHitVecIter tempIter = tempVec.begin();
3256  eTofHitVecIter erasedIter = erasedVec.begin();
3257 
3258 
3259  if(tempVec.size() < 1 ) return;
3260 
3261  while( tempVec.size() != 0 ) {
3262 
3263  tempIter = tempVec.begin();
3264  erasedIter = erasedVec.begin();
3265 
3266 
3267  tempMMVec.push_back(*tempIter);
3268  erasedVec.erase( erasedIter );
3269 
3270  // int sizeOld = tempMMVec.size();
3271  int count =0;
3272  int countwhile = 0;
3273 
3274  for(unsigned int s =0; s < tempMMVec.size(); s++){
3275 
3276  count++;
3277 
3278  erasedIter = erasedVec.begin();
3279 
3280  if(erasedVec.size() <= 0 ) continue;
3281 
3282  countwhile = 0;
3283 
3284  while( erasedIter != erasedVec.end() ) {
3285 
3286  countwhile++;
3287 
3288  if(tempMMVec.at(s).trackId == erasedIter->trackId && tempMMVec.at(s).index2ETofHit == erasedIter->index2ETofHit){
3289 
3290  erasedVec.erase( erasedIter );
3291  erasedIter++;
3292  continue;}
3293  if(tempMMVec.at(s).trackId == erasedIter->trackId || tempMMVec.at(s).index2ETofHit == erasedIter->index2ETofHit){
3294  if(!mIsSim){
3295  // erasedIter->matchFlag = 0;
3296  }
3297  tempMMVec.push_back(*erasedIter);
3298 
3299  erasedVec.erase( erasedIter );
3300 
3301  }
3302  if( erasedVec.size() <= 0 ) break;
3303  if( erasedIter == erasedVec.end()) break;
3304  erasedIter++;
3305 
3306  } //while inner
3307 
3308  }// for
3309 
3310  // deltaSize = sizeOld - tempMMVec.size();
3311 
3312  MMMap[tempMMVec.begin()->trackId] = tempMMVec;
3313  tempMMVec.clear();
3314 
3315  tempVec = erasedVec;
3316  }
3317 
3318 
3319  outputMap = MMMap;
3320 
3321 
3322 }
3323 //---------------------------------------------------------------------------
3324 void
3325 StETofMatchMaker::sortandcluster(eTofHitVec& matchCandVec , eTofHitVec& detectorHitVec , eTofHitVec& intersectionVec , eTofHitVec& finalMatchVec){
3326 
3327 
3328  //flag Overlap-Hits -------------------------------------------------------------------
3329  std::map< Int_t, eTofHitVec > overlapHitMap;
3330  eTofHitVec overlapHitVec;
3331  eTofHitVec tempVecOL = matchCandVec;
3332  eTofHitVecIter detHitIter;
3333  eTofHitVecIter detHitIter2;
3334 
3335  for( auto detHitIter = tempVecOL.begin(); detHitIter != tempVecOL.end(); ) {
3336 
3337  detHitIter = tempVecOL.begin();
3338  detHitIter2 = tempVecOL.begin();
3339 
3340  bool isOverlap = false;
3341  int counterId1 = (detHitIter->sector*100) + (detHitIter->plane*10) + (detHitIter->counter);
3342 
3343  for( auto detHitIter2 = tempVecOL.begin(); detHitIter2 != tempVecOL.end(); ) {
3344 
3345  int counterId2 = (detHitIter2->sector*100) + (detHitIter2->plane*10) + (detHitIter2->counter);
3346 
3347  if(counterId1 != counterId2 && detHitIter->trackId == detHitIter2->trackId){
3348 
3349  int mf2 = counterId2 ;
3350 
3351  detHitIter2->matchFlag = mf2;
3352 
3353  matchCandVec.at(detHitIter2 - tempVecOL.begin()).matchFlag = 1;
3354 
3355  overlapHitVec.push_back(*detHitIter2);
3356  tempVecOL.erase(detHitIter2);
3357 
3358  isOverlap = true;
3359  }
3360 
3361  if( tempVecOL.size() <= 0 ) break;
3362  if( detHitIter2 == tempVecOL.end()) break;
3363  detHitIter2++;
3364 
3365  }
3366 
3367  if(isOverlap){
3368 
3369  detHitIter->matchFlag = counterId1;
3370 
3371  matchCandVec.at(detHitIter - tempVecOL.begin()).matchFlag = 1;
3372 
3373  overlapHitVec.push_back(*detHitIter);
3374 
3375  //fill map
3376  overlapHitMap[overlapHitVec.begin()->trackId] = overlapHitVec;
3377 
3378  overlapHitVec.clear();
3379 
3380  }
3381  tempVecOL.erase(detHitIter);
3382 
3383  if( tempVecOL.size() <= 0 ) break;
3384  if( detHitIter == tempVecOL.end()) break;
3385  detHitIter++;
3386 
3387  }
3388 
3389  // fill match cand vec counter wise
3390  std::vector< eTofHitVec > matchCandVecCounter(108);
3391 
3392  for(int i =0; i < 108; i++){
3393 
3394  for(unsigned int n = 0; n < matchCandVec.size(); n++){
3395 
3396  int sector = matchCandVec.at(n).sector;
3397  int plane = matchCandVec.at(n).plane;
3398  int counter = matchCandVec.at(n).counter;
3399 
3400  int counterId = 9*(sector - 13) + 3*(plane - 1) + (counter -1);
3401 
3402  if(counterId == i ) {
3403  matchCandVecCounter.at(i).push_back(matchCandVec.at(n));
3404  }
3405  }//loop over hits
3406  }//loop over counters
3407 
3408 
3409  // loop over counters
3410  for(int counterNr = 0; counterNr < 108; counterNr++){
3411 
3412  // sort & flag Match candidates
3413  eTofHitVec tempVec = matchCandVecCounter.at(counterNr);
3414  eTofHitVec tempVec2 = matchCandVecCounter.at(counterNr);
3415  std::map< Int_t, eTofHitVec > MMMap;
3416 
3417  sortMatchCases(tempVec, MMMap);
3418 
3419  // final containers
3420  std::map< Int_t, eTofHitVec > MultMultMap;
3421  std::map< Int_t, eTofHitVec > SingleHitMap;
3422  std::map< Int_t, eTofHitVec > SingleTrackMap;
3423  eTofHitVec ssVec;
3424 
3425  map<Int_t, eTofHitVec >::iterator it;
3426 
3427  for (it = MMMap.begin(); it != MMMap.end(); it++)
3428  {
3429  int nTracks = 1;
3430  int nHits = 1;
3431 
3432  for(unsigned int l =0; l< it->second.size(); l++){
3433  for(unsigned int j = l; j< it->second.size(); j++){
3434 
3435  if( it->second.at(l).trackId != it->second.at(j).trackId) nTracks++;
3436  if( it->second.at(l).index2ETofHit != it->second.at(j).index2ETofHit ) nHits++;
3437 
3438  } // for inner
3439  } //for outer
3440 
3441 
3442  // cases::
3443  // Single Hit - Single Track
3444  if(nTracks == 1 && nHits == 1) {
3445 
3446  ssVec.push_back(it->second.front() );
3447 
3448  int isMerged = 10; // 10 codes for normal hit
3449  int isOl = it->second.front().matchFlag;
3450 
3451  if( it->second.front().clusterSize > 999 ) {
3452 
3453  isMerged = 20; // 20 codes for single hit
3454  it->second.front().clusterSize = 1;
3455  }
3456 
3457  it->second.front().matchFlag = 100 + isMerged + isOl;
3458  finalMatchVec.push_back(it->second.front());
3459  }
3460 
3461 
3462  // Single Hit - Multi Track
3463  if( nTracks > 1 && nHits == 1) {
3464 
3465  double dr = 0.0;
3466  double dr_best = 99999.0; // dy for SHs at 27
3467  unsigned int ind = 0;
3468  unsigned int ind_best = 0;
3469 
3470  for(unsigned int l =0; l < it->second.size(); l++){
3471 
3472  dr = (it->second.at(l).deltaX * it->second.at(l).deltaX) + (it->second.at(l).deltaY * it->second.at(l).deltaY);
3473  ind = l;
3474 
3475  if(dr <= dr_best){
3476  dr_best = dr;
3477  ind_best = ind;
3478  }
3479  }
3480 
3481  SingleHitMap[it->first] = it->second;
3482 
3483  //pick closest track and push to finalMatchVec
3484  int isMerged = 10;
3485  int isOl = it->second.at(ind_best).matchFlag;
3486 
3487  if( it->second.at(ind_best).clusterSize > 999 ){
3488 
3489  isMerged = 20;
3490  it->second.at(ind_best).clusterSize = 1;
3491  }
3492 
3493  it->second.at(ind_best).matchFlag = 300 + isMerged + isOl;
3494  finalMatchVec.push_back(it->second.at(ind_best));
3495  }
3496 
3497 
3498  // Multi Hit - Single Track
3499  if( nTracks ==1 && nHits > 1) {
3500 
3501  bool isN = false;
3502  bool isS = false;
3503 
3504  for(unsigned int l =0; l < it->second.size(); l++){
3505 
3506  if(it->second.at(l).clusterSize < 999){
3507  isN = true;
3508  }else{
3509  isS = true;
3510  }
3511  }
3512 
3513  SingleTrackMap[it->first] = it->second;
3514 
3515  // sort by merge cases :: SS, NN, SN
3516  //NN
3517  if(isN && (!isS)){
3518 
3519  std::vector< std::vector<int> > mergeIndVec(it->second.size());
3520 
3521  double dr_sum=0;
3522  double dr_diff = 0;
3523  double dr_mean=0;
3524 
3525  double dr = 0.0;
3526  double dr_best = 99999.0; // dy for SHs at 27
3527  unsigned int ind = 0;
3528  unsigned int ind_best = 0;
3529 
3530  for(unsigned int l =0; l < it->second.size(); l++){
3531 
3532  dr = sqrt((it->second.at(l).deltaX * it->second.at(l).deltaX) + (it->second.at(l).deltaY * it->second.at(l).deltaY));
3533  ind = l;
3534 
3535  dr_sum += abs(dr);
3536 
3537  if(dr <= dr_best){
3538  dr_best = dr;
3539  ind_best = ind;
3540  }
3541  }
3542 
3543  dr_mean = dr_sum / it->second.size();
3544 
3545  for(unsigned int c =0; c < it->second.size(); c++){
3546 
3547  dr = sqrt((it->second.at(c).deltaX * it->second.at(c).deltaX) + (it->second.at(c).deltaY * it->second.at(c).deltaY));
3548  dr_diff += abs(dr - dr_mean);
3549  }
3550 
3551  // NN Hits already merged in HitMaker
3552 
3553  int mergedCluSz = 0;
3554  int mergedMatchFlag = 0;
3555  int isMerged = 0;
3556  int isOl = it->second.at(ind_best).matchFlag;
3557 
3558  if(it->second.at(ind_best).clusterSize > 100 && it->second.at(ind_best).clusterSize < 200){
3559 
3560  mergedCluSz = it->second.at(ind_best).clusterSize % 100;
3561 
3562  }else{
3563 
3564  mergedCluSz = it->second.at(ind_best).clusterSize;
3565  }
3566 
3567  if(mergedCluSz > 1){ isMerged = 30; // 30 codes for normal-normal-merge
3568  }else{
3569  isMerged = 10;
3570  }
3571 
3572  mergedMatchFlag = 200 + isMerged + isOl; // 200 codes for SingleTrackMultiHit
3573  it->second.at(ind_best).matchFlag = mergedMatchFlag;
3574 
3575  finalMatchVec.push_back(it->second.at(ind_best));
3576  }
3577 
3578  //SS
3579  if(isS && (!isN)){
3580 
3581  std::vector< std::vector<int> > mergeIndVec(it->second.size());
3582 
3583  for(unsigned int l =0; l < it->second.size(); l++){
3584  mergeIndVec.at(l).push_back(0);
3585  }
3586 
3587  double dr = 0.0;
3588  double dr_best = 99999.0; // dy for SHs at 27
3589  unsigned int ind = 0;
3590  unsigned int ind_best = 0;
3591 
3592  for(unsigned int l =0; l < it->second.size(); l++){
3593 
3594  // localY doesnt contain any ETOF information -> not usefull for merging single sided hits
3595  dr = it->second.at(l).deltaX;
3596  ind = l;
3597 
3598  if(dr <= dr_best){
3599  dr_best = dr;
3600  ind_best = ind;
3601  }
3602  }
3603 
3604  // merge MatchCands
3605 
3606  eTofHitVec hitVec = it->second ;
3607 
3608  double mergedTime = it->second.at(ind_best).hitTime;
3609  double mergedToT = it->second.at(ind_best).tot;
3610  double mergedPosY = it->second.at(ind_best).localY;
3611  double mergedPosX = it->second.at(ind_best).localX;
3612  int mergedCluSz = 1;
3613  int mergedMatchFlag = 0;
3614  int mergedIdTruth = it->second.at(ind_best).IdTruth;
3615 
3616  for(unsigned int j=0; j < hitVec.size(); j++) {
3617 
3618  if( j == ind_best) continue;
3619 
3620  double dx = it->second.at(ind_best).localX - hitVec.at(j).localX;
3621  double dy = it->second.at(ind_best).localY - hitVec.at(j).localY;
3622  double dt = abs( it->second.at(ind_best).hitTime - hitVec.at(j).hitTime);
3623 
3624  // merge
3625  if( abs(dx) < dx_3sig && abs(dy) < dy_3sig && abs(dt) < dt_3sig ){
3626 
3627  mergedTime += hitVec.at(j).hitTime;
3628  mergedToT += hitVec.at(j).tot;
3629  mergedPosY += hitVec.at(j).localY;
3630  mergedPosX += hitVec.at(j).localX;
3631  mergedCluSz++;
3632 
3633  if(mergedIdTruth != hitVec.at(j).IdTruth) mergedIdTruth =0;
3634 
3635  }
3636  }
3637 
3638  // create mergend hit and MC;
3639  mergedTime /= mergedCluSz;
3640  mergedToT /= mergedCluSz;
3641  mergedPosY /= mergedCluSz;
3642  mergedPosX /= mergedCluSz;
3643  int isMerged = 0;
3644  int isOl = it->second.at(ind_best).matchFlag;
3645 
3646  if(mergedCluSz > 1){ isMerged = 40; // codes for sigle-single-merge
3647  }else{
3648  isMerged = 20;
3649  }
3650 
3651  mergedMatchFlag = 200 + isMerged + isOl; // 200 codes for SingleTrackMultiHit
3652 
3653  // use only the floating point remainder of the time with respect the the bTof clock range
3654  mergedTime = fmod( mergedTime, eTofConst::bTofClockCycle );
3655  if( mergedTime < 0 ) mergedTime += eTofConst::bTofClockCycle;
3656 
3657  it->second.at(ind_best).hitTime = mergedTime;
3658  it->second.at(ind_best).tot = mergedToT;
3659  it->second.at(ind_best).localX = mergedPosX;
3660  it->second.at(ind_best).localY = mergedPosY;
3661  it->second.at(ind_best).IdTruth = mergedIdTruth;
3662  it->second.at(ind_best).matchFlag = mergedMatchFlag;
3663  it->second.at(ind_best).clusterSize = mergedCluSz;
3664 
3665 
3666  finalMatchVec.push_back(it->second.at(ind_best));
3667  }
3668 
3669  //SN
3670  if(isN && isS){
3671 
3672  std::vector< std::vector<int> > mergeIndVec(it->second.size());
3673 
3674  for(unsigned int l =0; l < it->second.size(); l++){
3675  mergeIndVec.at(l).push_back(0);
3676  }
3677 
3678  double dr = 0.0;
3679  double dr_best = 99999.0; // dy for SHs at 27
3680  unsigned int ind = 0;
3681  unsigned int ind_best = 0;
3682 
3683  for(unsigned int l =0; l < it->second.size(); l++){
3684 
3685  if(it->second.at(l).clusterSize > 999) continue;
3686 
3687  // localY doesnt contain any ETOF information for singleSidedHits-> not usefull for merging later on
3688  dr = it->second.at(l).deltaX*it->second.at(l).deltaX + it->second.at(l).deltaY*it->second.at(l).deltaY;
3689  ind = l;
3690 
3691  if(dr <= dr_best){
3692  dr_best = dr;
3693  ind_best = ind;
3694  }
3695  }
3696 
3697 
3698  // merge MatchCands
3699  eTofHitVec hitVec = it->second ;
3700 
3701  double mergedTime = it->second.at(ind_best).hitTime;
3702  double mergedToT = it->second.at(ind_best).tot;
3703  double mergedPosY = it->second.at(ind_best).localY;
3704  double mergedPosX = it->second.at(ind_best).localX;
3705  int mergedCluSz = 1;
3706  int mergedMatchFlag = 0;
3707  int mergedIdTruth = it->second.at(ind_best).IdTruth;
3708 
3709  for(unsigned int j=0; j < hitVec.size(); j++) {
3710 
3711  if( j == ind_best) continue;
3712 
3713  double dx = it->second.at(ind_best).localX - hitVec.at(j).localX;
3714  double dy = it->second.at(ind_best).localY - hitVec.at(j).localY;
3715  double dt = abs( it->second.at(ind_best).hitTime - hitVec.at(j).hitTime);
3716 
3717  // merge
3718  if( abs(dx) < dx_3sig && abs(dy) < dy_3sig && abs(dt) < dt_3sig ){
3719 
3720  mergedTime += hitVec.at(j).hitTime;
3721  mergedToT += hitVec.at(j).tot;
3722  mergedPosY += hitVec.at(j).localY;
3723  mergedPosX += hitVec.at(j).localX;
3724  mergedCluSz++;
3725 
3726  if(mergedIdTruth != hitVec.at(j).IdTruth) mergedIdTruth =0;
3727  }
3728  }
3729 
3730  // create mergend hit and MC
3731  mergedTime /= mergedCluSz;
3732  mergedToT /= mergedCluSz;
3733  mergedPosY /= mergedCluSz;
3734  mergedPosX /= mergedCluSz;
3735  int isMerged = 0;
3736  int isOl = it->second.at(ind_best).matchFlag;
3737 
3738  if(mergedCluSz > 1){ isMerged = 50; // codes for sigle-normal-merge
3739  }else{
3740  isMerged = 10;
3741  }
3742 
3743  mergedMatchFlag = 200 + isMerged + isOl; // 200 codes for SingleTrackMultiHit
3744 
3745  // use only the floating point remainder of the time with respect the the bTof clock range
3746  mergedTime = fmod( mergedTime, eTofConst::bTofClockCycle );
3747  if( mergedTime < 0 ) mergedTime += eTofConst::bTofClockCycle;
3748 
3749  it->second.at(ind_best).hitTime = mergedTime;
3750  it->second.at(ind_best).tot = mergedToT;
3751  it->second.at(ind_best).localX = mergedPosX;
3752  it->second.at(ind_best).localY = mergedPosY;
3753  it->second.at(ind_best).IdTruth = mergedIdTruth;
3754  it->second.at(ind_best).matchFlag = mergedMatchFlag;
3755  it->second.at(ind_best).clusterSize = mergedCluSz ;
3756 
3757  finalMatchVec.push_back(it->second.at(ind_best));
3758  }
3759  } // multi-hit-single-track
3760 
3761 
3762  // Multi Hit - Multi Track
3763  if(nTracks > 1 && nHits > 1) {
3764 
3765  // for each track pick closest hit
3766  eTofHitVec hitVec = it->second ;
3767  eTofHitVec bestMatchVec;
3768  eTofHitVec mergeCandVec;
3769  eTofHitVec ambigVec;
3770  std::map< Int_t, StructETofHit > bestMatchMap;
3771  std::map< Int_t, eTofHitVec > mergeCandMap;
3772  std::map< Int_t, eTofHitVec > mergeCandMap2;
3773  std::map< Int_t, eTofHitVec > ambigMap;
3774  std::vector<int> indVec;
3775 
3776  for(unsigned int l =0; l < it->second.size(); l++){
3777 
3778  double dr = it->second.at(l).deltaX*it->second.at(l).deltaX + it->second.at(l).deltaY*it->second.at(l).deltaY;
3779  double dr_best = 99999.0; // dy for SHs at 27
3780  unsigned int ind = 0;
3781  unsigned int ind_best = l;
3782  int trackId = it->second.at(l).trackId;
3783  int vcnt = 0;
3784 
3785  if(std::find(indVec.begin(), indVec.end(), trackId) != indVec.end()) continue;
3786 
3787  for(unsigned int n = 0; n < it->second.size(); n++){
3788 
3789  if(it->second.at(n).trackId != trackId) continue;
3790 
3791  // localY doesnt contain any ETOF information for sHits-> take nHit if possible
3792  dr = it->second.at(n).deltaX*it->second.at(n).deltaX + it->second.at(n).deltaY*it->second.at(n).deltaY;
3793  ind = n;
3794 
3795  if(dr < dr_best){
3796 
3797  if(vcnt){
3798  mergeCandVec.push_back(it->second.at(ind_best));
3799  }else{
3800  vcnt++;
3801  }
3802  dr_best = dr;
3803  ind_best = ind;
3804 
3805  }else{
3806 
3807  mergeCandVec.push_back(it->second.at(n));
3808  }
3809  }
3810 
3811  indVec.push_back(trackId);
3812  bestMatchMap[trackId] = it->second.at(ind_best);
3813  bestMatchVec.push_back(it->second.at(ind_best));
3814  }
3815 
3816 
3817  std::vector<int> indVecBMtrack;
3818  std::vector<int> indVecBMhit;
3819 
3820  for(unsigned int b =0; b < bestMatchVec.size() ; b++){
3821  indVecBMtrack.push_back(bestMatchVec.at(b).trackId);
3822  indVecBMhit.push_back(bestMatchVec.at(b).index2ETofHit);
3823  }
3824 
3825  std::vector<int> indVecUsedTrack;
3826  std::vector<int> indVecReplaceTrack;
3827  std::vector<int> indVecUsedHit;
3828  eTofHitVec MatchVecTemp = bestMatchVec;
3829  eTofHitVec finalbestMatchVec;
3830 
3831  while(MatchVecTemp.size() > 0){
3832 
3833  double dr = 0.0;
3834  double dr_best = 99999.0; // dy for SHs at 27
3835  unsigned int ind = 0;
3836  unsigned int ind_best = 0;
3837  for(unsigned int b =0; b < MatchVecTemp.size() ; b++){
3838 
3839  ind = b;
3840 
3841  dr = MatchVecTemp.at(b).deltaX * MatchVecTemp.at(b).deltaX + MatchVecTemp.at(b).deltaY * MatchVecTemp.at(b).deltaY;
3842  if(dr <= dr_best){
3843  dr_best = dr;
3844  ind_best = ind;
3845  }
3846  }
3847 
3848  finalbestMatchVec.push_back(MatchVecTemp.at(ind_best));
3849  indVecUsedTrack.push_back(MatchVecTemp.at(ind_best).trackId);
3850  indVecUsedHit.push_back(MatchVecTemp.at(ind_best).index2ETofHit);
3851  MatchVecTemp.erase(MatchVecTemp.begin() + ind_best);
3852 
3853  //remove all matches with same hit id
3854  for(unsigned int b =0; b < MatchVecTemp.size() ; b++){
3855 
3856  if(std::find(indVecUsedHit.begin(), indVecUsedHit.end(), MatchVecTemp.at(b).index2ETofHit) != indVecUsedHit.end()) {
3857 
3858  indVecReplaceTrack.push_back(MatchVecTemp.at(b).trackId);
3859  MatchVecTemp.erase(MatchVecTemp.begin() + b);
3860  b = -1;
3861  }
3862  }
3863 
3864  //check for replacement
3865  std::sort( indVecReplaceTrack.begin(), indVecReplaceTrack.end() );
3866  indVecReplaceTrack.erase( unique( indVecReplaceTrack.begin(), indVecReplaceTrack.end() ), indVecReplaceTrack.end() );
3867 
3868  bool found1 = false;
3869  double dx1 = 0;
3870  double dy1 = 0;
3871  double dx_best1 = 99999.0;
3872  double dy_best1 = 99999.0;
3873  unsigned int ind1 = 0;
3874  unsigned int ind_best1 = 0;
3875  for(unsigned int i = 0; i < mergeCandVec.size();i++){
3876 
3877  ind1 = i;
3878 
3879  if(!(std::find(indVecReplaceTrack.begin(), indVecReplaceTrack.end(), mergeCandVec.at(i).trackId) != indVecReplaceTrack.end())) continue;
3880  if(std::find(indVecUsedTrack.begin(), indVecUsedTrack.end(), mergeCandVec.at(i).index2ETofHit) != indVecUsedTrack.end()) continue;
3881  if(std::find(indVecUsedHit.begin(), indVecUsedHit.end(), mergeCandVec.at(i).index2ETofHit) != indVecUsedHit.end()) continue;
3882  if(std::find(indVecBMhit.begin(), indVecBMhit.end(), mergeCandVec.at(i).index2ETofHit) != indVecBMhit.end()) continue;
3883 
3884  dx1 = mergeCandVec.at(i).deltaX;
3885  dy1 = mergeCandVec.at(i).deltaY;
3886 
3887  if(dy1 < dy_best1){ dy_best1 = dy1;}
3888 
3889  if(dx1 < dx_best1){
3890  dx_best1 = dx1;
3891  ind_best1 = ind1;
3892  found1 = true;
3893  } else if(dx1 == dx_best1){
3894 
3895  if(dy1 < dy_best1 && dy1 < dy_max && dy1 != 27.0){
3896  ind_best1 = ind1;
3897  found1 = true;
3898  } else if( dy1 == 27.0){
3899  ind_best1 = ind1;
3900  found1 = true;
3901  }
3902  }
3903  }
3904 
3905  if(found1){
3906 
3907  finalbestMatchVec.push_back(mergeCandVec.at(ind_best1));
3908  indVecUsedTrack.push_back(mergeCandVec.at(ind_best1).trackId);
3909  indVecUsedHit.push_back(mergeCandVec.at(ind_best1).index2ETofHit);
3910  mergeCandVec.erase(mergeCandVec.begin() + ind_best1);
3911  }
3912 
3913  bestMatchVec = finalbestMatchVec;
3914 
3915  for(unsigned int i=0;i< bestMatchVec.size();i++){
3916 
3917  if(bestMatchVec.at(i).clusterSize < 999 ){
3918  bestMatchVec.at(i).matchFlag = 410;
3919  }else{
3920  bestMatchVec.at(i).matchFlag = 420;
3921  bestMatchVec.at(i).clusterSize -= 1000;
3922  }
3923  finalMatchVec.push_back(bestMatchVec.at(i));
3924  }
3925  }
3926  }
3927  }// loop over MMMap
3928  }//loop over counters
3929 }
StMuDst::primaryVertex
static StMuPrimaryVertex * primaryVertex()
return pointer to current primary vertex
Definition: StMuDst.h:322
StMuETofHit::localY
double localY() const
Y-position.
Definition: StMuETofHit.h:203
StMuETofDigi::associatedHitId
int associatedHitId() const
Definition: StMuETofDigi.h:222
StMuETofHit::zPlane
unsigned int zPlane() const
ZPlane.
Definition: StMuETofHit.h:194
StMuDst::globalTracks
static TObjArray * globalTracks()
returns pointer to the global tracks list
Definition: StMuDst.h:303
StMuETofDigi::sector
unsigned int sector() const
Sector.
Definition: StMuETofDigi.h:212
StETofPidTraits::setETofHit
void setETofHit(StETofHit *hit)
PID functions – to be added (?)
Definition: StETofPidTraits.cxx:70
StMuTrack::vertexIndex
Int_t vertexIndex() const
Returns index of associated primary vertex.
Definition: StMuTrack.cxx:600
StMuETofHit::clusterSize
unsigned int clusterSize() const
Cluster size.
Definition: StMuETofHit.h:200
StMuETofHit::time
double time() const
Time.
Definition: StMuETofHit.h:197
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
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
StMuTrack::nHitsDedx
UShort_t nHitsDedx() const
Return number of hits used for dEdx.
Definition: StMuTrack.h:238
StMuETofHit::counter
unsigned int counter() const
Counter.
Definition: StMuETofHit.h:195
StETofGeometry::helixCrossCounter
void helixCrossCounter(const StPhysicalHelixD &helix, std::vector< int > &idVec, std::vector< StThreeVectorD > &crossVec, std::vector< StThreeVectorD > &localVec, std::vector< double > &thetaVec, std::vector< double > &pathLenVec)
Definition: StETofGeometry.cxx:1080
StMuETofHit::sector
unsigned int sector() const
Sector.
Definition: StMuETofHit.h:193
StETofHit::time
double time() const
Time.
Definition: StETofHit.h:193
StMuDst::etofHit
static StMuETofHit * etofHit(int i)
returns pointer to the i-th StMuETofHit
Definition: StMuDst.h:429
tof
Definition: tof.h:15
StMuTrack::etofPidTraits
const StMuETofPidTraits & etofPidTraits() const
dongx
Definition: StMuTrack.h:265
StETofHit::clusterSize
unsigned int clusterSize() const
Cluster size.
Definition: StETofHit.h:196
ETofTrack
ETOF track class.
Definition: StETofMatchMaker.h:58
StMuETofDigi::zPlane
unsigned int zPlane() const
ZPlane.
Definition: StMuETofDigi.h:213
StHelix::pathLength
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
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
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
StMuTrack::index2ETofHit
Int_t index2ETofHit() const
dongx
Definition: StMuTrack.h:235
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
StETofHit::counter
unsigned int counter() const
Counter.
Definition: StETofHit.h:191
StMuDst::etofArray
static TClonesArray * etofArray(int type)
returns pointer to the n-th TClonesArray from the etof arrays // FS
Definition: StMuDst.h:289
StMuETofHit::localX
double localX() const
X-position.
Definition: StMuETofHit.h:202
StMuTrack::nHitsPoss
UShort_t nHitsPoss() const
Return number of possible hits on track.
Definition: StMuTrack.cxx:242
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
StETofHit::associatedTrack
StTrack * associatedTrack()
pointer to the track which has been matched to this hit
Definition: StETofHit.h:201
StETofHit::zPlane
unsigned int zPlane() const
ZPlane.
Definition: StETofHit.h:190
StMuTrack::momentum
const StThreeVectorF & momentum() const
Returns 3-momentum at dca to primary vertex.
Definition: StMuTrack.h:260
StMuDst::etofDigi
static StMuETofDigi * etofDigi(int i)
returns pointer to the i-th StMuEtofDigi
Definition: StMuDst.h:427
StETofPidTraits::timeOfFlight
float timeOfFlight() const
timing for PID
Definition: StETofPidTraits.h:99
StETofHit::localX
double localX() const
X-position.
Definition: StETofHit.h:198
StMuETofHit::totalTot
double totalTot() const
Total Tot.
Definition: StMuETofHit.h:198
StMuTrack::dEdx
Double_t dEdx() const
Returns measured dE/dx value.
Definition: StMuTrack.h:248
StETofHit::totalTot
double totalTot() const
Total Tot.
Definition: StETofHit.h:194
StETofHit::localY
double localY() const
Y-position.
Definition: StETofHit.h:199
StMuDst::btofHeader
static StBTofHeader * btofHeader()
returns pointer to the btofHeader - dongx
Definition: StMuDst.h:423
StMuDst::array
static TClonesArray * array(int type)
returns pointer to the n-th TClonesArray
Definition: StMuDst.h:262
StMuTrack::primaryTrack
const StMuTrack * primaryTrack() const
Returns pointer to associated primary track. Null pointer if no global track available.
Definition: StMuTrack.cxx:578
StETofHit::sector
unsigned int sector() const
Sector.
Definition: StETofHit.h:189
StMuTrack::nSigmaProton
Double_t nSigmaProton() const
Returns Craig&#39;s distance to the calculated dE/dx band for protons in units of sigma.
Definition: StMuTrack.h:247
StETofPidTraits::matchFlag
unsigned short matchFlag() const
matching information
Definition: StETofPidTraits.h:93
StMuETofPidTraits::timeOfFlight
float timeOfFlight() const
timing for PID
Definition: StMuETofPidTraits.h:95
StMuTrack::nSigmaKaon
Double_t nSigmaKaon() const
Returns Craig&#39;s distance to the calculated dE/dx band for kaons in units of sigma.
Definition: StMuTrack.h:246
StMuTrack::setETofPidTraits
void setETofPidTraits(const StMuETofPidTraits &pid)
dongx
Definition: StMuTrack.h:270
StMuTrack::outerHelix
StPhysicalHelixD outerHelix() const
Returns outer helix (last measured point)
Definition: StMuTrack.cxx:412
StMuETofDigi::calibTot
double calibTot() const
Getter for calibrated Tot.
Definition: StMuETofDigi.h:210
StMuETofDigi::counter
unsigned int counter() const
Counter.
Definition: StMuETofDigi.h:214
kStOk
Definition: Stypes.h:41
TDataSet::Find
virtual TDataSet * Find(const char *path) const
Definition: TDataSet.cxx:362
StMuTrack::setIndex2ETofHit
void setIndex2ETofHit(Int_t i)
dongx
Definition: StMuTrack.h:146
StMuETofPidTraits::matchFlag
unsigned short matchFlag() const
matching information
Definition: StMuETofPidTraits.h:89
Generated by   doxygen 1.8.5