StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StETofHitMaker.cxx
1 /***************************************************************************
2  *
3  * $Id: StETofHitMaker.cxx,v 1.9 2021/01/30 19:19:16 weidenkaff Exp $
4  *
5  * Author: Philipp Weidenkaff & Florian Seck, April 2018
6  ***************************************************************************
7  *
8  * Description: StETofHitMaker - class to read the eTofCollection from
9  * StEvent and combine digis on both sides of each read-out strip into a hit.
10  * The hits on each strip are further merger into clusters.
11  * The eTOF collection is filled with hits and written to StEvent.
12  *
13  ***************************************************************************
14  *
15  * $Log: StETofHitMaker.cxx,v $
16  * Revision 1.9 2021/01/30 19:19:16 weidenkaff
17  * fixed deleting of created hits with StEvent data
18  *
19  * Revision 1.8 2021/01/29 15:08:31 weidenkaff
20  * fixed memory leak in StEtofHitMaker.cxx by adding a delete to merged hits
21  *
22  * Revision 1.7 2020/01/16 03:40:23 fseck
23  * add possibility to calculate VPD start time + updated deadtime handling for negative hit times
24  *
25  * Revision 1.6 2019/12/17 03:27:51 fseck
26  * update to histograms for .hist.root files
27  *
28  * Revision 1.5 2019/12/12 02:27:09 fseck
29  * enable clock jump correction by default
30  *
31  * Revision 1.4 2019/12/10 15:58:33 fseck
32  * ignore digis in dead time software-wise + possibility to correct clock jumps based on hit position via setting a flag
33  *
34  * Revision 1.3 2019/03/25 01:07:40 fseck
35  * added more correlation & average hit time histograms for offline QA
36  *
37  * Revision 1.2 2019/03/08 19:07:20 fseck
38  * introduced dead time handling + fixed clustering to only pick up hits on adjacent strips + moved QA histograms for clustered hits into separate function + added correlation plots to bTOF hits
39  *
40  * Revision 1.1 2019/02/19 19:52:28 jeromel
41  * Reviewed code provided by F.Seck
42  *
43  *
44  ***************************************************************************/
45 #include <algorithm>
46 #include <vector>
47 #include <cmath>
48 
49 #include "TFile.h"
50 #include "TH1F.h"
51 #include "TH2F.h"
52 #include "TH3F.h"
53 
54 #include "StEvent.h"
55 #include "StETofCollection.h"
56 #include "StETofDigi.h"
57 #include "StETofHit.h"
58 
59 #include "StBTofCollection.h"
60 #include "StBTofHeader.h"
61 #include "StBTofHit.h"
62 
63 #include "StEpdCollection.h"
64 #include "StEpdHit.h"
65 
66 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
67 #include "StMuDSTMaker/COMMON/StMuDst.h"
68 #include "StMuDSTMaker/COMMON/StMuETofDigi.h"
69 #include "StMuDSTMaker/COMMON/StMuETofHit.h"
70 #include "StMuDSTMaker/COMMON/StMuBTofHit.h"
71 #include "StMuDSTMaker/COMMON/StMuEpdHit.h"
72 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
73 
74 #include "StChain/StChainOpt.h" // for renaming the histogram file
75 
76 #include "StETofHitMaker.h"
77 #include "StETofUtil/StETofConstants.h"
78 #include "StETofUtil/StETofGeometry.h"
79 
80 #include "tables/St_etofHitParam_Table.h"
81 #include "tables/St_etofSignalVelocity_Table.h"
82 #include "tables/St_etofModCounter_Table.h"
83 
84 #include "StarClassLibrary/SystemOfUnits.h"
85 #include "StarClassLibrary/PhysicalConstants.h"
86 #include "StBTofUtil/tofPathLength.hh"
87 
88 
89 //_____________________________________________________________
90 StETofHitMaker::StETofHitMaker( const char* name )
91 : StMaker( "etofHit", name ),
92  mEvent( nullptr ),
93  mMuDst( nullptr ),
94  mETofGeom( nullptr ),
95  mFileNameHitParam( "" ),
96  mFileNameSignalVelocity( "" ),
97  mFileNameModMatrix( "" ),
98  mFileNameAlignParam( "" ),
99  mStoreDigi(),
100  mMapDigiIndex(),
101  mStoreHit(),
102  mMapHitDigiIndices(),
103  mMapHitIndexDigiIndices(),
104  mMaxYPos( 15. ),
105  mMergingRadius( 1. ),
106  mSigVel(),
107  mSoftwareDeadTime( 150. ),
108  mDoClockJumpShift( true ),
109  mDoDoubleClockJumpShift( true ),
110  mClockJumpDirection(),
111  mModMatrix(),
112  mGet4doublejumpTmin(-1.0),
113  mGet4doublejumpFlag(),
114  mGet4doublejumpTimes(),
115  mIsSim( false ),
116  mDoQA( false ),
117  mDebug( false ),
118  mApCorr(false),
119  mHistFileName( "" ),
120  mHistograms(),
121  mCounterActive()
122 {
123  LOG_DEBUG << "StETofHitMaker::ctor" << endm;
124 }
125 
126 //_____________________________________________________________
127 StETofHitMaker::~StETofHitMaker()
128 {
129  /* no op */
130 }
131 
132 //_____________________________________________________________
133 Int_t
134 StETofHitMaker::Init()
135 {
136  LOG_INFO << "StETofHitMaker::Init()" << endm;
137 
138  bookHistograms();
139 
140  return kStOk;
141 }
142 
143 //_____________________________________________________________
144 Int_t
145 StETofHitMaker::InitRun( Int_t runnumber )
146 {
147  LOG_INFO << "StETofHitMaker::InitRun()" << endm;
148 
149  TDataSet* dbDataSet = nullptr;
150  std::ifstream paramFile;
151 
152  // --------------------------------------------------------------------------------------------
153  // initialize hit building parameters from parameter file (if filename is provided) or database:
154  // -- hit param
155  // -- signal velocity
156  // --------------------------------------------------------------------------------------------
157 
158  // hit param
159  if( mFileNameHitParam.empty() ) {
160  LOG_INFO << "etofHitParam: no filename provided --> load database table" << endm;
161 
162  dbDataSet = GetDataBase( "Calibrations/etof/etofHitParam" );
163 
164  St_etofHitParam* etofHitParam = static_cast< St_etofHitParam* > ( dbDataSet->Find( "etofHitParam" ) );
165  if( !etofHitParam ) {
166  LOG_ERROR << "unable to get the hit params from the database" << endm;
167  return kStFatal;
168  }
169 
170  etofHitParam_st* hitParamTable = etofHitParam->GetTable();
171 
172  mMaxYPos = hitParamTable->maxLocalY;
173  mMergingRadius = hitParamTable->clusterMergeRadius;
174  }
175  else {
176  LOG_INFO << "etofHitParam: filename provided --> use parameter file: " << mFileNameHitParam.c_str() << endm;
177 
178  paramFile.open( mFileNameHitParam.c_str() );
179 
180  if( !paramFile.is_open() ) {
181  LOG_ERROR << "unable to get the 'etofHitParam' parameters from file --> file does not exist" << endm;
182  return kStFatal;
183  }
184 
185  float temp1 = 0;
186  float temp2 = 0;
187  if( paramFile.good() ) {
188  paramFile >> temp1 >> temp2;
189  }
190 
191  paramFile.close();
192 
193  if( temp1 > 0. ) {
194  mMaxYPos = temp1;
195  }
196  if( temp2 > 0. ) {
197  mMergingRadius = temp2;
198  }
199  }
200 
201  LOG_INFO << " maximum local Y: " << mMaxYPos << " , cluster merging radius: " << mMergingRadius << endm;
202 
203  // --------------------------------------------------------------------------------------------
204 
205  // signal velocities
206  mSigVel.clear();
207 
208  if( mFileNameSignalVelocity.empty() ) {
209  LOG_INFO << "etofSignalVelocity: no filename provided --> load database table" << endm;
210 
211  dbDataSet = GetDataBase( "Calibrations/etof/etofSignalVelocity" );
212 
213  St_etofSignalVelocity* etofSignalVelocity = static_cast< St_etofSignalVelocity* > ( dbDataSet->Find( "etofSignalVelocity" ) );
214 
215 
216  if( !etofSignalVelocity ) {
217  LOG_ERROR << "unable to get the signal velocity from the database" << endm;
218  return kStFatal;
219  }
220 
221  etofSignalVelocity_st* velocityTable = etofSignalVelocity->GetTable();
222 
223  for( size_t i=0; i<eTofConst::nCountersInSystem; i++ ) {
224  if( velocityTable->signalVelocity[ i ] > 0 ) {
225  mSigVel[ detectorToKey( i ) ] = velocityTable->signalVelocity[ i ];
226  }
227  }
228  }
229  else {
230  LOG_INFO << "etofSignalVelocity: filename provided --> use parameter file: " << mFileNameSignalVelocity.c_str() << endm;
231 
232  paramFile.open( mFileNameSignalVelocity.c_str() );
233 
234  if( !paramFile.is_open() ) {
235  LOG_ERROR << "unable to get the 'etofSignalVelocity' parameters from file --> file does not exist" << endm;
236  return kStFatal;
237  }
238 
239  std::vector< float > velocity;
240  float temp;
241  while( paramFile >> temp ) {
242  velocity.push_back( temp );
243  }
244 
245  paramFile.close();
246 
247  if( velocity.size() != eTofConst::nCountersInSystem ) {
248  LOG_ERROR << "parameter file for 'etofSignalVelocity' has not the right amount of entries: ";
249  LOG_ERROR << velocity.size() << " instead of " << eTofConst::nCountersInSystem << " !!!!" << endm;
250  return kStFatal;
251  }
252 
253  for( size_t i=0; i<eTofConst::nCountersInSystem; i++ ) {
254  if( velocity.at( i ) > 0 ) {
255  mSigVel[ detectorToKey( i ) ] = velocity.at( i );
256  }
257  }
258  }
259 
260  for( const auto& kv : mSigVel ) {
261  LOG_INFO << "counter key: " << kv.first << " --> signal velocity = " << kv.second << " cm / ns" << endm;
262  }
263 
264 
265  // --------------------------------------------------------------------------------------------
266  //initialize Counter Modification map (flip local xy etc.)
267 
268  if(!mFileNameModMatrix.empty()){ // "from file"
269 
270  LOG_INFO << "etofModMatrix: filename provided --> use parameter file: " << mFileNameSignalVelocity.c_str() << endm;
271 
272  paramFile.open( mFileNameModMatrix.c_str() );
273 
274  if( !paramFile.is_open() ) {
275  LOG_ERROR << "unable to get the 'etofModMatrix' parameters from file --> file does not exist" << endm;
276  return kStFatal;
277  }
278 
279  std::vector< int > modMatrix;
280  int temp;
281  while( paramFile >> temp ) {
282  modMatrix.push_back( temp );
283  }
284 
285  paramFile.close();
286 
287 
288  if( modMatrix.size() != eTofConst::nCountersInSystem ) {
289  LOG_ERROR << "parameter file for 'etofModMatrix' has not the right amount of entries: ";
290  LOG_ERROR << modMatrix.size() << " instead of " << eTofConst::nCountersInSystem << " !!!!" << endm;
291  return kStFatal;
292  }
293 
294  for( size_t i=0; i<eTofConst::nCountersInSystem; i++ ) {
295  if( modMatrix.at( i ) > 0 ) {
296  mModMatrix[ detectorToKey( i ) ] = modMatrix.at( i );
297  }
298  }
299 
300  }else{ //from Db
301 
302  dbDataSet = GetDataBase( "Geometry/etof/etofModCounter" );
303 
304  // TDataSet* dbDataSet = StMaker::GetChain()->GetDataBase("Geometry/etof/etofModCounter");
305  if( !dbDataSet ) {
306  LOG_ERROR << "unable to get the dataset from the database" << endm;
307  return kStFatal;
308  }
309 
310  St_etofModCounter* etofModCounter = static_cast< St_etofModCounter* > ( dbDataSet->Find("etofModCounter") );
311 
312  if( !etofModCounter ) {
313  LOG_WARN << "unable to get the ModMap from the database" << endm;
314  return kStFatal;
315 
316  }else{
317 
318  etofModCounter_st* ModCounterTable = etofModCounter->GetTable();
319 
320  for( size_t i=0; i<eTofConst::nCountersInSystem; i++ ) {
321  mModMatrix[ detectorToKey( i ) ] = ModCounterTable->detectorModFlag[ i ];
322  }
323  }
324  }
325 
326  // --------------------------------------------------------------------------------------------
327  for( int i=0; i<eTofConst::nCountersInSystem * 8; i++ ) {
328  int key = detectorToKey( i / 8 ) * 10 + ( i % 8 ) + 1;
329  mClockJumpDirection[ key ] = -1; //changed default to -1, read: backwards in time, to reduce losses before algoritm identifies direction.
330 
331  LOG_DEBUG << key << " " << mClockJumpDirection.at( key ) << endm;
332 
333  mGet4doublejumpFlag[ key ] = 0;
334 
335  for(int k=0; k<3; k++){
336  mGet4doublejumpTimes[key].push_back(-999);
337  }
338  }
339 
340  // --------------------------------------------------------------------------------------------
341  for(int i=0; i<eTofConst::nCountersInSystem; i++ ) {
342  mCounterActive.push_back( false );
343  }
344  // --------------------------------------------------------------------------------------------
345  // initializie etof geometry
346  // --------------------------------------------------------------------------------------------
347 
348  if( !mETofGeom ) {
349  LOG_INFO << " creating a new eTOF geometry . . . " << endm;
350  mETofGeom = new StETofGeometry( "etofGeometry", "etofGeometry in HitMaker" );
351  }
352 
353  if( mETofGeom && !mETofGeom->isInitDone() ) {
354  LOG_INFO << " eTOF geometry initialization ... " << endm;
355 
356  if( !gGeoManager ) GetDataBase( "VmcGeometry" );
357 
358  if( !gGeoManager ) {
359  LOG_ERROR << "Cannot get GeoManager" << endm;
360  return kStFatal;
361  }
362 
363  LOG_DEBUG << " gGeoManager: Should set alignment file now! " << mFileNameAlignParam <<" ! "<< endm;
364  if (mFileNameAlignParam != ""){
365  LOG_INFO << " gGeoManager: Setting alignment file: " << mFileNameAlignParam << endm;
366  mETofGeom->setFileNameAlignParam(mFileNameAlignParam);
367  }
368  mETofGeom->init(gGeoManager);
369  LOG_DEBUG << " init done " << endm;
370  }
371 
372  return kStOk;
373 }
374 
375 //_____________________________________________________________
376 Int_t
377 StETofHitMaker::FinishRun( Int_t runnumber )
378 {
379  LOG_INFO << "StETofHitMaker::FinishRun()" << endm;
380  if( mDoQA ) {
381  for( int iCounter = 0; iCounter < eTofConst::nCountersInSystem; iCounter++ ) {
382  if( mCounterActive.at( iCounter ) ) {
383  int day = ( runnumber % 1000000 ) / 1000;
384  int run = runnumber % 1000;
385 
386  mHistograms.at( "counter_active" )->Fill( 100 * day + run, iCounter );
387  }
388  }
389  }
390 
391  if( mETofGeom ) {
392  mETofGeom->reset();
393  }
394 
395  return kStOk;
396 }
397 
398 //_____________________________________________________________
399 Int_t
400 StETofHitMaker::Finish()
401 {
402  LOG_INFO << "StETofHitMaker::Finish()" << endm;
403 
404  if( mDoQA ) {
405  LOG_INFO << "Finish() - writing *.etofHit.root ..." << endm;
406  setHistFileName();
407  writeHistograms();
408  }
409 
410  return kStOk;
411 }
412 
413 //_____________________________________________________________
414 Int_t
415 StETofHitMaker::Make()
416 {
417  LOG_DEBUG << "StETofHitMaker::Make(): starting ..." << endm;
418 
419  mEvent = ( StEvent* ) GetInputDS( "StEvent" );
420  // mEvent = NULL; //don't check for StEvent for genDst.C testing. PW
421 
422  if ( mEvent ) {
423  LOG_DEBUG << "Make() - running on StEvent" << endm;
424  StETofCollection* etofCollection = mEvent->etofCollection();
425 
426  if( !etofCollection ) { //additional check for empty StEvents structures produced by other Makers. Needed for genDst.C
427  LOG_WARN << "Make() - Found StEvent data structure, but no eTOF collection. Try MuDst processing instead" << endm;
428  mMuDst = ( StMuDst* ) GetInputDS( "MuDst" );
429 
430  if( mMuDst ) {
431  LOG_DEBUG << "Make() - running on MuDsts" << endm;
432 
433  processMuDst();
434 
435  return kStOk;
436  }
437  }
438 
439  processStEvent();
440 
441  return kStOk;
442  }
443  else {
444  LOG_DEBUG << "Make(): no StEvent found" << endm;
445 
446  mMuDst = ( StMuDst* ) GetInputDS( "MuDst" );
447 
448  if( mMuDst ) {
449  LOG_DEBUG << "Make() - running on MuDsts" << endm;
450 
451  processMuDst();
452 
453  return kStOk;
454  }
455  else {
456  LOG_WARN << "Make() - no StMuDst or StEvent" << endm;
457  return kStOk;
458  }
459  }
460 }
461 
462 //_____________________________________________________________
463 void
464 StETofHitMaker::processStEvent()
465 {
466  StETofCollection* etofCollection = mEvent->etofCollection();
467 
468  if( !etofCollection ) {
469  LOG_WARN << "processStEvent() - no etof collection" << endm;
470  return;
471  }
472 
473  if( !etofCollection->digisPresent() ) {
474  LOG_WARN << "processStEvent() - no digis present" << endm;
475  return;
476  }
477 
478  const StSPtrVecETofDigi& etofDigis = etofCollection->etofDigis();
479 
480  //---------------------------------
481 
482  size_t nDigis = etofDigis.size();
483  if( mDebug ) {
484  LOG_DEBUG << "processStEvent() - # fired eTOF digis : " << nDigis << endm;
485  }
486 
487  bool isMuDst = false;
488 
489  // clear existing hits from eTofCollection (important for afterburner mode )
490  clearHits( isMuDst );
491 
492  // clear all storage containers
493  //(to be sure that there are no digis left from the previous event)
494  clearStorage();
495 
496  //sort digis into the storage (one vector of digis for each strip)
497  size_t nDigisInStore = 0;
498 
499  for( size_t i = 0; i<nDigis; i++ ) {
500  StETofDigi* aDigi = etofDigis[ i ];
501 
502  if( fillStorage( aDigi, i ) ) {
503  nDigisInStore++;
504  }
505  }
506  // LOG_INFO << "processStEvent() - storage is filled with " << nDigisInStore << " digis" << endm;
507 
508  matchSides();
509 
510  double tstart = startTime();
511 
512  if( mDoQA ) {
513  fillUnclusteredHitQA( tstart, isMuDst );
514  }
515 
516  mergeClusters( isMuDst );
517 
518  assignAssociatedHits( isMuDst );
519 
520 
521  if( mDoQA ) {
522  mHistograms.at( "multiplicity_etofDigis_etofHits" )->Fill( nDigisInStore, etofCollection->etofHits().size() );
523  }
524 
525  if( etofCollection->hitsPresent() ) {
526  fillHitQA( isMuDst, tstart );
527  }
528 }
529 
530 //_____________________________________________________________
531 void
532 StETofHitMaker::processMuDst()
533 {
534  if( !mMuDst->etofArray( muETofDigi ) ) {
535  LOG_WARN << "processMuDst() - no digi array" << endm;
536  return;
537  }
538 
539  if( !mMuDst->numberOfETofDigi() ) {
540  LOG_WARN << "processMuDst() - no digis present" << endm;
541  return;
542  }
543 
544  //---------------------------------
545 
546  size_t nDigis = mMuDst->numberOfETofDigi();
547  if( mDebug ) {
548  LOG_DEBUG << "processMuDst() - # fired eTOF digis : " << nDigis << endm;
549  }
550  bool isMuDst = true;
551 
552  // clear existing hits from eTofCollection (important for afterburner mode )
553  clearHits( isMuDst );
554 
555  // clear all storage containers
556  //(to be sure that there are no digis left from the previous event)
557  clearStorage();
558 
559  //sort digis into the storage (one vector of digis for each strip)
560  size_t nDigisInStore = 0;
561 
562  for( size_t i = 0; i<nDigis; i++ ) {
563  StMuETofDigi* aDigi = mMuDst->etofDigi( i );
564 
565  if( fillStorage( aDigi, i ) ) {
566  nDigisInStore++;
567  }
568  }
569  // LOG_INFO << "processMuDst() - storage is filled with " << nDigisInStore << " digis" << endm;
570 
571  matchSides();
572 
573  double tstart = startTime();
574 
575  if( mDoQA ) {
576  fillUnclusteredHitQA( tstart, isMuDst );
577  }
578 
579  mergeClusters( isMuDst );
580 
581  assignAssociatedHits( isMuDst );
582 
583 
584  if( mDoQA ) {
585  mHistograms.at( "multiplicity_etofDigis_etofHits" )->Fill( nDigisInStore, mMuDst->numberOfETofHit() );
586  }
587 
588 
589  if( mMuDst->numberOfETofHit() ) {
590  fillHitQA( isMuDst, tstart );
591  }
592 }
593 //_____________________________________________________________
594 
595 
596 //_____________________________________________________________
597 // get the start time -- from bTOF header for now
598 double
599 StETofHitMaker::startTime()
600 {
601  if( mDebug ) {
602  LOG_INFO << "startTime(): -- loading start time from bTOF header" << endm;
603  }
604 
605  if(mIsSim){
606 
607  LOG_INFO << "mIsSim = true --> startTime set to 0" << endm;
608 
609  return 0.;
610  }
611 
612  StBTofHeader* btofHeader = nullptr;
613 
614  if( mEvent ) {
615  StBTofCollection* btofCollection = ( StBTofCollection* ) mEvent->btofCollection();
616 
617  if ( btofCollection ) {
618  btofHeader = btofCollection->tofHeader();
619  }
620  else {
621  LOG_DEBUG << "no StBTofCollection found by getTstart" << endm;
622  return -9999.;
623  }
624  }
625  else if( mMuDst ) {
626  btofHeader = mMuDst->btofHeader();
627  }
628 
629  if( !btofHeader ) {
630  LOG_DEBUG << "startTime(): -- no bTOF header --> no start time avaiable" << endm;
631  return -9999.;
632  }
633 
634  double tstart = btofHeader->tStart();
635 
636  if( !isfinite( tstart ) ) {
637  LOG_DEBUG << "startTime(): -- from bTOF header is NaN" << endm;
638  return -9999.;
639  }
640 
641  if( tstart != -9999. ) {
642  tstart = fmod( tstart, eTofConst::bTofClockCycle );
643  if( tstart < 0. ) tstart += eTofConst::bTofClockCycle;
644  }
645 
646  if( mDebug ) {
647  LOG_INFO << "startTime(): -- start time: " << tstart << endm;
648  }
649 
650  return tstart;
651 }
652 
653 //_____________________________________________________________
654 // get the start time (and vertex) -- from VPD
655 void
656 StETofHitMaker::startTimeVpd( double& tstart, double& vertexVz )
657 {
658  tstart = -9999.;
659  vertexVz = -999.;
660 
661  if( mDebug ) {
662  LOG_INFO << "startTimeVpd(): -- calculating VPD start time from bTOF header" << endm;
663  }
664 
665  StBTofHeader* btofHeader = nullptr;
666 
667  if( mEvent ) {
668  StBTofCollection* btofCollection = ( StBTofCollection* ) mEvent->btofCollection();
669 
670  if ( btofCollection ) {
671  btofHeader = btofCollection->tofHeader();
672  }
673  else {
674  LOG_DEBUG << "no StBTofCollection found by getTstart" << endm;
675  return;
676  }
677  }
678  else if( mMuDst ) {
679  btofHeader = mMuDst->btofHeader();
680  }
681 
682  if( !btofHeader ) {
683  LOG_DEBUG << "startTimeVpd(): -- no bTOF header --> no start time avaiable" << endm;
684  return;
685  }
686 
687  const int nVpd = 19; // number of VPD tubes on each side of STAR
688 
689  int nWest = btofHeader->numberOfVpdHits( west );
690  int nEast = btofHeader->numberOfVpdHits( east );
691 
692  double vpdLeTime[ 2 * nVpd ];
693 
694 
695  // check if bTof header is filled with useful information
696  if( fabs( btofHeader->vpdVz() ) > 200. ) {
697  if( mDoQA ) {
698  LOG_INFO << "startTimeVpd(): no valid Vpd data in the bTOF header " << endm;
699  }
700  return;
701  }
702  else {
703  vertexVz = btofHeader->vpdVz();
704  LOG_DEBUG << "startTimeVpd(): Vpd vertex is at: " << vertexVz << endm;
705  }
706 
707  double tMean = 0.;
708  int nTubes = 0;
709  int nTubesWest = 0;
710  int nTubesEast = 0;
711 
712  // west side
713  for( int i=0; i< nVpd; i++ ) {
714  vpdLeTime[ i ] = btofHeader->vpdTime( west, i+1 );
715  if( vpdLeTime[ i ] > 0. ) {
716  updateCyclicRunningMean( vpdLeTime[ i ], tMean, nTubes, eTofConst::bTofClockCycle );
717  nTubesWest++;
718  LOG_DEBUG << "startTimeVpd(): loading VPD west tubeId = " << i+1 << " time " << vpdLeTime[ i ] << endm;
719  }
720  }
721 
722  // east side
723  for( int i=0; i< nVpd; i++ ) {
724  vpdLeTime[ i + nVpd ] = btofHeader->vpdTime( east, i+1 );
725  if( vpdLeTime[ i + nVpd ] > 0. ) {
726  updateCyclicRunningMean( vpdLeTime[ i + nVpd ], tMean, nTubes, eTofConst::bTofClockCycle );
727  nTubesEast++;
728  LOG_DEBUG << "startTimeVpd(): loading VPD east tubeId = " << i+1 << " time " << vpdLeTime[ i + nVpd ] << endm;
729  }
730  }
731 
732  if( nTubesEast >= 2 && nTubesWest >= 2 ) {
733  tstart = tMean;
734  }
735 
736  if( tstart != -9999 ) {
737  tstart = fmod( tstart, eTofConst::bTofClockCycle );
738  if( tstart < 0. ) tstart += eTofConst::bTofClockCycle;
739  }
740 
741  if( mDoQA ) {
742  LOG_INFO << "startTimeVpd(): -- sum: " << tMean << " nWest: " << nWest << " nEast: " << nEast;
743  LOG_INFO << " --> compare (if bTofHeader is filled): " << btofHeader->tStart() << endm;
744  LOG_INFO << "startTimeVpd(): -- start time: " << tstart << endm;
745  }
746 }
747 
748 
749 //_____________________________________________________________
750 bool
751 StETofHitMaker::fillStorage( StETofDigi* aDigi, unsigned int index )
752 {
753  if( !aDigi ) {
754  LOG_WARN << "No digi found" << endm;
755  return false;
756  }
757 
758  if( fabs( aDigi->calibTime() ) < 1e-5 && aDigi->calibTot() < 0 ) {
759  if( mDebug ) {
760  LOG_DEBUG << "fillStorage() - digi not calibrated, most likely since it is outside the trigger window or pulser. Ignore." << endm;
761  }
762  return false;
763  }
764 
765  if( aDigi->sector() == 0 || aDigi->zPlane() == 0 || aDigi->counter() == 0 || aDigi->strip() == 0 ) {
766  LOG_WARN << "fillStorage() - sector / zPlane / counter / strip was not assigned to the digi" << endm;
767  return false;
768  }
769 
770  StETofDigi* pDigi = new StETofDigi( *aDigi );
771 
772  if( mDebug ) {
773  LOG_DEBUG << "fillStorage() -- sector plane counter strip: " << pDigi->sector() << " ";
774  LOG_DEBUG << pDigi->zPlane() << " " << pDigi->counter() << " " << pDigi->strip() << endm;
775  LOG_DEBUG << "fillStorage() -- calibTime: " << pDigi->calibTime();
776  LOG_DEBUG << " calibToT: " << pDigi->calibTot() << endm;
777  }
778 
779  unsigned int key = pDigi->sector() * 10000 +
780  pDigi->zPlane() * 1000 +
781  pDigi->counter() * 100 +
782  pDigi->strip();
783 
784  mStoreDigi[ key ].push_back( pDigi );
785 
786  mMapDigiIndex[ pDigi ] = index;
787 
788 
789  if( mDoQA ) {
790  std::string histName_digi_tot = "digi_tot_s" + std::to_string( pDigi->sector() ) + "m" + std::to_string( pDigi->zPlane() ) + "c" + std::to_string( pDigi->counter() );
791 
792  mHistograms[ histName_digi_tot ]->Fill( pDigi->strip() - 1 + pDigi->side() * 0.3, pDigi->calibTot() );
793  }
794 
795  return true;
796 }//::fillStorage
797 
798 //_____________________________________________________________
799 bool
800 StETofHitMaker::fillStorage( StMuETofDigi* aDigi, unsigned int index )
801 {
802  return fillStorage( ( StETofDigi* ) aDigi, index );
803 }//::fillStorage
804 
805 
806 //_____________________________________________________________
807 void
808 StETofHitMaker::clearHits( const bool isMuDst )
809 {
810  if( !isMuDst ) {
811  if( !( mEvent->etofCollection()->hitsPresent() ) ) {
812  return;
813  }
814 
815  LOG_DEBUG << "clearHits() - number of hits present (before clear): " << mEvent->etofCollection()->etofHits().size() << endm;
816 
817  // clear hits (if there are any when running in afterburner mode)
818  StSPtrVecETofHit& etofHits = mEvent->etofCollection()->etofHits();
819  etofHits.clear();
820 
821  LOG_DEBUG << "clearHits() - number of hits present (after clear): " << mEvent->etofCollection()->etofHits().size() << endm;
822 
823  // and remove pointers to associated hit of the digis
824  StSPtrVecETofDigi& etofDigis = mEvent->etofCollection()->etofDigis();
825  size_t nDigis = etofDigis.size();
826 
827  for( size_t i=0; i<nDigis; i++ ) {
828  StETofDigi* aDigi = etofDigis[ i ];
829 
830  if( !aDigi ) continue;
831 
832  aDigi->setAssociatedHit( nullptr );
833  }
834  LOG_DEBUG << "clearHits() - associated hits of digis set to nullptr" << endm;
835  }
836  else {
837  if( mMuDst->numberOfETofHit() == 0 ) {
838  return;
839  }
840 
841  LOG_DEBUG << "clearHits() - number of hits present (before clear): " << mMuDst->numberOfETofHit() << endm;
842 
843  // clear hits (if there are any when running in afterburner mode)
844  mMuDst->etofArray( muETofHit )->Clear( "C" );
845 
846  LOG_DEBUG << "clearHits() - number of hits present (after clear): " << mMuDst->numberOfETofHit() << endm;
847 
848  // and remove pointers to associated hit of the digis
849  size_t nDigis = mMuDst->numberOfETofDigi();
850 
851  for( size_t i=0; i<nDigis; i++ ) {
852  StMuETofDigi* aDigi = mMuDst->etofDigi( i );
853 
854  if( !aDigi ) continue;
855 
856  aDigi->setAssociatedHitId( -1 );
857  }
858  LOG_DEBUG << "clearHits() - associated hit id of digis set to -1" << endm;
859  }
860 }
861 
862 
863 //_____________________________________________________________
864 void
865 StETofHitMaker::clearStorage()
866 {
867  // clear mStoreDigi -- if any digis are left
868  for( auto kv = mStoreDigi.begin(); kv != mStoreDigi.end(); kv++ ) {
869  size_t remainingDigis = kv->second.size();
870 
871  if( mDebug ) {
872  LOG_DEBUG << "strip key " << kv->first << " has " << remainingDigis << " left" << endm;
873  }
874 
875  for( unsigned int i=0; i<remainingDigis; i++ ) {
876  delete kv->second.at( i );
877  }
878  kv->second.clear();
879  }
880  mStoreDigi.clear();
881 
882 
883  // clear mStoreHit -- if any hits are left
884  for( auto kv = mStoreHit.begin(); kv != mStoreHit.end(); kv++ ) {
885  size_t remainingHits = kv->second.size();
886 
887  if( mDebug ) {
888  LOG_DEBUG << "detector key " << kv->first << " has " << remainingHits << " left" << endm;
889  }
890 
891  for( size_t i=0; i<remainingHits; i++ ) {
892  delete kv->second.at( i );
893  }
894  kv->second.clear();
895  }
896  mStoreHit.clear();
897 
898 
899  // clear mMapDigiIndex
900  mMapDigiIndex.clear();
901 
902  // clear mMapHitDigiIndices
903  mMapHitDigiIndices.clear();
904 
905  // clear mMapHitIndexDigiIndices
906  mMapHitIndexDigiIndices.clear();
907 }//::clearStorage
908 
909 //_____________________________________________________________
910 void
911 StETofHitMaker::matchSides()
912 {
913  std::vector< StETofDigi* >::iterator iterDigi;
914 
915  std::string histNameDigisErased;
916 
917  for( auto kv = mStoreDigi.begin(); kv != mStoreDigi.end(); kv++ ) {
918  unsigned int stripIndex = kv->first;
919  std::vector< StETofDigi* > *digiVec = &( kv->second );
920 
921  // timeorder digis from both sides via lambda functions of C++11
922  std::sort( digiVec->begin(), digiVec->end(), [] ( StETofDigi* lhs, StETofDigi* rhs ) {
923  return lhs->calibTime() < rhs->calibTime();
924  }
925  );
926 
927  int nDigisOnStrip = digiVec->size();
928  //--------------------------------------------------------------------------------
929  // print out for testing
930  if( mDebug ) {
931  LOG_INFO << stripIndex << " size: " << nDigisOnStrip << endm;
932 
933  for( size_t i=0; i<digiVec->size(); i++ ) {
934  LOG_INFO << "matchSides() - DIGI: " << digiVec->at( i ) << " ";
935  LOG_INFO << "calibTime=" << setprecision( 16 ) << digiVec->at( i )->calibTime() << " " << endm;
936  LOG_INFO << "calibTot=" << setprecision( 4 ) << digiVec->at( i )->calibTot() << " ";
937  LOG_INFO << "side=" << setprecision( 4 ) << digiVec->at( i )->side() << endm;
938  }
939  }
940  //--------------------------------------------------------------------------------
941 
942  int detIndex = stripIndex / 100;
943  int sector = detIndex / 100;
944  int plane = ( detIndex % 100 ) / 10;
945  int counter = detIndex % 10;
946  int strip = stripIndex % 100;
947 
948  if( mDoQA ) {
949  std::string histNameDigisPerStrip = "digisPerStrip_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
950  mHistograms.at( histNameDigisPerStrip )->Fill( strip, nDigisOnStrip );
951 
952  histNameDigisErased = "digisErased_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
953  }
954 
955 
956  //--------------------------------------------------------------------------------
957  // remove digis from the vector that fall within the dead time caused by another digi on the same side of a strip
958  std::vector< double > deadTime( 2, -60. );
959 
960  for( auto it = digiVec->begin(); it != digiVec->end(); it++ ) {
961  if( (*it)->rawTime() - deadTime.at( (*it)->side() - 1 ) < mSoftwareDeadTime ) { //TODO: move to DB. use raw time
962 
963  if( mDebug ) {
964  LOG_INFO << "digi within dead time --> ignore ... ( geomId : " << stripIndex * 10 + (*it)->side();
965  LOG_INFO << " dead time: " << setprecision( 16 ) << deadTime.at( (*it)->side() - 1 );
966  LOG_INFO << " calib time: " << setprecision( 16 ) << (*it)->calibTime();
967  LOG_INFO << " difference: " << (*it)->calibTime() - deadTime.at( (*it)->side() - 1 ) << endm;
968  }
969 
970  delete *it;
971  digiVec->erase( it );
972  if( mDoQA ){
973  mHistograms[ histNameDigisErased ]->Fill(1);
974  }
975  it--;
976  }
977  else {
978  deadTime.at( (*it)->side() - 1 ) = (*it)->rawTime();
979  }
980  }
981  //--------------------------------------------------------------------------------
982 
983  std::vector< unsigned int > containedDigiIndices; //
984  double posX = 0.0;
985  double posY = 0.0;
986  double time = 0.0;
987  double timeDiff = 0.0;
988  double totSum = 0.0;
989  double t_corr_afterpulse = 0.0;
990 
991 
992  if( mDoQA && digiVec->size() == 1 ) {
993  mHistograms.at( histNameDigisErased )->Fill( 2 );
994  }
995 
996 
997  //single sided digi hit building
998  if( digiVec->size() == 1 ) {
999 
1000  // create the hit candidate:
1001  StETofDigi* xDigiA = digiVec->at( 0 );
1002  StETofDigi* xDigiB = digiVec->at( 0 );
1003 
1004  //get get4flag statistics
1005  // StMuETofHeader* etofHeader = mMuDst->etofHeader();
1006  // TClass* headerClass = etofHeader->IsA();
1007  // std::vector< Bool_t > vMissmatchVec = etofHeader->missMatchFlagVec();
1008  // std::vector< bool > goodEventFlagVec = mMuDst->etofHeader()->goodEventFlagVec();
1009 
1010  // the "strip" time is the mean time between each end
1011  time = 0.5 * ( xDigiA->calibTime() + xDigiB->calibTime() );
1012  //TODO: Afterpulse handling: correct hit time by the time difference between the first and second digi on the same side
1013  if(!mIsSim && mApCorr){//merge skip corrections for simulation
1014  time += t_corr_afterpulse;
1015  }//merge
1016  // weight of merging of hits (later) is the total charge => sum of both ends ToT
1017  totSum = xDigiA->calibTot() + xDigiB->calibTot();
1018 
1019  if(xDigiA->side() == 1){
1020  posY = 1;
1021  }else{
1022  posY = -1;
1023  }
1024 
1025 
1026  // use local coordinates... (0,0,0) is in the center of counter
1027  posX = ( -1 * eTofConst::nStrips / 2. + strip - 0.5 ) * eTofConst::stripPitch;
1028 
1029  unsigned int clusterSize = 1000;
1030 
1031  StETofHit* constructedHit = new StETofHit( sector, plane, counter, time, totSum, clusterSize, posX, posY );
1032 
1033  mStoreHit[ detIndex ].push_back( constructedHit );
1034 
1035  containedDigiIndices.push_back( mMapDigiIndex.at( xDigiA ) );
1036  containedDigiIndices.push_back( mMapDigiIndex.at( xDigiB ) );
1037 
1038  mMapHitDigiIndices[ constructedHit ] = containedDigiIndices;
1039 
1040  }
1041 
1042 
1043  // loop over digis on the same strip
1044  while( digiVec->size() > 1 ) {
1045  if( mDebug ) {
1046  LOG_DEBUG << stripIndex << " -- digiVec->size() -- " << digiVec->size() << endm;
1047  }
1048  // treat consecutive digis on the same side:
1049  // we want to have the first and second digi to be on different sides
1050  // of the strip in order to build hits out of them
1051  while( digiVec->at( 0 )->side() == digiVec->at( 1 )->side() ) {
1052 
1053  if( digiVec->size() > 2 ) { //more than 2 digis left on strip
1054 
1055  // test for three (or more) consecutive digis on the same side
1056  if( digiVec->at( 2 )->side() == digiVec->at( 0 )->side() ) {
1057 
1058  // delete first digi
1059  iterDigi = digiVec->begin();
1060  delete *iterDigi;
1061  digiVec->erase( iterDigi );
1062  if( mDoQA ) {
1063  mHistograms.at( histNameDigisErased )->Fill( 3 );
1064  }
1065  }
1066  else { // --> third digi is on the other side compared to first and second digi
1067  if( digiVec->at( 2 )->calibTime() - digiVec->at( 0 )->calibTime() >
1068  digiVec->at( 2 )->calibTime() - digiVec->at( 1 )->calibTime() ) {
1069 
1070  // third digi is not same side and fits better with second digi
1071  // --> delete first digi
1072  iterDigi = digiVec->begin();
1073  delete *iterDigi;
1074  digiVec->erase( iterDigi );
1075  //TODO: Afterpulse handling: save time difference between digi 1 and digi 2 and substract from hit time!
1076  if(mApCorr) t_corr_afterpulse = digiVec->at( 0 )->calibTime() - digiVec->at( 1 )->calibTime();
1077  if( mDoQA ) {
1078  mHistograms.at( histNameDigisErased )->Fill( 4 );
1079  }
1080  }
1081  else {
1082  // third digi is not same side and fits better with first digi
1083  // --> delete second digi
1084  iterDigi = digiVec->begin() + 1;
1085  delete *iterDigi;
1086  digiVec->erase( iterDigi );
1087  if( mDoQA ) {
1088  mHistograms.at( histNameDigisErased )->Fill( 7 ); //TODO: Seperate in QA. Missed afterpulse side?! Should be rare
1089  }
1090  }
1091  }
1092  }
1093  else{ // --> 2 or less digis left on the strip (on the same side of the strip)
1094  // delete the remaining digi
1095  iterDigi = digiVec->begin();
1096  delete *iterDigi;
1097  digiVec->erase( iterDigi );
1098  if( mDoQA && digiVec->size() == 1 ){
1099  mHistograms.at( histNameDigisErased )->Fill( 5 );
1100  }
1101  }
1102 
1103  if( digiVec->size() < 2 ) { //only one digi left on strip. break loop.
1104  if(mDoQA && digiVec->size() == 1) {
1105  mHistograms.at( histNameDigisErased )->Fill( 5 );
1106  }
1107  break;
1108  }
1109  } // first and second digi in the vector are on different sides
1110 
1111  if( mDebug ) {
1112  LOG_DEBUG << "matchSides() - digi processing for sector " << stripIndex / 10000;
1113  LOG_DEBUG << " plane " << ( stripIndex % 10000 ) / 1000 << " counter " << ( stripIndex % 1000 ) / 100;
1114  LOG_DEBUG << " strip " << stripIndex % 100;
1115  LOG_DEBUG << " size: " << digiVec->size() << endm;
1116  }
1117 
1118  if( digiVec->size() < 2 ) {
1119  // only one digi left on strip. break loop.
1120  if( mDoQA ) {
1121  mHistograms.at( histNameDigisErased )->Fill( 5 );
1122  }
1123  break;
1124  }
1125 
1126  // two digis --> both sides present
1127  StETofDigi* xDigiA = digiVec->at( 0 );
1128  StETofDigi* xDigiB = digiVec->at( 1 );
1129 
1130  timeDiff = xDigiA->calibTime() - xDigiB->calibTime();
1131 
1132  if( mDebug ) {
1133  LOG_DEBUG << "matchSides() - time difference in ns: " << timeDiff << endm;
1134  }
1135 
1136  // side 1 is the top, side 2 is bottom
1137  if( xDigiA->side() == 2 ) {
1138  posY = mSigVel.at( detIndex ) * timeDiff * 0.5;
1139  }
1140  else {
1141  posY = -1 * mSigVel.at( detIndex ) * timeDiff * 0.5;
1142  }
1143 
1144 
1145  // check for a better match if the local y position is outside the detector bounds
1146  if( fabs( posY ) > mMaxYPos && digiVec->size() > 2 ) {
1147  if( mDebug ) {
1148  LOG_DEBUG << "matchSides() - hit candidate outside correlation window, check for better possible digis" << endm;
1149  LOG_DEBUG << "size of digi vector: " << digiVec->size() << endm;
1150  }
1151 
1152  StETofDigi* xDigiC = digiVec->at( 2 );
1153 
1154  double posYNew = 0.;
1155  double timeDiffNew = 0.;
1156 
1157  if( xDigiC->side() == xDigiA->side() ) {
1158  timeDiffNew = xDigiC->calibTime() - xDigiB->calibTime();
1159  }
1160  else {
1161  timeDiffNew = xDigiA->calibTime() - xDigiC->calibTime();
1162  }
1163 
1164  if( xDigiA->side() == 2 ) {
1165  posYNew = mSigVel.at( detIndex ) * timeDiffNew * 0.5;
1166  }
1167  else {
1168  posYNew = -1 * mSigVel.at( detIndex ) * timeDiffNew * 0.5;
1169  }
1170 
1171  if( fabs( posYNew ) < fabs( posY ) ) {
1172  if( mDebug ) {
1173  LOG_DEBUG << "matchSides() - found better match for hit candidate -> changing out digis" << endm;
1174  }
1175 
1176  timeDiff = timeDiffNew;
1177  posY = posYNew;
1178 
1179  if( xDigiC->side() == xDigiA->side() ) {
1180  //TODO: Afterpulse handling: save time difference between digi 1 and digi 2 and substract from hit time!
1181  if(mApCorr) t_corr_afterpulse = xDigiA->calibTime() - xDigiC->calibTime();
1182  xDigiA = xDigiC;
1183  iterDigi = digiVec->begin();
1184  delete *iterDigi;
1185  digiVec->erase( iterDigi );
1186  if( mDoQA ){
1187  mHistograms.at( histNameDigisErased )->Fill( 4 );
1188  }
1189  }
1190  else {
1191  //TODO: Afterpulse handling: save time difference between digi 1 and digi 2 and substract from hit time!
1192  if(mApCorr) t_corr_afterpulse = xDigiB->calibTime() - xDigiC->calibTime();
1193  xDigiB = xDigiC;
1194  iterDigi = digiVec->begin() + 1;
1195  delete *iterDigi;
1196  digiVec->erase( iterDigi );
1197  if( mDoQA ){
1198  mHistograms.at( histNameDigisErased )->Fill( 4 );
1199  }
1200  }
1201  }
1202  else { // --> keeps candidate even if it is outside correlation window
1203  if( mDebug ) {
1204  LOG_DEBUG << "matchSides() - no better match -> keep this hit candidate" << endm;
1205  }
1206  }
1207  } // check for better match with third digi done
1208 
1209 
1210  if( xDigiA->side() == xDigiB->side() ) {
1211  LOG_ERROR << "matchSides() - wrong combinations of digis:" << endm;
1212  LOG_ERROR << *xDigiA << endm;
1213  LOG_ERROR << *xDigiB << endm;
1214  }
1215 
1216 
1217 
1218  // create the hit candidate:
1219  // the "strip" time is the mean time between each end
1220  time = 0.5 * ( xDigiA->calibTime() + xDigiB->calibTime() );
1221  //TODO: Afterpulse handling: correct hit time by the time difference between the first and second digi on the same side
1222  if(!mIsSim && mApCorr){//merge skip corrections for simulation
1223  time += t_corr_afterpulse;
1224  }
1225  // weight of merging of hits (later) is the total charge => sum of both ends ToT
1226  totSum = xDigiA->calibTot() + xDigiB->calibTot();
1227 
1228  // use local coordinates... (0,0,0) is in the center of counter
1229  posX = ( -1 * eTofConst::nStrips / 2. + strip - 0.5 ) * eTofConst::stripPitch;
1230 
1231  // check if the hit position (and hence time) is likely wrong due to a clock jumps in one of the Get4s
1232  bool hasClockJump = false;
1233  if( fabs( posY ) > 0.5 * ( eTofConst::coarseClockCycle * mSigVel.at( detIndex ) - eTofConst::stripLength ) * 0.9 ) {
1234  hasClockJump = true;
1235  }
1236 
1237 
1238  //check for position jumps -> get clock jumps
1239  bool leftjump = false;
1240  bool outsider = false;
1241  double dt = 0;
1242 
1243  if(xDigiA->side() == 2 ){
1244  dt = xDigiA->calibTime() - xDigiB->calibTime();
1245  }else{
1246  dt = xDigiB->calibTime() - xDigiA->calibTime();
1247  }
1248  if(abs(dt)> 8.5 ){
1249  outsider = true;
1250  }
1251  if(dt < 8.5 && dt > 0){
1252  leftjump = true;
1253  }
1254 
1255  //---------------------------------------------------------
1256  // correct for single side clock jumps. Sync signal recovers time jumps, so no double jumps *should* occur
1257  // it seems more likely that one Get4 misses one clock pulse and is 6.25ns behind, however jumps in both directions seem to be seen in the data
1258  // strategy: subtract half a clock tick from the hit time which shifts it to the right time or 6.25 ns too early
1259  // --> deal with too early hits in the matchMaker --> feed-back via mClockJumpDirection map
1260  //---------------------------------------------------------
1261  if( mDoClockJumpShift && hasClockJump ) {
1262  if( mDoQA ) {
1263  LOG_INFO << "shifted hit time in direction: " << mClockJumpDirection.at( detIndex * 10 + ( strip - 1 ) / 4 + 1 ) << endm;
1264  }
1265 
1266  time -= eTofConst::coarseClockCycle * 0.5 * mClockJumpDirection.at( detIndex * 10 + ( strip - 1 ) / 4 + 1 );
1267  timeDiff -= eTofConst::coarseClockCycle * ( ( timeDiff < 0 ) ? -1 : ( timeDiff > 0 ) );
1268 
1269  if(leftjump && mDoDoubleClockJumpShift){
1270  time += 2*(eTofConst::coarseClockCycle * 0.5 * mClockJumpDirection.at( detIndex * 10 + ( strip - 1 ) / 4 + 1 ));
1271  }
1272 
1273  // shift "uncorrectable" (improper RbR offsets, doublejumps, missmatches etc. ) hits far off in time and position to be sorted out later
1274  if(outsider && mDoDoubleClockJumpShift){
1275  time += 100*(eTofConst::coarseClockCycle * 0.5 * mClockJumpDirection.at( detIndex * 10 + ( strip - 1 ) / 4 + 1 ));
1276  timeDiff -= 100*(eTofConst::coarseClockCycle * ( ( timeDiff < 0 ) ? -1 : ( timeDiff > 0 ) )) + 2.0;
1277  }
1278 
1279  if( mDoQA ) {
1280  LOG_INFO << "shifted hit on: " << sector << "-" << plane << "-" << counter << " -- " << strip << " Direction " << mClockJumpDirection.at( detIndex * 10 + ( strip - 1 ) / 4 + 1 ) << endm;
1281  }
1282 
1283  if( xDigiA->side() == 2 ) { // recalculate Y-Position based on new time.
1284  posY = mSigVel.at( detIndex ) * timeDiff * 0.5;
1285  }
1286  else {
1287  posY = -1 * mSigVel.at( detIndex ) * timeDiff * 0.5;
1288  }
1289  }
1290 
1291 
1292  if( mDebug ) {
1293  LOG_DEBUG << "detIndex=" << detIndex << "posX=" << posX << " posY=" << posY << " time= " << time << " totSum=" << totSum << endm;
1294  }
1295 
1296  // build a hit (clustersize is one strip at this point)
1297  unsigned int clusterSize = 1;
1298  if( hasClockJump ) {
1299  // add 100 to the cluster size to identify jumped hits in the matchMaker
1300  clusterSize += 100;
1301  }
1302 
1303 
1304  //Modify individual counters if needed (e.g. flip in local y due to switched cables)
1305  if(mModMatrix.at(detIndex) > 0){
1306  int mode = mModMatrix.at(detIndex);
1307  modifyHit(mode, posX , posY , time);
1308  }
1309 
1310  StETofHit* constructedHit = new StETofHit( sector, plane, counter, time, totSum, clusterSize, posX, posY );
1311 
1312  //Check for "same direction double clockjumps" and update FlagMap
1313  if(mDoDoubleClockJumpShift){
1314  double starttime = startTime();
1315 
1316  if( starttime > 0){
1317 
1318  double tof = fmod( constructedHit->time() , eTofConst::bTofClockCycle ) - starttime;
1319  double exptof = 0;
1320 
1321  if(mETofGeom){
1322 
1323  StMuPrimaryVertex* pVtx = mMuDst->primaryVertex( 0 );
1324  StThreeVectorD posGlo = mETofGeom->hitLocal2Master( constructedHit );
1325 
1326  if( pVtx ) {
1327  StThreeVectorF vtxPos = pVtx->position();
1328  exptof = tofPathLength(&vtxPos , &posGlo , 0) / ( nanosecond * c_light );
1329  }
1330  }
1331 
1332  int get4Nr = detIndex * 10 + ( strip - 1 ) / 4 + 1;
1333 
1334  double tMin = mGet4doublejumpTmin;
1335  double t1 = mGet4doublejumpTimes.at(get4Nr).at(0);
1336  double t2 = mGet4doublejumpTimes.at(get4Nr).at(1);
1337 
1338  mGet4doublejumpTimes.at(get4Nr).erase(mGet4doublejumpTimes.at(get4Nr).begin());
1339  mGet4doublejumpTimes.at(get4Nr).push_back(tof - exptof);
1340  if(mGet4doublejumpTimes.at(get4Nr).size() > 2 ){
1341  mGet4doublejumpTimes.at(get4Nr).erase(mGet4doublejumpTimes.at(get4Nr).begin());
1342  }
1343 
1344  t1 = mGet4doublejumpTimes.at(get4Nr).at(0);
1345  t2 = mGet4doublejumpTimes.at(get4Nr).at(1);
1346 
1347  if(t2 < tMin && t1 < tMin && t2 > -999 && t1 > -999){
1348  mGet4doublejumpFlag.at(get4Nr) = 1;
1349  }
1350 
1351  if(t2 > tMin && t1 > tMin ){
1352  mGet4doublejumpFlag.at(get4Nr) = 0;
1353  }
1354 
1355  if(mGet4doublejumpFlag.at(get4Nr) == 1){
1356  constructedHit->setTime(constructedHit->time() + eTofConst::coarseClockCycle);
1357  constructedHit->setClusterSize(constructedHit->clusterSize() + 200);
1358  tof += eTofConst::coarseClockCycle;
1359  }
1360  }
1361  }
1362 
1363  // push hit into intermediate collection
1364  mStoreHit[ detIndex ].push_back( constructedHit );
1365 
1366  // fill pointer vector
1367  std::vector< unsigned int > containedDigiIndices;
1368 
1369  containedDigiIndices.push_back( mMapDigiIndex.at( xDigiA ) );
1370  containedDigiIndices.push_back( mMapDigiIndex.at( xDigiB ) );
1371 
1372  mMapHitDigiIndices[ constructedHit ] = containedDigiIndices;
1373 
1374  //reset afterpulse time correction!
1375  t_corr_afterpulse = 0.0;
1376 
1377  // pass IdTruth to hit
1378  if(mIsSim){
1379  int DigiIdA = xDigiA->rawTot();
1380  int DigiIdB = xDigiB->rawTot();
1381 
1382  if(DigiIdB==DigiIdA){constructedHit->setIdTruth(DigiIdA);}
1383  else{constructedHit->setIdTruth(9999);}
1384  }
1385 
1386  LOG_DEBUG << *( mStoreHit.at( detIndex ).back() ) << endm;
1387 
1388  // erase the two used digis!
1389  iterDigi = digiVec->begin();
1390  delete *iterDigi;
1391  delete *(iterDigi+1);
1392  digiVec->erase( iterDigi + 1 );
1393  digiVec->erase( iterDigi );
1394  if( mDoQA ){
1395  mHistograms.at( histNameDigisErased )->Fill( 6 );
1396  mHistograms.at( histNameDigisErased )->Fill( 6 );
1397  }
1398  } // end of loop over digis on the same strip
1399 
1400  } // end of loop over strips
1401 
1402  if( mDebug ) {
1403  LOG_DEBUG << "matchSides() - matching of sides done" << endm;
1404  }
1405 }//::matchSides
1406 
1407 
1408 //_____________________________________________________________
1409 void
1410 StETofHitMaker::fillUnclusteredHitQA( const double& tstart, const bool isMuDst )
1411 {
1412  if( !mDoQA ) {
1413  return;
1414  }
1415 
1416  // ---------------------------------------
1417  if( fabs( tstart + 9999. ) < 1.e-5 ) {
1418  LOG_WARN << "fillUnclusteredHitQA(): -- no valid start time available ... skip filling histograms with time of flight information" << endm;
1419  }
1420  // ---------------------------------------
1421 
1422  int nHitsPrinted = 0;
1423 
1424  for( const auto& kv : mStoreHit ) {
1425  unsigned int detIndex = kv.first;
1426 
1427  unsigned int sector = detIndex / 100;
1428  unsigned int plane = ( detIndex % 100 ) / 10;
1429  unsigned int counter = detIndex % 10;
1430 
1431  if( mDebug ) {
1432  LOG_DEBUG << sector << " " << plane << " " << counter << " size hitVec: " << kv.second.size() << endm;
1433  }
1434 
1435  for( const auto& hit : kv.second ) {
1436  if( mDebug ) {
1437  LOG_DEBUG << "matchSides() - HIT: ";
1438  LOG_DEBUG << "time=" << setprecision( 16 ) << hit->time() << " ";
1439  LOG_DEBUG << "totalTot=" << setprecision( 4 ) << hit->totalTot() << " ";
1440  LOG_DEBUG << "localX=" << setprecision( 4 ) << hit->localX() << " ";
1441  LOG_DEBUG << "localY=" << setprecision( 4 ) << hit->localY() << endm;
1442  }
1443 
1444  mCounterActive.at( ( sector - 13 ) * 9 + ( plane - 1 ) * 3 + ( counter - 1 ) ) = true;
1445 
1446  mHistograms.at( "unclusteredHit_tot" )->Fill( hit->totalTot() );
1447 
1448  std::string histNameTot = "unclusteredHit_tot_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
1449  mHistograms.at( histNameTot )->Fill( hit->totalTot(), hit->localX() );
1450 
1451 
1452 
1453  // ---------------------------------------
1454 
1455  if( mMapHitDigiIndices.at( hit ).size() >= 2 ) {
1456  int digiIndexA = mMapHitDigiIndices.at( hit ).at( 0 );
1457  int digiIndexB = mMapHitDigiIndices.at( hit ).at( 1 );
1458 
1459  float totA = 0.;
1460  float totB = 0.;
1461 
1462  bool sideSwitch = false;
1463 
1464  if( !isMuDst ) {
1465  StSPtrVecETofDigi& etofDigis = mEvent->etofCollection()->etofDigis();
1466  totA = etofDigis[ digiIndexA ]->calibTot();
1467  totB = etofDigis[ digiIndexB ]->calibTot();
1468 
1469  if( etofDigis[ digiIndexA ]->side() == 2 ) {
1470  sideSwitch = true;
1471  }
1472  }
1473  else {
1474  totA = mMuDst->etofDigi( digiIndexA )->calibTot();
1475  totB = mMuDst->etofDigi( digiIndexB )->calibTot();
1476 
1477  if( mMuDst->etofDigi( digiIndexA )->side() == 2 ) {
1478  sideSwitch = true;
1479  }
1480  }
1481 
1482  if( mDebug ) {
1483  LOG_DEBUG << "tot of digis in the hit: " << totA << " and " << totB << " sideSwitch: " << sideSwitch << endm;
1484  }
1485 
1486  float totDiff = totA - totB;
1487  if( sideSwitch ) totDiff *= -1;
1488 
1489  mHistograms.at( "unclusteredHit_tot_difference" )->Fill( totDiff );
1490 
1491  if( fabs( hit->localY() ) > mMaxYPos ) {
1492  mHistograms.at( "unclusteredHit_tail_totAsym" )->Fill( totDiff / ( totA + totB) );
1493  }
1494  }
1495 
1496  // ---------------------------------------
1497 
1498 
1499 
1500  std::string histNamePos = "unclusteredHit_pos_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
1501  mHistograms.at( histNamePos )->Fill( hit->localX(), hit->localY() );
1502 
1503  std::string histNamePosJump = "unclusteredHit_jump_pos_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
1504  if( hit->clusterSize() > 100 ) mHistograms.at( histNamePosJump )->Fill( hit->localX(), hit->localY() );
1505 
1506  // ---------------------------------------
1507  if( fabs( tstart + 9999. ) < 1.e-5 ) continue;
1508 
1509  double tof = fmod( hit->time(), eTofConst::bTofClockCycle ) - tstart;
1510 
1511  if( mDebug ) {
1512  LOG_DEBUG << "hit time, hit time mod bTOF clock cycle, start time, time difference: ";
1513  LOG_DEBUG << hit->time() << " , " << tof + tstart << " , " << tstart << " , " << tof << endm;
1514  LOG_DEBUG << "sector, plane, counter: " << hit->sector() << " , " << hit->zPlane() << " , " << hit->counter() << endm;
1515  }
1516  mHistograms.at( "unclusteredHit_tof_fullrange_all" )->Fill( tof );
1517 
1518 
1519  while( tof < 0. ) {
1520  tof += eTofConst::bTofClockCycle;
1521  }
1522  tof = fmod( tof, eTofConst::bTofClockCycle );
1523 
1524 
1525  mHistograms.at( "unclusteredHit_tof_fullrange" )->Fill( tof );
1526 
1527  mHistograms.at( "unclusteredHit_tof" )->Fill( tof );
1528 
1529 
1530  std::string histNameTof = "unclusteredHit_tof_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
1531  mHistograms.at( histNameTof )->Fill( hit->localX(), tof );
1532 
1533  std::string histNameTofJump = "unclusteredHit_jump_tof_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
1534  if( hit->clusterSize() > 100 ) mHistograms.at( histNameTofJump )->Fill( hit->localX(), tof );
1535 
1536 
1537  if( isMuDst ) {
1538  StMuPrimaryVertex* pVtx = mMuDst->primaryVertex( 0 );
1539  if( pVtx ) {
1540  StThreeVectorD vtxPos = pVtx->position();
1541 
1542  if( fabs( vtxPos.z() ) <= 10. && fabs( vtxPos.perp() < 2.5 ) ) {
1543  histNameTof = "unclusteredHit_pVtx_tof_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
1544  mHistograms.at( histNameTof )->Fill( hit->localX(), tof );
1545  }
1546  }
1547  }
1548 
1549 
1550 
1551  //if( fabs( hit->localY() ) > 25. ) {
1552  // histNameTof = "unclusteredHit_tail_tof_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
1553  // mHistograms.at( histNameTof )->Fill( hit->localX(), tof );
1554  //}
1555  //else {
1556  // histNameTof = "unclusteredHit_good_tof_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
1557  // mHistograms.at( histNameTof )->Fill( hit->localX(), tof );
1558  //}
1559 
1560 
1561  if( mDebug ) {
1562  if( fabs( tof ) > 1000. && nHitsPrinted < 5 ) {
1563  LOG_INFO << "TOF UNNORMALLY LARGE: " << tof << " !!! " << endm;
1564  nHitsPrinted++;
1565  }
1566  }
1567 
1568  // ---------------------------------------
1569  } // end of loop over hits in hitVec
1570  } // end of loop over mStoreHit
1571 }//::fillUnclusteredHitQA
1572 
1573 //_____________________________________________________________
1574 void
1575 StETofHitMaker::mergeClusters( const bool isMuDst )
1576 {
1577  std::vector< StETofHit* >::iterator iterHit;
1578  std::vector< std::vector< StETofDigi* > >::iterator iterDigi;
1579 
1580  for( auto kv = mStoreHit.begin(); kv != mStoreHit.end(); kv++ ) {
1581  unsigned int detIndex = kv->first;
1582 
1583  unsigned int sector = detIndex / 100;
1584  unsigned int plane = ( detIndex % 100 ) / 10;
1585  unsigned int counter = detIndex % 10;
1586 
1587  std::vector< StETofHit* > *hitVec = &( kv->second );
1588 
1589  size_t nHitsOnDet = 0;
1590 
1591  while( hitVec->size() > 0 ) {
1592  if( mDebug ) {
1593  LOG_DEBUG << "mergeClusters() - hit vector size: " << hitVec->size();
1594  LOG_DEBUG << " for sector " << sector << " plane " << plane << " counter " << counter << endm;
1595  LOG_DEBUG << "mergeClusters() - checking hit vector for possible hits to merge with..." << endm;
1596  }
1597 
1598  bool isMissmatch = false;
1599  int idTruth = 0;
1600 
1601  StETofHit* pHit = hitVec->at( 0 );
1602 
1603  // scale with tot for weigthed average
1604  double weight = pHit->totalTot();
1605 
1606  if( weight==0 ) {
1607  weight = 0.001;
1608  }
1609 
1610  double weightedTime = pHit->time() * weight;
1611  double weightedPosX = pHit->localX() * weight;
1612  double weightedPosY = pHit->localY() * weight;
1613  double weightsTotSum = 1. * weight;
1614 
1615  // currently only one-strip clusters --> lowest and highest strip identical
1616  unsigned int clusterSize = 1;
1617  int lowestStrip = ceil( pHit->localX() / eTofConst::stripPitch );
1618  int highestStrip = lowestStrip;
1619 
1620  bool hasClockJump = false;
1621  if( pHit->clusterSize() > 100 && pHit->clusterSize() < 999) {
1622  hasClockJump = true;
1623  }
1624 
1625  bool isSingleSided = false;
1626  if(pHit->clusterSize() > 999){
1627  isSingleSided = true;
1628  }
1629 
1630  unsigned int index = 1;
1631  while( hitVec->size() > 1 ) {
1632  if( mDebug ) {
1633  LOG_DEBUG << "mergeClusters() - index in hit vector = " << index << endm;
1634  }
1635  if( index >= hitVec->size() ) {
1636  if( mDebug ) {
1637  LOG_DEBUG << "mergeClusters() - loop index is exceeds size of hit vector -> stop looping" << endm;
1638  }
1639  break;
1640  }
1641 
1642  StETofHit* pMergeHit = hitVec->at( index );
1643 
1644  //calculate distance measures
1645  double timeDiff = pHit->time() - pMergeHit->time();
1646  double posYDiff = ( pHit->localY() - pMergeHit->localY() ) / mSigVel.at( detIndex ); // divide by signal velocity
1647 
1648  bool isHigherAdjacentStip = false;
1649  bool isLowerAdjacentStip = false;
1650 
1651  if( ceil( pMergeHit->localX() / eTofConst::stripPitch ) - highestStrip == 1 ) {
1652  isHigherAdjacentStip = true;
1653  }
1654  else if( ceil( pMergeHit->localX() / eTofConst::stripPitch ) - lowestStrip == -1 ) {
1655  isLowerAdjacentStip = true;
1656  }
1657 
1658  double MergingRadius = 0;
1659 
1660  // dont merge single sided matches here!! has to happen after matching!!
1661  if(pMergeHit->clusterSize() > 500 || pHit->clusterSize() > 500){
1662  MergingRadius = 0;
1663  }else{
1664  MergingRadius = mMergingRadius;
1665  }
1666 
1667  // check merging condition: X is not convoluted into the clusterbuilding radius
1668  // since it is not supposed to be zero --> check if X position is on a adjacent strip
1669  if( ( isHigherAdjacentStip || isLowerAdjacentStip ) &&
1670  ( sqrt( timeDiff * timeDiff + posYDiff * posYDiff ) ) < MergingRadius ) //
1671  {
1672  if( mDebug ) {
1673  LOG_DEBUG << "mergeClusters() - merging is going on" << endm;
1674  }
1675 
1676  if( mDoQA && clusterSize == 1 ) {
1677  std::string histName_unclusteredHit_delT_pos = "unclusteredHit_delT_pos_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
1678  std::string histName_unclusteredHit_delT_tot = "unclusteredHit_delT_tot_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
1679  std::string histName_unclusteredHit_delT = "unclusteredHit_delT_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
1680  mHistograms.at( histName_unclusteredHit_delT )->Fill( pHit->localX(), timeDiff );
1681  ((TH3F&) *mHistograms[ histName_unclusteredHit_delT_tot ]).TH3F::Fill( pHit->localX(), timeDiff, pHit->totalTot() );
1682  ((TH3F&) *mHistograms[ histName_unclusteredHit_delT_pos ]).TH3F::Fill( pHit->localX(), timeDiff, posYDiff );
1683  }
1684  //merge hit into cluster
1685  double hitWeight = pMergeHit->totalTot();
1686 
1687  if( hitWeight==0 ) {
1688  hitWeight = 0.001;
1689  }
1690 
1691  weightedTime += ( pMergeHit->time() * hitWeight );
1692  weightedPosX += ( pMergeHit->localX() * hitWeight );
1693  weightedPosY += ( pMergeHit->localY() * hitWeight );
1694  weightsTotSum += hitWeight;
1695 
1696  clusterSize++;
1697  if( pMergeHit->clusterSize() > 100 && pMergeHit->clusterSize() < 200) {
1698  hasClockJump = true;
1699  }
1700 
1701  if( isHigherAdjacentStip ) {
1702  highestStrip++;
1703  }
1704  else if( isLowerAdjacentStip ) {
1705  lowestStrip--;
1706  }
1707 
1708  if( mDebug ) {
1709  LOG_DEBUG << "mergeClusters() - detector: " << detIndex << " seed hit localX: " << pHit->localX();
1710  LOG_DEBUG << " <> hit to merge localX: " << pMergeHit->localX();
1711  LOG_DEBUG << " <> weighted localX: " << weightedPosX / weightsTotSum << endm;
1712  }
1713 
1714  // merge contained digi index vectors
1715  for( unsigned int i : mMapHitDigiIndices.at( pMergeHit ) ) {
1716  mMapHitDigiIndices.at( pHit ).push_back( i );
1717  }
1718 
1719 
1720  // merge IdTruth information
1721  if(mIsSim){
1722 
1723  if(pMergeHit->idTruth() !=pHit->idTruth()){
1724  idTruth = -99999;
1725  isMissmatch = 1;
1726  }
1727  else if(!isMissmatch){
1728  idTruth = pHit->idTruth();
1729  }
1730  }
1731 
1732 
1733  // erase the hit that was merged
1734  iterHit = hitVec->begin() + index;
1735  delete *iterHit;
1736  hitVec->erase( iterHit );
1737 
1738  if( mDebug ) {
1739  LOG_DEBUG << "mergeClusters() - deleting merged hit from Map" << endm;
1740  }
1741 
1742  mMapHitDigiIndices.erase( pMergeHit );
1743 
1744  }
1745  else { // merging condition not fulfilled
1746  if( mDebug ) {
1747  LOG_DEBUG << "mergeClusters() - merging condition not fulfilled -- check the next hit" << endm;
1748  }
1749  index++;
1750  } // check next hit
1751 
1752  } // end of loop over hits for merging
1753 
1754  // set idTruth if nothing to merge
1755  if(mIsSim){
1756  if(idTruth == 0){
1757  idTruth = pHit->idTruth();
1758  }
1759  }
1760 
1761  // renormalize with the total ToT
1762  weightedTime /= weightsTotSum;
1763  weightedPosX /= weightsTotSum;
1764  weightedPosY /= weightsTotSum;
1765 
1766  // use only the floating point remainder of the time with respect the the bTof clock range
1767  weightedTime = fmod( weightedTime, eTofConst::bTofClockCycle );
1768  if( weightedTime < 0 ) weightedTime += eTofConst::bTofClockCycle;
1769 
1770  if( hasClockJump ) {
1771  clusterSize += 100;
1772  }
1773 
1774  if(isSingleSided){
1775  clusterSize += 1000;
1776  }
1777 
1778  if( mDebug ) {
1779  LOG_DEBUG << "mergeClusters() - MERGED HIT: ";
1780  LOG_DEBUG << "sector: " << sector << " plane: " << plane << " counter: " << counter << "\n";
1781  LOG_DEBUG << "time = " << setprecision( 16 ) << weightedTime << " ";
1782  LOG_DEBUG << "totalTot = " << setprecision( 4 ) << weightsTotSum << " ";
1783  LOG_DEBUG << "clusterSize = " << setprecision( 1 ) << clusterSize % 100 << " ";
1784  LOG_DEBUG << "localX = " << setprecision( 4 ) << weightedPosX << " ";
1785  LOG_DEBUG << "localY = " << setprecision( 4 ) << weightedPosY << endm;
1786  }
1787 
1788  // create combined hit
1789  StETofHit* combinedHit = new StETofHit( sector, plane, counter, weightedTime, weightsTotSum, clusterSize, weightedPosX, weightedPosY );
1790 
1791  if(mIsSim){
1792  combinedHit->setIdTruth(idTruth);
1793  }
1794 
1795 
1796  // fill hit into the eTOF collection or the eTOf hit array depending on StEvent or MuDst input
1797  if( !isMuDst ) {
1798  mEvent->etofCollection()->addHit( combinedHit );
1799 
1800  // copy contained digis vector map over to the ETofCollection address
1801  mMapHitDigiIndices[ combinedHit ] = mMapHitDigiIndices.at( pHit );
1802 
1803  if( mDebug ) {
1804  LOG_DEBUG << "mergeClusters(): size of digi vector for combined hit " << nHitsOnDet << " on the counter: ";
1805  LOG_DEBUG << mMapHitDigiIndices.at( combinedHit ).size() << " copied over from merged hit vector of size "<< mMapHitDigiIndices.at( pHit ).size() << endm;
1806  }
1807  }
1808  else{
1809  mMuDst->addETofHit( ( StMuETofHit* ) combinedHit );
1810 
1811  int lastHitIndex = mMuDst->numberOfETofHit() - 1;
1812 
1813  // copy contained digis vector map over to the ETofCollection address
1814  mMapHitIndexDigiIndices[ lastHitIndex ] = mMapHitDigiIndices.at( pHit );
1815  if( mDebug ) {
1816  LOG_DEBUG << "mergeClusters(): size of digi vector for combined hit " << nHitsOnDet << " on the counter: ";
1817  LOG_DEBUG << mMapHitIndexDigiIndices.at( lastHitIndex ).size() << " copied over from merged hit vector of size "<< mMapHitDigiIndices.at( pHit ).size() << endm;
1818  }
1819  }
1820 
1821 
1822  // delete hit from the internal storage
1823  iterHit = hitVec->begin();
1824  delete *iterHit;
1825  hitVec->erase( iterHit );
1826 
1827  mMapHitDigiIndices.erase( pHit );
1828 
1829  if( mDebug ) {
1830  if( !isMuDst && mMapHitDigiIndices.count( combinedHit ) ) {
1831  LOG_DEBUG << "mergeClusters(): size of digi vector for combined hit " << nHitsOnDet << " on the counter: " << mMapHitDigiIndices.at( combinedHit ).size() << " after deletion of merged hit" << endm;
1832  }
1833  if( isMuDst && mMapHitIndexDigiIndices.count( nHitsOnDet ) ) {
1834  LOG_DEBUG << "mergeClusters(): size of digi vector for combined hit " << nHitsOnDet << " on the counter: " << mMapHitIndexDigiIndices.at( mMuDst->numberOfETofHit() - 1 ).size() << endm;
1835  }
1836  }
1837 
1838  if( isMuDst ) {
1839  //MuDst copies the combined hit over into new format. Not deleting the hits afterwards causes a memory leak as the created hits are never deleted. In StEvent, etofCollection->Clear() takes care of this.
1840  delete combinedHit;
1841  }
1842  nHitsOnDet++;
1843  } // end of loop over hits
1844 
1845  } // end of loop over detectors
1846 
1847  if( mDebug ) {
1848  LOG_DEBUG << "mergeClusters() - merging into clusters done" << endm;
1849  }
1850 }//::mergeClusters
1851 
1852 //_____________________________________________________________
1853 void
1854 StETofHitMaker::assignAssociatedHits( const bool isMuDst )
1855 {
1856  if( !isMuDst ) {
1857  StSPtrVecETofHit& etofHits = mEvent->etofCollection()->etofHits();
1858  StSPtrVecETofDigi& etofDigis = mEvent->etofCollection()->etofDigis();
1859 
1860  for( size_t ihit = 0; ihit < etofHits.size(); ihit++ ) {
1861  StETofHit* aHit = etofHits[ ihit ];
1862 
1863  if( mDebug ) {
1864  LOG_DEBUG << "assignAssociatedHits(): size of digi vector for hit " << ihit << ": "<< mMapHitDigiIndices.at( aHit ).size() << endm;
1865  }
1866 
1867  for ( size_t idigi = 0; idigi < mMapHitDigiIndices.at( aHit ).size(); idigi++ ) {
1868  if( mDebug ) {
1869  LOG_DEBUG << "assignAssociatedHits(): link points to digi " << etofDigis[ mMapHitDigiIndices.at( aHit ).at( idigi ) ] << endm;
1870  }
1871  etofDigis[ mMapHitDigiIndices.at( aHit ).at( idigi ) ]->setAssociatedHit( aHit );
1872  }
1873  }
1874 
1875  }
1876  else {
1877  for( const auto& kv : mMapHitIndexDigiIndices ) {
1878  if( mDebug ) {
1879  LOG_DEBUG << "assignAssociatedHits(): size of digi vector for hit index " << kv.first << ": "<< kv.second.size() << endm;
1880  }
1881 
1882  for( const auto& v: kv.second ) {
1883  if( mDebug ) {
1884  LOG_DEBUG << "assignAssociatedHits(): link to digi index " << v << endm;
1885  }
1886  mMuDst->etofDigi( v )->setAssociatedHitId( kv.first );
1887  }
1888  }
1889  }
1890 
1891  LOG_DEBUG << "assignAssociatedHits() - association between digis and hits done" << endm;
1892 }//::assignAssociatedHits
1893 
1894 
1895 
1896 //_____________________________________________________________
1897 void
1898 StETofHitMaker::fillHitQA( const bool isMuDst, const double& tstart )
1899 {
1900  int bTofCentral = 700;
1901 
1902  double vpdStart = -9999.;
1903  double vertexVz = -999.;
1904 
1905  if( mDoQA ) {
1906  startTimeVpd( vpdStart, vertexVz );
1907  }
1908 
1909  if( !isMuDst ) {
1910  // --------------------------------------------------------
1911  // analyze hits in eTOF
1912  // --------------------------------------------------------
1913  const StSPtrVecETofHit& etofHits = mEvent->etofCollection()->etofHits();
1914 
1915  int nHitsETof = 0;
1916  double averageETofHitTime = 0.;
1917 
1918  for( size_t i=0; i<etofHits.size(); i++ ) {
1919  StETofHit* aHit = etofHits[ i ];
1920 
1921  if( !aHit ) {
1922  continue;
1923  }
1924 
1925  if( mDebug ) {
1926  LOG_DEBUG << *aHit << endm;
1927  }
1928 
1929  updateCyclicRunningMean( aHit->time(), averageETofHitTime, nHitsETof, eTofConst::bTofClockCycle );
1930 
1931  // fill histogram to be saved in .hist.root file
1932 
1933  string histNamePos = "etofHit_pos_s" + std::to_string( aHit->sector() ) + "m" + std::to_string( aHit->zPlane() ) + "c" + std::to_string( aHit->counter() );
1934  mHistograms.at( histNamePos )->Fill( aHit->localX(), aHit->localY() );
1935 
1936  if( mDoQA ) {
1937  std::string histNameClustersize = "clustersize_s" + std::to_string( aHit->sector() ) + "m" + std::to_string( aHit->zPlane() );
1938  mHistograms.at( histNameClustersize )->Fill( aHit->clusterSize() % 100 );
1939  }
1940 
1941  // if tstart exists
1942  if( (fabs( tstart ) > 0.001 && fabs( tstart - ( eTofConst::bTofClockCycle - 9999. ) ) > 0.001) || mIsSim ) {
1943  double tof = aHit->time() - tstart;
1944  if( tof < -800 ) {
1945  tof += eTofConst::bTofClockCycle;
1946  }
1947 
1948  mHistograms.at( "etofHit_tof" )->Fill( tof );
1949  mHistograms.at( "etofHit_tof_fullrange" )->Fill( tof );
1950  }
1951 
1952  if( mDoQA ) {
1953  if( fabs( vpdStart - ( eTofConst::bTofClockCycle - 9999. ) ) > 0.001 ) {
1954  double tofVpd = aHit->time() - vpdStart;
1955  if( tofVpd < -800 ) {
1956  tofVpd += eTofConst::bTofClockCycle;
1957  }
1958  mHistograms.at( "etofHit_vpdVz_tof" )->Fill( vertexVz, tofVpd );
1959  }
1960  }
1961  }
1962 
1963  // --------------------------------------------------------
1964  // analyze hits in bTOF to get the eTOF-bTOF correlation
1965  // --------------------------------------------------------
1966  StBTofCollection* btofCollection = mEvent->btofCollection();
1967 
1968  if( !btofCollection || !btofCollection->hitsPresent() ) {
1969  LOG_WARN << "fillHitQA - no btof collection or no bTof hits present" << endm;
1970  return;
1971  }
1972 
1973  const StSPtrVecBTofHit& btofHits = btofCollection->tofHits();
1974 
1975  int nHitsBTof = 0;
1976  double averageBTofHitTime = 0.;
1977 
1978  for( size_t i=0; i<btofHits.size(); i++ ) {
1979  StBTofHit* aHit = btofHits[ i ];
1980 
1981  if( !aHit ) {
1982  continue;
1983  }
1984 
1985  updateCyclicRunningMean( aHit->leadingEdgeTime(), averageBTofHitTime, nHitsBTof, eTofConst::bTofClockCycle );
1986 
1987  // if doQA && tstart exists
1988  if( mDoQA && ((fabs( tstart ) > 0.001 && fabs( tstart - ( eTofConst::bTofClockCycle - 9999. ) ) > 0.001) || mIsSim) ) {
1989  double tof = aHit->leadingEdgeTime() - tstart;
1990  if( tof < 0 ) {
1991  tof += eTofConst::bTofClockCycle;
1992  }
1993  mHistograms.at( "btofHit_tof_fullrange" )->Fill( tof );
1994  }
1995  }
1996 
1997  float diff = averageETofHitTime - averageBTofHitTime;
1998  if( diff < -800 ) diff += eTofConst::bTofClockCycle;
1999  mHistograms.at( "averageTimeDiff_etofHits_btofHits" )->Fill( diff );
2000  mHistograms.at( "multiplicity_etofHits_btofHits" )->Fill( nHitsETof, nHitsBTof );
2001 
2002  if( mDoQA && nHitsBTof > bTofCentral ) {
2003  std::vector< int > etofHitsPerModule( eTofConst::nModules );
2004  for( size_t i=0; i<etofHits.size(); i++ ) {
2005  StETofHit* aHit = etofHits[ i ];
2006 
2007  if( !aHit ) {
2008  continue;
2009  }
2010  etofHitsPerModule.at( ( aHit->sector() - 13 ) * 3 + aHit->zPlane() - 1 ) += 1;
2011  }
2012 
2013  for( size_t i=0; i<eTofConst::nModules; i++ ) {
2014  mHistograms.at( "hitMultiplicityPerModuleCentral" )->Fill( i, etofHitsPerModule.at( i ) );
2015  }
2016  }
2017 
2018  // --------------------------------------------------------
2019  // analyze correlation with EPD East
2020  // --------------------------------------------------------
2021  StEpdCollection* epdCollection = mEvent->epdCollection();
2022  if( !epdCollection || !epdCollection->hitsPresent() ) {
2023  LOG_WARN << "fillHitQA - no epd collection or no epd hits present" << endm;
2024  return;
2025  }
2026 
2027  const StSPtrVecEpdHit& epdHits = epdCollection->epdHits();
2028 
2029  float nHitsEpdEast = 0.;
2030 
2031  for( size_t i=0; i<epdHits.size(); i++ ) {
2032  StEpdHit* epdHit = epdHits[ i ];
2033  if( !epdHit ) {
2034  continue;
2035  }
2036 
2037  if( epdHit->nMIP() < 0.3 ) continue;
2038  if( epdHit->id() > 0 ) continue; //positive id is the west side
2039 
2040  if( epdHit->nMIP() < 5 ) {
2041  nHitsEpdEast += epdHit->nMIP();
2042  }
2043  else {
2044  nHitsEpdEast += 5; //high hits are dominated by landau fluctuations
2045  }
2046  }
2047  mHistograms.at( "multiplicity_etofHits_epdEast" )->Fill( nHitsETof, nHitsEpdEast );
2048  if( mDoQA ) {
2049  mHistograms.at( "multiplicity_btofHits_epdEast" )->Fill( nHitsBTof, nHitsEpdEast );
2050  }
2051 
2052  }
2053  else {
2054  // --------------------------------------------------------
2055  // analyze hits in eTOF
2056  // --------------------------------------------------------
2057  int nHitsETof = 0;
2058  double averageETofHitTime = 0.;
2059 
2060  for( size_t i=0; i<mMuDst->numberOfETofHit(); i++ ) {
2061  StMuETofHit* aHit = mMuDst->etofHit( i );
2062 
2063  if( !aHit ) {
2064  continue;
2065  }
2066 
2067  if( mDebug ) {
2068  LOG_DEBUG << "hit (" << i << "): sector,plane,counter=" << aHit->sector() << ",";
2069  LOG_DEBUG << aHit->zPlane() << "," << aHit->counter() << " time=" << aHit->time();
2070  LOG_DEBUG << " localX=" << aHit->localX() << " localY=" << aHit->localY();
2071  LOG_DEBUG << " clustersize=" << aHit->clusterSize() % 100 << endm;
2072  }
2073 
2074  updateCyclicRunningMean( aHit->time(), averageETofHitTime, nHitsETof, eTofConst::bTofClockCycle );
2075 
2076  // fill histogram to be saved in .hist.root file
2077  string histNamePos = "etofHit_pos_s" + std::to_string( aHit->sector() ) + "m" + std::to_string( aHit->zPlane() ) + "c" + std::to_string( aHit->counter() );
2078  mHistograms.at( histNamePos )->Fill( aHit->localX(), aHit->localY() );
2079 
2080  if( mDoQA ) {
2081  std::string histNameClustersize = "clustersize_s" + std::to_string( aHit->sector() ) + "m" + std::to_string( aHit->zPlane() );
2082  mHistograms.at( histNameClustersize )->Fill( aHit->clusterSize() % 100 );
2083  }
2084 
2085  // if tstart exists
2086  if( (fabs( tstart ) > 0.001 && fabs( tstart - ( eTofConst::bTofClockCycle - 9999. ) ) > 0.001) || mIsSim ) {
2087  double tof = aHit->time() - tstart;
2088  if( tof < -800 ) {
2089  tof += eTofConst::bTofClockCycle;
2090  }
2091 
2092  mHistograms.at( "etofHit_tof" )->Fill( tof );
2093  mHistograms.at( "etofHit_tof_fullrange" )->Fill( tof );
2094  }
2095 
2096  if( mDoQA ) {
2097  if( fabs( vpdStart - ( eTofConst::bTofClockCycle - 9999. ) ) > 0.001 ) {
2098  double tofVpd = aHit->time() - vpdStart;
2099  if( tofVpd < -800 ) {
2100  tofVpd += eTofConst::bTofClockCycle;
2101  }
2102  mHistograms.at( "etofHit_vpdVz_tof" )->Fill( vertexVz, tofVpd );
2103  }
2104  }
2105  }
2106 
2107  // --------------------------------------------------------
2108  // analyze hits in bTOF to get the eTOF-bTOF correlation
2109  // --------------------------------------------------------
2110  if( !mMuDst->btofArray( muBTofHit ) || !mMuDst->numberOfBTofHit() ) {
2111  LOG_WARN << "fillHitQA - no btof hit array or no btof hits present" << endm;
2112  return;
2113  }
2114 
2115  int nHitsBTof = 0;
2116  double averageBTofHitTime = 0.;
2117 
2118  for( size_t i=0; i<mMuDst->numberOfBTofHit(); i++ ) {
2119  StMuBTofHit* aHit = mMuDst->btofHit( i );
2120 
2121  if( !aHit ) {
2122  continue;
2123  }
2124 
2125  updateCyclicRunningMean( aHit->leadingEdgeTime(), averageBTofHitTime, nHitsBTof, eTofConst::bTofClockCycle );
2126 
2127  // if doQA && tstart exists
2128  if( mDoQA && ((fabs( tstart ) > 0.001 && fabs( tstart - ( eTofConst::bTofClockCycle - 9999. ) ) > 0.001) || mIsSim) ) {
2129  double tof = aHit->leadingEdgeTime() - tstart;
2130  if( tof < -800 ) {
2131  tof += eTofConst::bTofClockCycle;
2132  }
2133 
2134  mHistograms.at( "btofHit_tof_fullrange" )->Fill( tof );
2135  }
2136  }
2137 
2138  double diff = averageETofHitTime - averageBTofHitTime;
2139  if( diff < -800 ) diff += eTofConst::bTofClockCycle;
2140 
2141  mHistograms.at( "averageTimeDiff_etofHits_btofHits" )->Fill( diff );
2142  mHistograms.at( "multiplicity_etofHits_btofHits" )->Fill( nHitsETof, nHitsBTof );
2143 
2144  if( mDoQA && nHitsBTof > bTofCentral ) {
2145  std::vector< int > etofHitsPerModule( eTofConst::nModules );
2146  for( size_t i=0; i<mMuDst->numberOfETofHit(); i++ ) {
2147  StMuETofHit* aHit = mMuDst->etofHit( i );
2148 
2149  if( !aHit ) {
2150  continue;
2151  }
2152  etofHitsPerModule.at( ( aHit->sector() - 13 ) * 3 + aHit->zPlane() - 1 ) += 1;
2153  }
2154 
2155  for( size_t i=0; i<eTofConst::nModules; i++ ) {
2156  mHistograms.at( "hitMultiplicityPerModuleCentral" )->Fill( i, etofHitsPerModule.at( i ) );
2157  }
2158  }
2159 
2160  // --------------------------------------------------------
2161  // analyze correlation with EPD East
2162  // --------------------------------------------------------
2163  if( !mMuDst->epdHits() || !mMuDst->numberOfEpdHit() ) {
2164  LOG_WARN << "fillHitQA - no epd hit array or no epd hits present" << endm;
2165  return;
2166  }
2167 
2168  size_t nHitsEpd = mMuDst->numberOfEpdHit();
2169  float nHitsEpdEast = 0.;
2170 
2171  for( size_t i=0; i<nHitsEpd; i++ ) {
2172  StMuEpdHit* epdHit = mMuDst->epdHit( i );
2173  if( !epdHit ) {
2174  continue;
2175  }
2176 
2177  if( epdHit->nMIP() < 0.3 ) continue;
2178  if( epdHit->id() > 0 ) continue; //positive id is the west side
2179 
2180  if( epdHit->nMIP() < 5 ) {
2181  nHitsEpdEast += epdHit->nMIP();
2182  }
2183  else {
2184  nHitsEpdEast += 5; //high hits are dominated by landau fluctuations
2185  }
2186  }
2187 
2188  mHistograms.at( "multiplicity_etofHits_epdEast" )->Fill( nHitsETof, nHitsEpdEast );
2189  if( mDoQA ) {
2190  mHistograms.at( "multiplicity_btofHits_epdEast" )->Fill( nHitsBTof, nHitsEpdEast );
2191  }
2192 
2193  }
2194 
2195  LOG_DEBUG << "fillHitQA() - histograms filled" << endm;
2196 }//::fillHitQA
2197 
2198 
2199 //_____________________________________________________________
2200 // setHistFileName uses the string argument from the chain being run to set
2201 // the name of the output histogram file.
2202 void
2203 StETofHitMaker::setHistFileName()
2204 {
2205  string extension = ".etofHit.root";
2206 
2207  if( GetChainOpt()->GetFileOut() != nullptr ) {
2208  TString outFile = GetChainOpt()->GetFileOut();
2209 
2210  mHistFileName = ( string ) outFile;
2211 
2212  // get rid of .root
2213  size_t lastindex = mHistFileName.find_last_of( "." );
2214  mHistFileName = mHistFileName.substr( 0, lastindex );
2215 
2216  // get rid of .MuDst or .event if necessary
2217  lastindex = mHistFileName.find_last_of( "." );
2218  mHistFileName = mHistFileName.substr( 0, lastindex );
2219 
2220  // get rid of directories
2221  lastindex = mHistFileName.find_last_of( "/" );
2222  mHistFileName = mHistFileName.substr( lastindex + 1, mHistFileName.length() );
2223 
2224  mHistFileName = mHistFileName + extension;
2225  } else {
2226  LOG_ERROR << "Cannot set the output filename for histograms" << endm;
2227  }
2228 }
2229 
2230 //_____________________________________________________________
2231 void
2232 StETofHitMaker::bookHistograms()
2233 {
2234  LOG_INFO << "bookHistograms() ... " << endm;
2235 
2236  mHistograms[ "etofHit_tof" ] = new TH1F( "etofHit_tof", "eTOF hit time of flight;time of flight (ns);# hits", 4000, -100., 150 );
2237  mHistograms[ "etofHit_tof_fullrange" ] = new TH1F( "etofHit_tof_fullrange", "eTOF hit time of flight;time of flight (ns);# hits", 5000, -800., eTofConst::bTofClockCycle );
2238  mHistograms[ "averageTimeDiff_etofHits_btofHits" ] = new TH1F( "averageTimeDiff_etofHits_btofHits", "difference between average times in bTOF and eTOF hits;#DeltaT (ns);# events", 4000, -500, 500 );
2239  mHistograms[ "multiplicity_etofHits_btofHits" ] = new TH2F( "multiplicity_etofHits_btofHits", "multiplicity correlation between bTOF and eTOF;# eTOF hits;# bTOF hits", 300, 0, 300, 500, 0, 1000 );
2240  mHistograms[ "multiplicity_etofHits_epdEast" ] = new TH2F( "multiplicity_etofHits_epdEast", "multiplicity correlation between eTOF and east EPD;# eTOF hits;# hits east EPD", 300, 0, 300, 200, 0, 1000 );
2241 
2242  AddHist( mHistograms.at( "etofHit_tof" ) );
2243  AddHist( mHistograms.at( "etofHit_tof_fullrange" ) );
2244  AddHist( mHistograms.at( "averageTimeDiff_etofHits_btofHits" ) );
2245  AddHist( mHistograms.at( "multiplicity_etofHits_btofHits" ) );
2246  AddHist( mHistograms.at( "multiplicity_etofHits_epdEast" ) );
2247 
2248  if( mDoQA ) {
2249  mHistograms[ "etofHit_vpdVz_tof" ] = new TH2F( "etofHit_vpdVz_tof", "eTOF hit time of flight;VPD Vz (cm);time of flight (ns)", 100, -200, 200, 1000, -50., 50 );
2250  mHistograms[ "btofHit_tof_fullrange" ] = new TH1F( "btofHit_tof_fullrange", "bTOF hit time of flight;time of flight (ns);# hits", 5000, -800., eTofConst::bTofClockCycle );
2251  mHistograms[ "multiplicity_btofHits_epdEast" ] = new TH2F( "multiplicity_btofHits_epdEast", "multiplicity correlation between bTOF and east EPD;# bTOF hits;# hits east EPD", 200, 0, 1000, 200, 0, 1000 );
2252  mHistograms[ "hitMultiplicityPerModuleCentral" ] = new TH2F( "hitMultiplicityPerModuleCentral", "hit multiplicity per module in central bTOF events", 36, 0, 36, 50, 0, 50 );
2253  mHistograms[ "multiplicity_etofDigis_etofHits" ] = new TH2F( "multiplicity_etofDigis_etofHits", "multiplicity correlation between eTOF digis and hits;# eTOF digis;# eTOF hits", 500, 0, 1000, 500, 0, 1000 );
2254  }
2255 
2256  for( int sector = eTofConst::sectorStart; sector <= eTofConst::sectorStop; sector++ ) {
2257  for( int plane = eTofConst::zPlaneStart; plane <= eTofConst::zPlaneStop; plane++ ) {
2258  for( int counter = eTofConst::counterStart; counter <= eTofConst::counterStop; counter++ ) {
2259  std::string histName_etofHit_pos = "etofHit_pos_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2260 
2261  mHistograms[ histName_etofHit_pos ] = new TH2F( Form( "etofHit_pos_s%d_m%d_c%d", sector, plane, counter ), "eTOF hit localXY;local X (cm);local Y (cm)", 64, -16., 16., 400, -500., 500. );
2262  AddHist( mHistograms.at( histName_etofHit_pos ) );
2263 
2264  if( mDoQA ) {
2265  std::string histNameDigisPerStrip = "digisPerStrip_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2266  mHistograms[ histNameDigisPerStrip ] = new TH2F( Form( "digisPerStrip_s%d_m%d_c%d", sector, plane, counter ), "digis per strip;strip;# digis per event", 32, 0.5, 32.5, 20, 0., 20. );
2267 
2268  std::string histNameDigisErased = "digisErased_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2269  mHistograms[ histNameDigisErased ] = new TH1F( Form( "digisErased_s%d_m%d_c%d", sector, plane, counter ), "digis erased;;# digis", 6, 0.5, 6.5 );
2270  mHistograms[ histNameDigisErased ]->GetXaxis()->SetBinLabel( 1, "digi inside dead time" );
2271  mHistograms[ histNameDigisErased ]->GetXaxis()->SetBinLabel( 2, "single digi on strip" );
2272  mHistograms[ histNameDigisErased ]->GetXaxis()->SetBinLabel( 3, "3 consecutive same side digis" );
2273  mHistograms[ histNameDigisErased ]->GetXaxis()->SetBinLabel( 4, "better match for partner" );
2274  mHistograms[ histNameDigisErased ]->GetXaxis()->SetBinLabel( 5, "single remaining digi" );
2275  mHistograms[ histNameDigisErased ]->GetXaxis()->SetBinLabel( 6, "matched into pair" );
2276  }
2277  }
2278 
2279  if( mDoQA ) {
2280  std::string histName_clustersize = "clustersize_s" + std::to_string( sector ) + "m" + std::to_string( plane );
2281  mHistograms[ histName_clustersize ] = new TH1F( Form( "clustersize_s%d_m%d", sector, plane ), "clustersize;clustersize;# events", 32, 0., 32. );
2282  }
2283  }
2284  }
2285 
2286  if( mDoQA ) {
2287  mHistograms[ "unclusteredHit_tot" ] = new TH1F( "unclusteredHit_tot", "unclustered hit tot; tot (ns);# unclustered hits", 1000, 0., 80. );
2288  mHistograms[ "unclusteredHit_tof" ] = new TH1F( "unclusteredHit_tof", "unclustered hit time of flight; time of flight (ns);# unclustered hits", 5000, 0., 1000. );
2289 
2290  mHistograms[ "unclusteredHit_tot_difference" ] = new TH1F( "unclusteredHit_tot_difference", "unclustered hit tot difference of digis; #Delta tot (ns);# unclustered hits", 1000, -100., 100. );
2291  mHistograms[ "unclusteredHit_tail_totAsym" ] = new TH1F( "unclusteredHit_tail_totAsym", "unclustered hit tail tot asymmetry of digis; #Delta tot/ tot sum;# unclustered hits", 200, -1., 1. );
2292 
2293  mHistograms[ "unclusteredHit_tof_fullrange" ] = new TH1F( "unclusteredHit_tof_fullrange", "unclustered hit time of flight; time of flight (ns);# unclustered hits", 5000, -1. * eTofConst::bTofClockCycle, eTofConst::bTofClockCycle );
2294  mHistograms[ "unclusteredHit_tof_fullrange_all" ] = new TH1F( "unclusteredHit_tof_fullrange_all", "unclustered hit time of flight; time of flight (ns);# unclustered hits", 5000, -1. * eTofConst::bTofClockCycle, eTofConst::bTofClockCycle );
2295 
2296  for( int sector = eTofConst::sectorStart; sector <= eTofConst::sectorStop; sector++ ) {
2297  for( int plane = eTofConst::zPlaneStart; plane <= eTofConst::zPlaneStop; plane++ ) {
2298  for( int counter = eTofConst::counterStart; counter <= eTofConst::counterStop; counter++ ) {
2299  std::string histName_unclusteredHit_tot = "unclusteredHit_tot_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2300  std::string histName_unclusteredHit_pos = "unclusteredHit_pos_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2301  std::string histName_unclusteredHit_tof = "unclusteredHit_tof_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2302  std::string histName_unclusteredHit_delT = "unclusteredHit_delT_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2303  std::string histName_unclusteredHit_delT_pos = "unclusteredHit_delT_pos_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2304  std::string histName_unclusteredHit_delT_tot = "unclusteredHit_delT_tot_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2305 
2306  std::string histName_unclusteredHit_tail_tof = "unclusteredHit_tail_tof_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2307  std::string histName_unclusteredHit_good_tof = "unclusteredHit_good_tof_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2308  std::string histName_unclusteredHit_pVtx_tof = "unclusteredHit_pVtx_tof_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2309 
2310  mHistograms[ histName_unclusteredHit_tot ] = new TH2F( Form( "unclusteredHit_tot_s%d_m%d_c%d", sector, plane, counter ), "unclustered hit tot; tot (ns);local X (cm)", 1000, 0., 80., 32, -16., 16. );
2311  mHistograms[ histName_unclusteredHit_pos ] = new TH2F( Form( "unclusteredHit_pos_s%d_m%d_c%d", sector, plane, counter ), "unclustered hit localXY; local X (cm);local Y (cm)", 32, -16., 16., 400, -50., 50. );
2312 
2313  mHistograms[ histName_unclusteredHit_tof ] = new TH2F( Form( "unclusteredHit_tof_s%d_m%d_c%d", sector, plane, counter ), "unclustered hit time of flight;local X (cm);time of flight (ns)", 32, -16., 16., 5000, 0., 1000. );
2314  mHistograms[ histName_unclusteredHit_delT ] = new TH2F( Form( "unclusteredHit_delT_s%d_m%d_c%d", sector, plane, counter ), "unclustered hit time to cluster neighbor;local X (cm); delT_cluster (ns)", 32, -16., 16., 200, -1., 1. );
2315 
2316  mHistograms[ histName_unclusteredHit_delT_pos ] = new TH3F( Form( "unclusteredHit_delT_DelPos_s%d_m%d_c%d", sector, plane, counter ), "unclustered hit time to cluster neighbor;local X (cm); delT_cluster (ns); y-pos (cm)", 32, -16., 16., 60, -0.3, 0.3, 20, -15., 15. );
2317 
2318  mHistograms[ histName_unclusteredHit_delT_tot ] = new TH3F( Form( "unclusteredHit_delT_tot_s%d_m%d_c%d", sector, plane, counter ), "unclustered hit time to cluster neighbor;local X (cm); delT_cluster (ns); tot (ns)", 32, -16., 16., 60, -0.3, 0.3, 25, 0., 10. );
2319 
2320  mHistograms[ histName_unclusteredHit_tail_tof ] = new TH2F( Form( "unclusteredHit_tof_tail_s%d_m%d_c%d", sector, plane, counter ), "unclustered hit (tail) time of flight;local X (cm);time of flight (ns)", 32, -16., 16., 5000, 0., 1000. );
2321  mHistograms[ histName_unclusteredHit_good_tof ] = new TH2F( Form( "unclusteredHit_tof_good_s%d_m%d_c%d", sector, plane, counter ), "unclustered hit (good) time of flight;local X (cm);time of flight (ns)", 32, -16., 16., 2500, 0., 50. );
2322  mHistograms[ histName_unclusteredHit_pVtx_tof ] = new TH2F( Form( "unclusteredHit_pVtx_tof_s%d_m%d_c%d", sector, plane, counter ), "unclustered hit time of flight;local X (cm);time of flight (ns)", 32, -16., 16., 2500, 0., 50. );
2323 
2324  std::string histName_digi_tot = "digi_tot_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2325 
2326  mHistograms[ histName_digi_tot ] = new TH2F( Form( "digi_tot_s%d_m%d_c%d", sector, plane, counter ), "digi tot;(strip-1)+0.5*(side-1);tot (arb. units)", 64, 0., 32., 200, 0., 40. );
2327 
2328  std::string histName_unclusteredHit_pos_time = "unclusteredHit_pos_time_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2329  mHistograms[ histName_unclusteredHit_pos_time ] = new TH2F( Form( "unclusteredHit_pos_time_s%d_m%d_c%d", sector, plane, counter ), "hit position over time;time (s);hit position (cm)", 3000, 56000., 59000., 200, -100., 100. );
2330 
2331 
2332  std::string histName_unclusteredHit_jump_pos = "unclusteredHit_jump_pos_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2333  std::string histName_unclusteredHit_jump_tof = "unclusteredHit_jump_tof_s" + std::to_string( sector ) + "m" + std::to_string( plane ) + "c" + std::to_string( counter );
2334 
2335  mHistograms[ histName_unclusteredHit_jump_pos ] = new TH2F( Form( "unclusteredHit_jump_pos_s%d_m%d_c%d", sector, plane, counter ), "unclustered jump hit localXY; local X (cm);local Y (cm)", 8, -16., 16., 200, -70., 70. );
2336  mHistograms[ histName_unclusteredHit_jump_tof ] = new TH2F( Form( "unclusteredHit_jump_tof_s%d_m%d_c%d", sector, plane, counter ), "unclustered jump hit time of flight;local X (cm);time of flight (ns)", 8, -16., 16., 400, 0., 50. );
2337  }
2338  }
2339  }
2340 
2341  mHistograms[ "counter_active" ] = new TH2F( "counter_active", "active counters by run;100*day + run number;counter", 15000, 5000., 20000., 108 , 0.5, 108.5 );
2342  }
2343 
2344  for ( auto& kv : mHistograms ) {
2345  kv.second->SetDirectory( 0 );
2346  }
2347 }
2348 
2349 //_____________________________________________________________
2350 void
2351 StETofHitMaker::writeHistograms()
2352 {
2353  if( mHistFileName != "" ) {
2354  LOG_INFO << "writing histograms to: " << mHistFileName.c_str() << endm;
2355 
2356  TFile histFile( mHistFileName.c_str(), "RECREATE", "etofHit" );
2357  histFile.cd();
2358 
2359  for ( const auto& kv : mHistograms ) {
2360  if( kv.second->GetEntries() > 0 ) kv.second->Write();
2361  }
2362 
2363  histFile.Close();
2364  }
2365  else {
2366  LOG_INFO << "histogram file name is empty string --> cannot write histograms" << endm;
2367  }
2368 }
2369 
2370 
2371 //_____________________________________________________________
2372 unsigned int
2373 StETofHitMaker::detectorToKey( const unsigned int detectorId )
2374 {
2375  unsigned int sector = ( detectorId / eTofConst::nCountersPerSector ) + eTofConst::sectorStart;
2376  unsigned int zPlane = ( ( detectorId % eTofConst::nCountersPerSector ) / eTofConst::nCounters ) + eTofConst::zPlaneStart;
2377  unsigned int counter = ( detectorId % eTofConst::nCounters ) + eTofConst::counterStart;
2378 
2379  return sector * 100 + zPlane * 10 + counter;
2380 }
2381 
2382 //_____________________________________________________________
2383 void
2384 StETofHitMaker::updateCyclicRunningMean( const double& value, double& mean, int& count, const double& range )
2385 {
2386  double valIn = value;
2387  if( mean - value < -0.9 * range ) {
2388  valIn -= range;
2389  }
2390  else if( mean - value > 0.9 * range ) {
2391  valIn += range;
2392  }
2393 
2394  count++;
2395 
2396  double scaling = 1. / count;
2397 
2398  if( mDebug ) {
2399  LOG_INFO << "old mean: " << mean << " scaling: " << scaling << " value in: " << valIn << endm;
2400  }
2401 
2402  mean = valIn * scaling + mean * ( 1. - scaling );
2403 
2404  if( mDebug ) {
2405  LOG_INFO << "new mean: " << mean << endm;
2406  }
2407 }
2408 
2409 //_____________________________________________________________
2410 void
2411 StETofHitMaker::updateClockJumpMap( const std::map< int, int >& clockJumpDir )
2412 {
2413  for( const auto& kv: clockJumpDir ) {
2414  mClockJumpDirection.at( kv.first ) = kv.second;
2415  }
2416 
2417  if( mDebug ) {
2418  for( const auto& kv : mClockJumpDirection ) {
2419  if( kv.second != 1 ) {
2420  LOG_INFO << "StETofHitMaker::updateClockJumpMap() -- entry in clock jump map: " << kv.first << " value: " << kv.second << endm;
2421  }
2422  }
2423  }
2424 }
2425 
2426 //_____________________________________________________________
2427 
2428 void
2429 StETofHitMaker::modifyHit( int modMode, double& localX,double& localY, double& time )
2430 {
2431  switch (modMode) {
2432 
2433  case 0:
2434  return;
2435 
2436  case 1:
2437  localX *= -1;
2438  localY *= -1;
2439  break;
2440 
2441  case 2:
2442  localX *= -1;
2443  break;
2444 
2445  case 3:
2446  localY *= -1;
2447  break;
2448 
2449  case 4:
2450  std::swap(localX, localY), localY = -localY;
2451  break;
2452 
2453  case 5:
2454  std::swap(localX, localY), localX = -localX;
2455  break;
2456 
2457  case 99:
2458  localX = 9999;
2459  localY = 9999;
2460  break;
2461  }
2462 }
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
StMuETofHit::zPlane
unsigned int zPlane() const
ZPlane.
Definition: StMuETofHit.h:194
StMuDst::btofHit
static StMuBTofHit * btofHit(int i)
returns pointer to the i-th muBTofHit
Definition: StMuDst.h:419
StMuETofHit::clusterSize
unsigned int clusterSize() const
Cluster size.
Definition: StMuETofHit.h:200
StMuETofHit::time
double time() const
Time.
Definition: StMuETofHit.h:197
StETofDigi::calibTime
double calibTime() const
calibrated time
Definition: StETofDigi.h:206
StETofDigi::zPlane
unsigned int zPlane() const
ZPlane.
Definition: StETofDigi.h:213
StMuDst::epdHits
static TClonesArray * epdHits()
returns pointer to the EpdHitCollection
Definition: StMuDst.h:297
StMuETofDigi::side
unsigned int side() const
Side.
Definition: StMuETofDigi.h:217
StMuETofHit::counter
unsigned int counter() const
Counter.
Definition: StMuETofHit.h:195
StMuEpdHit::id
Short_t id() const
Definition: StMuEpdHit.h:127
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
StMuDst::btofArray
static TClonesArray * btofArray(int type)
returns pointer to the n-th TClonesArray from the btof arrays // dongx
Definition: StMuDst.h:287
StETofHit::clusterSize
unsigned int clusterSize() const
Cluster size.
Definition: StETofHit.h:196
StETofDigi::sector
unsigned int sector() const
Sector.
Definition: StETofDigi.h:212
StEpdHit::nMIP
float nMIP() const
gain calibrated energy loss in tile, in units of Landau MPV for one MIP
Definition: StEpdHit.h:140
StETofDigi::strip
unsigned int strip() const
Strip.
Definition: StETofDigi.h:215
StETofDigi::rawTot
double rawTot() const
Getter for uncalibrated Tot.
Definition: StETofDigi.h:208
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
StEpdHit
Stores information for tiles in STAR Event Plane Detector.
Definition: StEpdHit.h:43
StETofHit::zPlane
unsigned int zPlane() const
ZPlane.
Definition: StETofHit.h:190
StMuDst::etofDigi
static StMuETofDigi * etofDigi(int i)
returns pointer to the i-th StMuEtofDigi
Definition: StMuDst.h:427
StETofHit::localX
double localX() const
X-position.
Definition: StETofHit.h:198
StETofHit::idTruth
unsigned int idTruth() const
mc-true associated track id in simulation
Definition: StETofHit.h:204
StMuDst::addETofHit
void addETofHit(const StMuETofHit *hit)
Definition: StMuDst.cxx:679
StETofHit::totalTot
double totalTot() const
Total Tot.
Definition: StETofHit.h:194
StETofHit::localY
double localY() const
Y-position.
Definition: StETofHit.h:199
StMuEpdHit::nMIP
Float_t nMIP()
gain calibrated energy loss in tile, in units of Landau MPV for one MIP
Definition: StMuEpdHit.h:81
StMuDst::btofHeader
static StBTofHeader * btofHeader()
returns pointer to the btofHeader - dongx
Definition: StMuDst.h:423
StETofHit::sector
unsigned int sector() const
Sector.
Definition: StETofHit.h:189
StETofDigi::counter
unsigned int counter() const
Counter.
Definition: StETofDigi.h:214
StETofDigi::calibTot
double calibTot() const
Getter for calibrated Tot.
Definition: StETofDigi.h:210
StETofDigi::side
unsigned int side() const
Side.
Definition: StETofDigi.h:217
StMuETofDigi::calibTot
double calibTot() const
Getter for calibrated Tot.
Definition: StMuETofDigi.h:210
StEpdHit::id
short id() const
Definition: StEpdHit.h:146
kStOk
Definition: Stypes.h:41
TDataSet::Find
virtual TDataSet * Find(const char *path) const
Definition: TDataSet.cxx:362
Generated by   doxygen 1.8.5