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