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