StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StETofCalibMaker.cxx
1  /***************************************************************************
2  *
3  * $Id: StETofCalibMaker.cxx,v 1.8 2020/01/16 03:10:33 fseck Exp $
4  *
5  * Author: Florian Seck, April 2018
6  ***************************************************************************
7  *
8  * Description: StETofCalibMaker - class to read the eTofCollection from
9  * StEvent, do calibration, fill the collection with the updated information
10  * and write back to StEvent
11  *
12  ***************************************************************************
13  *
14  * $Log: StETofCalibMaker.cxx,v $
15  * Revision 1.8 2020/01/16 03:10:33 fseck
16  * update handling of reset time for new cases in Run20
17  *
18  * Revision 1.7 2019/12/19 02:19:23 fseck
19  * use known pulser time differences inside one Gbtx to recover missing pulser signals
20  *
21  * Revision 1.6 2019/12/12 02:26:37 fseck
22  * ignore duplicate digis from stuck firmware
23  *
24  * Revision 1.5 2019/12/10 15:55:01 fseck
25  * added new database tables for pulsers, updated pulser handling and trigger time calculation
26  *
27  * Revision 1.4 2019/05/08 23:56:44 fseck
28  * change of default value for reference pulser
29  *
30  * Revision 1.3 2019/03/25 01:09:46 fseck
31  * added first version of pulser correction procedure + fix in reading parameters from db
32  *
33  * Revision 1.2 2019/03/08 19:01:07 fseck
34  * pick up the right trigger and reset time on event-by-event basis + fix to clearing of calibrated tot in afterburner mode + flag pulser digis
35  *
36  * Revision 1.1 2019/02/19 19:52:28 jeromel
37  * Reviewed code provided by F.Seck
38  *
39  *
40  ***************************************************************************/
41 #include <iostream>
42 #include <sstream>
43 #include <fstream>
44 #include <cmath>
45 #include <iterator>
46 #include <algorithm>
47 
48 #include "TString.h"
49 #include "TFile.h"
50 #include "TH1F.h"
51 #include "TH2F.h"
52 #include "TClass.h"
53 #include "TObject.h"
54 #include "TProfile.h"
55 
56 #include "StChain/StChainOpt.h" // for renaming the histogram file
57 
58 #include "StEvent.h"
59 
60 #include "StETofCollection.h"
61 #include "StETofHeader.h"
62 #include "StETofDigi.h"
63 
64 #include "StMuDSTMaker/COMMON/StMuDst.h"
65 #include "StMuDSTMaker/COMMON/StMuETofDigi.h"
66 #include "StMuDSTMaker/COMMON/StMuETofHeader.h"
67 
68 #include "StETofCalibMaker.h"
69 #include "StETofUtil/StETofHardwareMap.h"
70 #include "StETofUtil/StETofConstants.h"
71 
72 #include "tables/St_etofCalibParam_Table.h"
73 #include "tables/St_etofElectronicsMap_Table.h"
74 #include "tables/St_etofStatusMap_Table.h"
75 #include "tables/St_etofTimingWindow_Table.h"
76 #include "tables/St_etofSignalVelocity_Table.h"
77 #include "tables/St_etofDigiTotCorr_Table.h"
78 #include "tables/St_etofDigiTimeCorr_Table.h"
79 #include "tables/St_etofDigiSlewCorr_Table.h"
80 #include "tables/St_etofResetTimeCorr_Table.h"
81 #include "tables/St_etofPulserTotPeak_Table.h"
82 #include "tables/St_etofPulserTimeDiffGbtx_Table.h"
83 
84 namespace etofSlewing {
85  const unsigned int nTotBins = 30;
86 
87  // for space-efficient value in the database to nanoseconds
88  const float conversionFactor = 1. / 100.;
89 }
90 
91 
92 //_____________________________________________________________
94 : StMaker( "etofCalib", name ),
95  mEvent( nullptr ),
96  mMuDst( nullptr ),
97  mHwMap( nullptr ),
98  mFileNameCalibParam( "" ),
99  mFileNameElectronicsMap( "" ),
100  mFileNameStatusMap( "" ),
101  mFileNameTimingWindow( "" ),
102  mFileNameSignalVelocity( "" ),
103  mFileNameCalibHistograms( "" ),
104  mFileNameResetTimeCorr( "" ),
105  mFileNamePulserTotPeak( "" ),
106  mFileNamePulserTimeDiffGbtx( "" ),
107  mRunYear( 0 ),
108  mGet4TotBinWidthNs( 1. ),
109  mMinDigisPerSlewBin( 25 ),
110  mResetTimeCorr( 0. ),
111  mTriggerTime( 0. ),
112  mResetTime( 0. ),
113  mResetTs( 0 ),
114  mPulserPeakTime( 0. ),
115  mReferencePulserIndex( 0 ),
116  mStrictPulserHandling( 0 ),
117  mUsePulserGbtxDiff( true ),
118  mDoQA( false ),
119  mDebug( false ),
120  mHistFileName( "" )
121 {
123  LOG_DEBUG << "StETofCalibMaker::ctor" << endm;
124 
125  mStatus.clear();
126  mTimingWindow.clear();
127  mPulserWindow.clear();
128  mSignalVelocity.clear();
129  mDigiTotCorr.clear();
130  mDigiSlewCorr.clear();
131 
132  mPulserPeakTot.clear();
133  mPulserTimeDiff.clear();
134  mPulserTimeDiffGbtx.clear();
135  mNPulsersCounter.clear();
136  mNStatusBitsCounter.clear();
137  mPulserPresent.clear();
138 
139  mJumpingPulsers.clear();
140 
141  mUnlockPulserState.clear();
142 
143  mStuckFwDigi.clear();
144 
145  mHistograms.clear();
146 }
147 
148 
149 //_____________________________________________________________
150 StETofCalibMaker::~StETofCalibMaker()
151 { /* no op */
152 
153 }
154 
155 
156 //_____________________________________________________________
157 Int_t
158 StETofCalibMaker::Init()
159 {
160  LOG_INFO << "StETofCalibMaker::Init()" << endm;
161 
162  bookHistograms();
163 
164  return kStOk;
165 }
166 
167 
168 //_____________________________________________________________
169 Int_t
170 StETofCalibMaker::InitRun( Int_t runnumber )
171 {
172  mRunYear = ( runnumber + 727000 ) / 1000000 + 1999;
173  LOG_INFO << "runnumber: " << runnumber << " --> year: " << mRunYear << endm;
174 
175  TDataSet* dbDataSet = nullptr;
176  std::ifstream paramFile;
177 
178  // --------------------------------------------------------------------------------------------
179  // initialize calibration parameters from parameter file (if filename is provided) or database:
180  // -- electronics-to-hardware map
181  // -- status map
182  // -- timing window
183  // -- calib param
184  // -- signal velocity
185  // -- calib TOT factor per channel
186  // -- calib time offset per channel: position + T0
187  // -- slewing corrections
188  // -- reset time correction
189  // --------------------------------------------------------------------------------------------
190 
191  // electronics-to-hardware map
192  if( mFileNameElectronicsMap.empty() ) {
193  LOG_INFO << "etofElectronicsMap: no filename provided --> load database table" << endm;
194 
195  dbDataSet = GetDataBase( "Geometry/etof/etofElectronicsMap" );
196 
197  St_etofElectronicsMap* etofElectronicsMap = static_cast< St_etofElectronicsMap* > ( dbDataSet->Find( "etofElectronicsMap" ) );
198  if( !etofElectronicsMap ) {
199  LOG_ERROR << "unable to get the electronics map from the database" << endm;
200  return kStFatal;
201  }
202 
203  mHwMap = new StETofHardwareMap( etofElectronicsMap, mRunYear );
204  }
205  else {
206  LOG_INFO << "etofElectronicsMap: filename provided --> use parameter file: " << mFileNameElectronicsMap.c_str() << endm;
207 
208  mHwMap = new StETofHardwareMap( mFileNameElectronicsMap, mRunYear );
209  }
210 
211  // --------------------------------------------------------------------------------------------
212 
213  // status map
214  mStatus.clear();
215 
216  if( mFileNameStatusMap.empty() ) {
217  LOG_INFO << "etofStatusMap: no filename provided --> load database table" << endm;
218 
219  dbDataSet = GetDataBase( "Calibrations/etof/etofStatusMap" );
220 
221  St_etofStatusMap* etofStatusMap = static_cast< St_etofStatusMap* > ( dbDataSet->Find( "etofStatusMap" ) );
222  if( !etofStatusMap ) {
223  LOG_ERROR << "unable to get the status map from the database" << endm;
224  return kStFatal;
225  }
226 
227  etofStatusMap_st* statusMapTable = etofStatusMap->GetTable();
228 
229  for( size_t i=0; i<eTofConst::nChannelsInSystem; i++ ) {
230  mStatus[ channelToKey( i ) ] = (int) statusMapTable->status[ i ];
231  }
232  }
233  else {
234  LOG_INFO << "etofStatusMap: filename provided --> use parameter file: " << mFileNameStatusMap.c_str() << endm;
235 
236  paramFile.open( mFileNameStatusMap.c_str() );
237 
238  if( !paramFile.is_open() ) {
239  LOG_ERROR << "unable to get the 'etofStatusMap' parameters from file --> file does not exist" << endm;
240  return kStFatal;
241  }
242 
243  std::vector< int > status;
244  int temp;
245  while( paramFile >> temp ) {
246  status.push_back( temp );
247  }
248 
249  paramFile.close();
250 
251  if( status.size() != eTofConst::nChannelsInSystem ) {
252  LOG_ERROR << "parameter file for 'etofStatusMap' has not the right amount of entries: ";
253  LOG_ERROR << status.size() << " instead of " << eTofConst::nChannelsInSystem << " !!!!" << endm;
254  return kStFatal;
255  }
256 
257  for( size_t i=0; i<eTofConst::nChannelsInSystem; i++ ) {
258  mStatus[ channelToKey( i ) ] = status[ i ];
259  }
260  }
261 
262  for( const auto& kv : mStatus ) {
263  LOG_DEBUG << "channel key: " << kv.first << " --> status = " << kv.second << endm;
264  }
265 
266  // --------------------------------------------------------------------------------------------
267 
268  // timing window
269  mTimingWindow.clear();
270  mPulserWindow.clear();
271 
272  if( mFileNameTimingWindow.empty() ) {
273  LOG_INFO << "etofTimingWindow: no filename provided --> load database table" << endm;
274 
275  dbDataSet = GetDataBase( "Calibrations/etof/etofTimingWindow" );
276 
277  St_etofTimingWindow* etofTimingWindow = static_cast< St_etofTimingWindow* > ( dbDataSet->Find( "etofTimingWindow" ) );
278  if( !etofTimingWindow ) {
279  LOG_ERROR << "unable to get the timing window from the database" << endm;
280  return kStFatal;
281  }
282 
283  etofTimingWindow_st* timingWindowTable = etofTimingWindow->GetTable();
284 
285  for( size_t i=0; i<12 ; i++ ) {
286  int key = timingWindowTable->afckAddress[ i ];
287  if( key > 0 ) {
288  mTimingWindow[ key ] = std::make_pair( timingWindowTable->timingMin[ i ], timingWindowTable->timingMax[ i ] );
289  mPulserWindow[ key ] = std::make_pair( timingWindowTable->pulserMin[ i ], timingWindowTable->pulserMax[ i ] );
290 
291  mPulserPeakTime = timingWindowTable->pulserPeak[ i ];
292  }
293  }
294  }
295  else {
296  LOG_INFO << "etofTimingWindow: filename provided --> use parameter file: " << mFileNameTimingWindow.c_str() << endm;
297 
298  paramFile.open( mFileNameTimingWindow.c_str() );
299 
300  if( !paramFile.is_open() ) {
301  LOG_ERROR << "unable to get the 'etofTimingWindow' parameters from file --> file does not exist" << endm;
302  return kStFatal;
303  }
304 
305  std::string temp;
306  std::vector< float > times;
307  int address;
308  float time;
309 
310  while( std::getline( paramFile, temp ) ) {
311  std::istringstream istring( temp );
312  times.clear();
313 
314  istring >>std::hex >> address >> std::dec;
315 
316  while( istring >> time ) {
317  times.push_back( time );
318  }
319 
320  if( times.size() != 6 ) {
321  LOG_ERROR << "parameter file for 'etofTimingWindow' has not the right amount of entries: ";
322  LOG_ERROR << "times for address (0x" << std::hex << address << std::dec << ") has " << times.size() << " instead of " << 6 << " !!!!" << endm;
323  return kStFatal;
324  }
325 
326  if( address > 0 ) {
327  mTimingWindow[ address ] = std::make_pair( times.at( 0 ), times.at( 1 ) );
328  mPulserWindow[ address ] = std::make_pair( times.at( 3 ), times.at( 4 ) );
329 
330  mPulserPeakTime = times.at( 5 );
331  }
332  }
333  paramFile.close();
334 
335  if( mTimingWindow.size() > 12 ) {
336  LOG_ERROR << " too many entries in mTimingWindow map ...." << endm;
337  return kStFatal;
338  }
339  if( mPulserWindow.size() > 12 ) {
340  LOG_ERROR << " too many entries in mPulserWindow map ...." << endm;
341  return kStFatal;
342  }
343  }
344 
345  for( const auto& kv : mTimingWindow ) {
346  LOG_DEBUG << "AFCK address: 0x" << std::hex << kv.first << std::dec << " --> timing window from " << kv.second.first << " to " << kv.second.second << " ns" << endm;
347  }
348  for( const auto& kv : mPulserWindow ) {
349  LOG_DEBUG << "AFCK address: 0x" << std::hex << kv.first << std::dec << " --> pulser window from " << kv.second.first << " to " << kv.second.second << " ns" << endm;
350  }
351  LOG_DEBUG << "pulser time peak at " << mPulserPeakTime << " ns" << endm;
352  // --------------------------------------------------------------------------------------------
353 
354  // calib param
355  if( mFileNameCalibParam.empty() ) {
356  LOG_INFO << "etofCalibParam: no filename provided --> load database table" << endm;
357 
358  dbDataSet = GetDataBase( "Calibrations/etof/etofCalibParam" );
359 
360  St_etofCalibParam* etofCalibParam = static_cast< St_etofCalibParam* > ( dbDataSet->Find( "etofCalibParam" ) );
361  if( !etofCalibParam ) {
362  LOG_ERROR << "unable to get the calibration params from the database" << endm;
363  return kStFatal;
364  }
365 
366  etofCalibParam_st* calibParamTable = etofCalibParam->GetTable();
367 
368  mGet4TotBinWidthNs = calibParamTable->get4TotBinWidthNs;
369  mMinDigisPerSlewBin = calibParamTable->minDigisInSlewBin;
370 
371  // only set the reference pulser index if it is not alredy set from outside by a steering macro
372  if( mReferencePulserIndex == 0 ) {
373  mReferencePulserIndex = calibParamTable->referencePulserIndex;
374  }
375  else {
376  LOG_INFO << "--- reference pulser index is set manually ---" << endm;
377  }
378  }
379  else {
380  LOG_INFO << "etofCalibParam: filename provided --> use parameter file: " << mFileNameCalibParam.c_str() << endm;
381 
382  paramFile.open( mFileNameCalibParam.c_str() );
383 
384  if( !paramFile.is_open() ) {
385  LOG_ERROR << "unable to get the 'etofCalibParam' parameters from file --> file does not exist" << endm;
386  return kStFatal;
387  }
388 
389  std::vector< float > param;
390  float temp;
391  while( paramFile >> temp ) {
392  param.push_back( temp );
393  }
394 
395  paramFile.close();
396 
397  if( param.size() != 3 ) {
398  LOG_ERROR << "parameter file for 'etofCalibParam' has not the right amount of entries: ";
399  LOG_ERROR << param.size() << " instead of 3 !!!!" << endm;
400  return kStFatal;
401  }
402 
403  if( param.at( 0 ) > 0. ) {
404  mGet4TotBinWidthNs = param.at( 0 );
405  }
406  if( param.at( 1 ) > 0 ) {
407  mMinDigisPerSlewBin = param.at( 1 );
408  }
409 
410  // only set the reference pulser index if it is not alredy set from outside by a steering macro
411  if( param.at( 2 ) > 0 && mReferencePulserIndex == 0 ) {
412  mReferencePulserIndex = param.at( 2 );
413  }
414  else {
415  LOG_INFO << "--- reference pulser index is set manually ---" << endm;
416  }
417  }
418 
419  LOG_INFO << " Get4 TOT bin width to ns conversion factor: " << mGet4TotBinWidthNs << endm;
420  LOG_INFO << " minimal number of digis required in slewing bin: " << mMinDigisPerSlewBin << endm;
421  LOG_INFO << " reference pulser index: " << mReferencePulserIndex << endm;
422 
423 
424  // --------------------------------------------------------------------------------------------
425 
426  // signal velocities
427  mSignalVelocity.clear();
428 
429  if( mFileNameSignalVelocity.empty() ) {
430  LOG_INFO << "etofSignalVelocity: no filename provided --> load database table" << endm;
431 
432  dbDataSet = GetDataBase( "Calibrations/etof/etofSignalVelocity" );
433 
434  St_etofSignalVelocity* etofSignalVelocity = static_cast< St_etofSignalVelocity* > ( dbDataSet->Find( "etofSignalVelocity" ) );
435  if( !etofSignalVelocity ) {
436  LOG_ERROR << "unable to get the signal velocity from the database" << endm;
437  return kStFatal;
438  }
439 
440  etofSignalVelocity_st* velocityTable = etofSignalVelocity->GetTable();
441 
442  for( size_t i=0; i<eTofConst::nCountersInSystem; i++ ) {
443  if( velocityTable->signalVelocity[ i ] > 0 ) {
444  mSignalVelocity[ detectorToKey( i ) ] = velocityTable->signalVelocity[ i ];
445  }
446  }
447  }
448  else {
449  LOG_INFO << "etofSignalVelocity: filename provided --> use parameter file: " << mFileNameSignalVelocity.c_str() << endm;
450 
451  paramFile.open( mFileNameSignalVelocity.c_str() );
452 
453  if( !paramFile.is_open() ) {
454  LOG_ERROR << "unable to get the 'etofSignalVelocity' parameters from file --> file does not exist" << endm;
455  return kStFatal;
456  }
457 
458  std::vector< float > velocity;
459  float temp;
460  while( paramFile >> temp ) {
461  velocity.push_back( temp );
462  }
463 
464  paramFile.close();
465 
466  if( velocity.size() != eTofConst::nCountersInSystem ) {
467  LOG_ERROR << "parameter file for 'etofSignalVelocity' has not the right amount of entries: ";
468  LOG_ERROR << velocity.size() << " instead of " << eTofConst::nCountersInSystem << " !!!!" << endm;
469  return kStFatal;
470  }
471 
472  for( size_t i=0; i<eTofConst::nCountersInSystem; i++ ) {
473  if( velocity.at( i ) > 0 ) {
474  mSignalVelocity[ detectorToKey( i ) ] = velocity.at( i );
475  }
476  }
477  }
478 
479  for( const auto& kv : mSignalVelocity ) {
480  LOG_DEBUG << "counter key: " << kv.first << " --> signal velocity = " << kv.second << " cm / ns" << endm;
481  }
482 
483  // --------------------------------------------------------------------------------------------
484 
485  // digi tot correction factor, time correction offset and slewing corrections
486  mDigiTotCorr.clear();
487  mDigiTimeCorr.clear();
488  mDigiSlewCorr.clear();
489 
490  if( mFileNameCalibHistograms.empty() ) {
491  //-------------------
492  // digi tot corr
493  //-------------------
494  LOG_INFO << "etofDigiTotCorr: no filename provided --> load database table" << endm;
495 
496  dbDataSet = GetDataBase( "Calibrations/etof/etofDigiTotCorr" );
497 
498  St_etofDigiTotCorr* etofDigiTotCorr = static_cast< St_etofDigiTotCorr* > ( dbDataSet->Find( "etofDigiTotCorr" ) );
499  if( !etofDigiTotCorr ) {
500  LOG_ERROR << "unable to get the digi tot correction from the database" << endm;
501  return kStFatal;
502  }
503 
504  etofDigiTotCorr_st* digiTotCorrTable = etofDigiTotCorr->GetTable();
505 
506  for( size_t i=0; i<eTofConst::nChannelsInSystem; i++ ) {
507  unsigned int key = channelToKey( i );
508  unsigned int detector = key / 1000;
509 
510  if( mDigiTotCorr.count( detector ) == 0 ) {
511  TString name = Form( "digiTotCorr_%d", detector );
512  mDigiTotCorr[ detector ] = new TH1F( name, name, 2 * eTofConst::nStrips, 0, 2 * eTofConst::nStrips );
513  }
514 
515  unsigned int strip = ( key % 1000 ) / 10;
516  unsigned int side = key % 10;
517 
518  if( mDebug ) {
519  LOG_DEBUG << i << " " << detector << " " << strip << " " << side << " " << digiTotCorrTable->totCorr[ i ] << endm;
520  }
521 
522  mDigiTotCorr.at( detector )->SetBinContent( strip + eTofConst::nStrips * ( side - 1 ) , digiTotCorrTable->totCorr[ i ] );
523  }
524 
525  for( auto& kv : mDigiTotCorr ) {
526  kv.second->SetDirectory( 0 );
527  }
528 
529 
530  //-------------------
531  // digi time corr
532  //-------------------
533  LOG_INFO << "etofDigiTimeCorr: no filename provided --> load database table" << endm;
534 
535  // (1) position offset + (2) T0 offset
536  dbDataSet = GetDataBase( "Calibrations/etof/etofDigiTimeCorr" );
537 
538  St_etofDigiTimeCorr* etofDigiTimeCorr = static_cast< St_etofDigiTimeCorr* > ( dbDataSet->Find( "etofDigiTimeCorr" ) );
539  if( !etofDigiTimeCorr ) {
540  LOG_ERROR << "unable to get the digi time correction from the database" << endm;
541  return kStFatal;
542  }
543 
544  etofDigiTimeCorr_st* digiTimeCorrTable = etofDigiTimeCorr->GetTable();
545 
546  for( size_t i=0; i<eTofConst::nChannelsInSystem; i++ ) {
547  unsigned int key = channelToKey( i );
548  unsigned int detector = key / 1000;
549 
550  if( mDigiTimeCorr.count( detector ) == 0 ) {
551  TString name = Form( "digiTimeCorr_%d", detector );
552  mDigiTimeCorr[ detector ] = new TH1F( name, name, 2 * eTofConst::nStrips, 0, 2 * eTofConst::nStrips );
553  }
554 
555  unsigned int strip = ( key % 1000 ) / 10;
556  unsigned int side = key % 10;
557 
558  if( mDebug ) {
559  LOG_DEBUG << i << " " << detector << " " << strip << " " << side << " " << digiTimeCorrTable->timeCorr[ i ] << endm;
560  }
561 
562  mDigiTimeCorr.at( detector )->SetBinContent( strip + eTofConst::nStrips * ( side - 1 ) , digiTimeCorrTable->timeCorr[ i ] );
563  }
564 
565  for( auto& kv : mDigiTimeCorr ) {
566  kv.second->SetDirectory( 0 );
567  }
568 
569 
570  //-------------------
571  // digi slewing corr
572  //-------------------
573  LOG_INFO << "etofDigiSlewCorr: no filename provided --> load database table" << endm;
574 
575  dbDataSet = GetDataBase( "Calibrations/etof/etofDigiSlewCorr" );
576 
577  St_etofDigiSlewCorr* etofDigiSlewCorr = static_cast< St_etofDigiSlewCorr* > ( dbDataSet->Find( "etofDigiSlewCorr" ) );
578  if( !etofDigiSlewCorr ) {
579  LOG_ERROR << "unable to get the digi slewing correction from the database" << endm;
580  return kStFatal;
581  }
582 
583  etofDigiSlewCorr_st* digiSlewCorrTable = etofDigiSlewCorr->GetTable();
584 
585  for( size_t i=0; i<eTofConst::nChannelsInSystem; i++ ) {
586  unsigned int key = channelToKey( i );
587 
588  // fill for each channel the tot-bin edges and correction values into a vector
589  std::vector< float > binEdges;
590  std::vector< float > corr;
591 
592  binEdges.push_back( 0. );
593  for( size_t j=0; j<etofSlewing::nTotBins; j++ ) {
594  int tableIndex = i * etofSlewing::nTotBins + j;
595 
596  if( i != digiSlewCorrTable->channelId[ tableIndex ] ) {
597  LOG_ERROR << "something went wrong in reading the database table for eTOF slewing corrections" << endm;
598  return kStFatal;
599  }
600 
601  binEdges.push_back( digiSlewCorrTable->upperTotEdge[ tableIndex ] * etofSlewing::conversionFactor );
602  corr.push_back( digiSlewCorrTable->corr[ tableIndex ] * etofSlewing::conversionFactor );
603  }
604 
605  // create slewing correction TProfile histograms
606  if( mDigiSlewCorr.count( key ) == 0 ) {
607  TString name = Form( "digiSlewCorr_%d", key );
608  mDigiSlewCorr[ key ] = new TProfile( name, name, etofSlewing::nTotBins, binEdges.data() );
609  }
610 
611  // fill the histograms with weight, so that GetBinEnties( bin ) is larger that mMinDigisPerSlewBin
612  for( size_t j=0; j<etofSlewing::nTotBins; j++ ) {
613  float totCenter = mDigiSlewCorr.at( key )->GetBinCenter( j+1 );
614  mDigiSlewCorr.at( key )->Fill( totCenter, corr.at( j ), mMinDigisPerSlewBin + 1 );
615  }
616  }
617 
618  for( auto& kv : mDigiSlewCorr ) {
619  kv.second->SetDirectory( 0 );
620  }
621  }
622  else {//input from file
623  LOG_INFO << "etof calibration histograms: filename provided --> use parameter file: " << mFileNameCalibHistograms.c_str() << endm;
624 
625  if( !isFileExisting( mFileNameCalibHistograms ) ) {
626  LOG_ERROR << "unable to get the 'etofDigiTotCorr', 'etofDigiTimeCorr', 'etofDigiSlewCorr' parameters from file --> file does not exist" << endm;
627  }
628 
629  TFile* histFile = new TFile( mFileNameCalibHistograms.c_str(), "READ" );
630 
631  if( !histFile ) {
632  LOG_ERROR << "No calibration file found: " << mFileNameCalibHistograms.c_str() << endm;
633  LOG_INFO << "setting all parameters to default" << endm;
634  }else if( histFile->IsZombie() ){
635  LOG_ERROR << "Zombie calibration file: " << mFileNameCalibHistograms.c_str() << endm;
636  LOG_INFO << "stopping execution" << endm;
637  return kStFatal;
638  }
639 
640  TFile* histOffsetFile = new TFile( mFileNameOffsetHistograms.c_str(), "READ" ); //create setter!
641 
642  if( !histOffsetFile ) {
643  LOG_INFO << "No offset file found: " << mFileNameOffsetHistograms.c_str() << endm;
644  LOG_INFO << "setting all parameters to default" << endm;
645  }else if( histOffsetFile->IsZombie() ) {
646  LOG_ERROR << "Zombie offset file: " << mFileNameOffsetHistograms.c_str() << endm;
647  LOG_INFO << "stopping execution" << endm;
648  return kStFatal;
649  }else{
650  LOG_INFO << "Successfully opened RunOffset file "<< mFileNameOffsetHistograms.c_str() << endm;
651  }
652 
653 
654 
655  TString hPosOffsetName = Form( "calib_Run%d_PosOffsets_pfx", runnumber );
656  TString hT0OffsetName = Form( "calib_Run%d_T0Offsets_pfx", runnumber );
657  TProfile* hPosOffsetProfile = nullptr;
658  TProfile* hT0OffsetProfile = nullptr;
659 
660  if( histOffsetFile && !(histOffsetFile->IsZombie()) ){
661  LOG_INFO << "Getting run offset histograms Pos: "<< hPosOffsetName << " T0: "<< hT0OffsetName << endm;
662  hPosOffsetProfile = ( TProfile* ) histOffsetFile->Get( hPosOffsetName );
663  hT0OffsetProfile = ( TProfile* ) histOffsetFile->Get( hT0OffsetName );
664  }
665 
666  for( unsigned int sector = eTofConst::sectorStart; sector <= eTofConst::sectorStop; sector++ ) {
667 
668  if( mRunYear == 2018 && sector != 18 ) continue;
669 
670  for( unsigned int zPlane = eTofConst::zPlaneStart; zPlane <= eTofConst::zPlaneStop; zPlane++ ) {
671  for( unsigned int counter = eTofConst::counterStart; counter <= eTofConst::counterStop; counter++ ) {
672 
673  unsigned int key = sector * 100 + zPlane * 10 + counter;
674 
675  LOG_DEBUG << "detectorId key: " << sector << " " << zPlane << " " << counter << endm;
676 
677  TString hname;
678  TProfile* hProfile;
679 
680 
681  //-------------------
682  // digi tot corr
683  //-------------------
684  if( mDigiTotCorr.count( key ) == 0 ) {
685  TString name = Form( "digiTotCorr_%d", key );
686  mDigiTotCorr[ key ] = new TH1F( name, name, 2 * eTofConst::nStrips, 0, 2 * eTofConst::nStrips );
687  }
688 
689  hname = Form( "calib_Sector%02d_ZPlane%d_Det%d_TotDigi_pfx", sector, zPlane, counter );
690  hProfile = ( TProfile* ) histFile->Get( hname );
691 
692  if( hProfile ) {
693  for( size_t i=1; i<=2 * eTofConst::nStrips; i++ ) {
694  if( hProfile->GetBinContent( i ) != 0 ) {
695  mDigiTotCorr.at( key )->SetBinContent( i , hProfile->GetBinContent( i ) );
696  }
697  else {
698  mDigiTotCorr.at( key )->SetBinContent( i , 1. );
699  }
700 
701 
702  if( mDigiTotCorr.at( key )->GetBinContent( i ) < 0.05 || mDigiTotCorr.at( key )->GetBinContent( i ) > 10 ) {
703  mDigiTotCorr.at( key )->SetBinContent( i , 1. );
704  }
705  }
706  }
707  else{
708  if( isFileExisting( mFileNameCalibHistograms ) ) {
709  LOG_WARN << "unable to find histogram: " << hname << endm;
710  }
711 
712  for( size_t i=1; i<=2 * eTofConst::nStrips; i++ ) {
713  mDigiTotCorr.at( key )->SetBinContent( i , 1. );
714  }
715  }
716 
717  LOG_DEBUG << "histogram " << mDigiTotCorr.at( key )->GetName() << " filled" << endm;
718 
719 
720  //-------------------
721  // digi time corr
722  //-------------------
723  if( mDigiTimeCorr.count( key ) == 0 ) {
724  TString name = Form( "digiTimeCorr_%d", key );
725  mDigiTimeCorr[ key ] = new TH1F( name, name, 2 * eTofConst::nStrips, 0, 2 * eTofConst::nStrips );
726  }
727 
728  // (1) position offset
729  hname = Form( "calib_Sector%02d_ZPlane%d_Det%d_Pos_pfx", sector, zPlane, counter );
730  hProfile = ( TProfile* ) histFile->Get( hname );
731 
732  double dRunOffset = 0;
733  if (hPosOffsetProfile) {
734  int mCounterBin = hPosOffsetProfile->FindBin( 9*(sector -13) + 3 * (zPlane - 1) + counter );
735  dRunOffset = hPosOffsetProfile->GetBinContent( mCounterBin );
736  LOG_DEBUG << "setting run position offset to "<< dRunOffset<< " for counter "<< ( 9*(sector -13) + 3 * (zPlane - 1) + counter ) << endm;
737  }else{
738  LOG_INFO << "position offset histogram "<<hPosOffsetName <<" not found" << endm;
739  }
740 
741  if( hProfile ) {
742  for( size_t i=1; i<= eTofConst::nStrips; i++ ) {
743  mDigiTimeCorr.at( key )->AddBinContent( i , -1. * ( hProfile->GetBinContent( i ) + dRunOffset ) / mSignalVelocity.at( key ) );
744  mDigiTimeCorr.at( key )->AddBinContent( eTofConst::nStrips + i , ( hProfile->GetBinContent( i ) + dRunOffset ) / mSignalVelocity.at( key ) );
745  }
746  }
747  else{
748  if( isFileExisting( mFileNameCalibHistograms ) ) {
749  LOG_WARN << "unable to find histogram: " << hname << endm;
750  }
751  }
752 
753  // (2) T0 offset
754  hname = Form( "calib_Sector%02d_ZPlane%d_Det%d_T0corr", sector, zPlane, counter );
755  hProfile = ( TProfile* ) histFile->Get( hname );
756 
757  dRunOffset = 0;
758  if (hT0OffsetProfile) {
759  int mCounterBin = hT0OffsetProfile->FindBin( 9*(sector -13) + 3 * (zPlane - 1) + counter );
760  dRunOffset = hT0OffsetProfile->GetBinContent( mCounterBin );
761  }
762  LOG_DEBUG << "setting run time offset to "<< dRunOffset<< " for counter "<< ( 9*(sector -13) + 3 * (zPlane - 1) + counter ) << endm;
763 
764  if( hProfile && hProfile->GetNbinsX() == 1 ) {
765  LOG_DEBUG << "T0 histogram with 1 bin: " << key << endm;
766  for( size_t i=1; i<=2 * eTofConst::nStrips; i++ ) {
767  mDigiTimeCorr.at( key )->AddBinContent( i , hProfile->GetBinContent( 1 ) + dRunOffset );
768  }
769  }
770  else if( hProfile && hProfile->GetNbinsX() == eTofConst::nStrips ) {
771  LOG_DEBUG << "T0 histogram with 32 bins: " << key << endm;
772  for( size_t i=1; i<= eTofConst::nStrips; i++ ) {
773  mDigiTimeCorr.at( key )->AddBinContent( i , hProfile->GetBinContent( i ) + dRunOffset );
774  mDigiTimeCorr.at( key )->AddBinContent( i + eTofConst::nStrips , hProfile->GetBinContent( i ) + dRunOffset );
775  }
776  }
777  else{
778  if( isFileExisting( mFileNameCalibHistograms ) ) {
779  LOG_WARN << "unable to find histogram: " << hname << endm;
780  }
781  }
782 
783  LOG_DEBUG << "histogram " << mDigiTimeCorr.at( key )->GetName() << " filled" << endm;
784 
785 
786  //-------------------
787  // digi slewing corr
788  //-------------------
789  for( unsigned int chan = 0; chan < eTofConst::nChannelsPerCounter; chan++ ) {
790  unsigned int strip = ( chan % 32 ) + 1;
791  unsigned int side = ( chan / 32 ) + 1;
792  unsigned int key = sector * 100000 + zPlane * 10000 + counter * 1000 + strip * 10 + side;
793 
794  LOG_DEBUG << "channelId key: " << sector << " " << zPlane << " " << counter << " " << strip << " " << side << endm;
795 
796  hname = Form( "calib_Sector%02d_ZPlane%d_Det%d_Chan%d_Walk_pfx", sector, zPlane, counter, chan );
797  hProfile = ( TProfile* ) histFile->Get( hname );
798 
799  // check if channel-wise slewing parameters are available otherwise (due to low statistics) use detector-wise
800  if( !hProfile ) {
801  LOG_DEBUG << "unable to find histogram: " << hname << "--> check detector-wise" << endm;
802  hname = Form( "calib_Sector%02d_ZPlane%d_Det%d_AvWalk_pfx", sector, zPlane, counter );
803  hProfile = ( TProfile* ) histFile->Get( hname );
804  }
805 
806  if( hProfile ) {
807  unsigned int nbins = hProfile->GetNbinsX();
808 
809  if( mDigiSlewCorr.count( key ) == 0 ) {
810  // histogram could have variable bin size --> get vector of bin edges
811  std::vector< float > bins( nbins + 1 );
812 
813  for( size_t i=0; i<nbins; i++ ) {
814  bins.at( i ) = hProfile->GetXaxis()->GetBinLowEdge( i+1 );
815  }
816  bins.at( nbins ) = hProfile->GetXaxis()->GetBinUpEdge( nbins );
817 
818  TString name = Form( "digiSlewCorr_%d", key );
819  mDigiSlewCorr[ key ] = new TProfile( name, name, nbins, bins.data() );
820  }
821 
822  for( size_t iTotBin=1; iTotBin<=nbins; iTotBin++ ) {
823  mDigiSlewCorr.at( key )->Fill( hProfile->GetBinCenter( iTotBin ), hProfile->GetBinContent( iTotBin ), hProfile->GetBinEntries( iTotBin ) );
824  }
825  }
826  else{
827  LOG_DEBUG << "unable to find histogram: " << hname << endm;
828  }
829 
830  }
831  }
832  }
833  }
834 
835  for( auto& kv : mDigiTotCorr ) {
836  kv.second->SetDirectory( 0 );
837  }
838  for( auto& kv : mDigiTimeCorr ) {
839  kv.second->SetDirectory( 0 );
840  }
841  for( auto& kv : mDigiSlewCorr ) {
842  kv.second->SetDirectory( 0 );
843  }
844 
845  histFile->Close();
846  delete histFile;
847  histFile = nullptr;
848  }
849 
850  // --------------------------------------------------------------------------------------------
851 
852  // reset time correction
853  mResetTimeCorr = 0.;
854 
855  if( mFileNameResetTimeCorr.empty() ) {
856  LOG_INFO << "etofResetTimeCorr: no filename provided --> load database table" << endm;
857 
858  dbDataSet = GetDataBase( "Calibrations/etof/etofResetTimeCorr" );
859 
860  St_etofResetTimeCorr* etofResetTimeCorr = static_cast< St_etofResetTimeCorr* > ( dbDataSet->Find( "etofResetTimeCorr" ) );
861  if( !etofResetTimeCorr ) {
862  LOG_ERROR << "unable to get the reset time correction from the database" << endm;
863  return kStFatal;
864  }
865 
866  etofResetTimeCorr_st* resetTimeCorrTable = etofResetTimeCorr->GetTable();
867 
868  mResetTimeCorr = resetTimeCorrTable->resetTimeOffset;
869  }
870  else {
871  LOG_INFO << "etofResetTimeCorr: filename provided --> use parameter file: " << mFileNameResetTimeCorr.c_str() << endm;
872 
873  paramFile.open( mFileNameResetTimeCorr.c_str() );
874 
875  if( !paramFile.is_open() ) {
876  LOG_ERROR << "unable to get the 'etofResetTimeCorr' parameters from file --> file does not exist" << endm;
877  return kStFatal;
878  }
879 
880  std::map< unsigned int, float > param;
881  double temp;
882  double temp2 = 0;
883  while( paramFile >> temp ) {
884  paramFile >> temp2;
885  param[ ( unsigned int ) temp ] = temp2;
886  }
887 
888  paramFile.close();
889 
890  for( const auto& kv : param ) {
891  LOG_DEBUG << "run: " << kv.first << " --> reset time corr = " << kv.second << " ns" << endm;
892  }
893 
894  if( param.count( runnumber ) ) {
895  mResetTimeCorr = param.at( runnumber );
896  }
897  }
898 
899  LOG_INFO << "additional reset time offset correction: " << mResetTimeCorr << " ns" << endm;
900 
901  // --------------------------------------------------------------------------------------------
902 
903  // pulser peak tot
904  mPulserPeakTot.clear();
905 
906  if( mFileNamePulserTotPeak.empty() ) {
907  LOG_INFO << "etofPulserPeakTot: no filename provided --> load database table" << endm;
908 
909  dbDataSet = GetDataBase( "Calibrations/etof/etofPulserTotPeak" );
910 
911  St_etofPulserTotPeak* etofPulserTotPeak = static_cast< St_etofPulserTotPeak* > ( dbDataSet->Find( "etofPulserTotPeak" ) );
912  if( !etofPulserTotPeak ) {
913  LOG_ERROR << "unable to get the pulser tot peak parameters from the database" << endm;
914  return kStFatal;
915  }
916 
917  etofPulserTotPeak_st* pulserTotTable = etofPulserTotPeak->GetTable();
918 
919  for( size_t i=0; i<eTofConst::nCountersInSystem * 2; i++ ) {
920  if( pulserTotTable->pulserTot[ i ] > 0 ) {
921  mPulserPeakTot[ sideToKey( i ) ] = ( int ) pulserTotTable->pulserTot[ i ];
922  }
923  }
924  }
925  else {
926  LOG_INFO << "etofPulserPeakTot: filename provided --> use parameter file: " << mFileNamePulserTotPeak.c_str() << endm;
927 
928  paramFile.open( mFileNamePulserTotPeak.c_str() );
929 
930  if( !paramFile.is_open() ) {
931  LOG_ERROR << "unable to get the 'etofPulserTotPeak' parameters from file --> file does not exist" << endm;
932  return kStFatal;
933  }
934 
935  std::vector< float > param;
936  float temp;
937  while( paramFile >> temp ) {
938  param.push_back( temp );
939  }
940 
941  paramFile.close();
942 
943  if( param.size() != eTofConst::nCountersInSystem * 2 ) {
944  LOG_ERROR << "parameter file for 'etofPulserTotPeak' has not the right amount of entries: ";
945  LOG_ERROR << param.size() << " instead of " << eTofConst::nCountersInSystem * 2 << " !!!!" << endm;
946  return kStFatal;
947  }
948 
949  for( size_t i=0; i<eTofConst::nCountersInSystem * 2; i++ ) {
950  if( param.at( i ) > 0 ) {
951  mPulserPeakTot[ sideToKey( i ) ] = param.at( i );
952  }
953  }
954  }
955 
956  for( const auto& kv : mPulserPeakTot ) {
957  LOG_DEBUG << "side key: " << kv.first << " --> pulser peak tot = " << kv.second << " (bin)" << endm;
958  }
959 
960  // --------------------------------------------------------------------------------------------
961 
962  // pulser time difference (initialized to some useful value if pulser is not there for a whole run)
963  mPulserTimeDiffGbtx.clear();
964 
965  if( mFileNamePulserTimeDiffGbtx.empty() ) {
966  LOG_INFO << "etofPulserTimeDiffGbtx: no filename provided --> load database table" << endm;
967 
968  dbDataSet = GetDataBase( "Calibrations/etof/etofPulserTimeDiffGbtx" );
969 
970  St_etofPulserTimeDiffGbtx* etofPulserTimeDiffGbtx = static_cast< St_etofPulserTimeDiffGbtx* > ( dbDataSet->Find( "etofPulserTimeDiffGbtx" ) );
971  if( !etofPulserTimeDiffGbtx ) {
972  LOG_ERROR << "unable to get the pulser time difference parameters from the database" << endm;
973  return kStFatal;
974  }
975 
976  etofPulserTimeDiffGbtx_st* pulserTimeTable = etofPulserTimeDiffGbtx->GetTable();
977 
978  for( size_t i=0; i<eTofConst::nModules * eTofConst::nSides; i++ ) {
979  unsigned int sector = eTofConst::sectorStart + i / ( eTofConst::nPlanes * eTofConst::nSides );
980  unsigned int zPlane = eTofConst::zPlaneStart + ( i % ( eTofConst::nPlanes * eTofConst::nSides ) ) / eTofConst::nSides;
981  unsigned int side = eTofConst::counterStart + i % eTofConst::nSides;
982 
983  mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 1 * 10 + side ] = 0.;
984  mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 2 * 10 + side ] = pulserTimeTable->pulserTimeDiffGbtx[ i * 2 + 0 ];
985  mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 3 * 10 + side ] = pulserTimeTable->pulserTimeDiffGbtx[ i * 2 + 1 ];
986  }
987  }
988  else {
989  LOG_INFO << "etofPulserTimeDiff: filename provided --> use parameter file: " << mFileNamePulserTimeDiffGbtx.c_str() << endm;
990 
991  paramFile.open( mFileNamePulserTimeDiffGbtx.c_str() );
992 
993  if( !paramFile.is_open() ) {
994  LOG_ERROR << "unable to get the 'etotPulserTimeDiffGbtc' parameters from file --> file does not exist" << endm;
995  return kStFatal;
996  }
997 
998  std::vector< float > param;
999  float temp;
1000  while( paramFile >> temp ) {
1001  param.push_back( temp );
1002  }
1003 
1004  paramFile.close();
1005 
1006  if( param.size() != eTofConst::nModules * eTofConst::nSides * 2 ) {
1007  LOG_ERROR << "parameter file for 'etofPulserTimeDiffGbtx' has not the right amount of entries: ";
1008  LOG_ERROR << param.size() << " instead of " << eTofConst::nModules * eTofConst::nSides << " !!!!" << endm;
1009  return kStFatal;
1010  }
1011 
1012  for( size_t i=0; i<eTofConst::nModules * eTofConst::nSides; i++ ) {
1013  unsigned int sector = eTofConst::sectorStart + i / ( eTofConst::nPlanes * eTofConst::nSides );
1014  unsigned int zPlane = eTofConst::zPlaneStart + ( i % ( eTofConst::nPlanes * eTofConst::nSides ) ) / eTofConst::nSides;
1015  unsigned int side = eTofConst::counterStart + i % eTofConst::nSides;
1016 
1017  mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 1 * 10 + side ] = 0.;
1018  mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 2 * 10 + side ] = param.at( i * 2 + 0 );
1019  mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 3 * 10 + side ] = param.at( i * 2 + 1 );
1020  }
1021  }
1022 
1023 
1024  // check if all values in the map are zero --> if yes, disable pulser averaging inside a Gbtx
1025  bool allZero = true;
1026 
1027  for( const auto& kv : mPulserTimeDiffGbtx ) {
1028  if( fabs( kv.second ) > 1e-4 ) allZero = false;
1029 
1030  if( mDebug ) {
1031  LOG_INFO << "side key: " << kv.first << " --> pulser time diff inside Gbtx ( to counter 1 ) = " << kv.second << endm;
1032  }
1033  }
1034 
1035  if( allZero ) {
1036  mUsePulserGbtxDiff = false;
1037  LOG_INFO << "the use of pulser relations inside a Gbtx is turned off" << endm;
1038  }
1039 
1040  //reset histo keeping track of mod average distance to clock
1041  mHistograms.at( "pulserDigiTimeDiff_GbtxCorrProfMod" )->Reset();
1042 
1043  // --------------------------------------------------------------------------------------------
1044 
1045 
1046  return kStOk;
1047 }
1048 
1049 
1050 //_____________________________________________________________
1051 Int_t
1052 StETofCalibMaker::FinishRun( Int_t runnumber )
1053 {
1054  if( mHwMap ) {
1055  delete mHwMap;
1056  mHwMap = nullptr;
1057  }
1058 
1059  // delete histograms from the map
1060  for( auto kv = mDigiTotCorr.begin(); kv != mDigiTotCorr.end(); kv++ ) {
1061  delete kv->second;
1062  kv->second = nullptr;
1063  }
1064  mDigiTotCorr.clear();
1065 
1066  for( auto kv = mDigiTimeCorr.begin(); kv != mDigiTimeCorr.end(); kv++ ) {
1067  delete kv->second;
1068  kv->second = nullptr;
1069  }
1070  mDigiTimeCorr.clear();
1071 
1072  for( auto kv = mDigiSlewCorr.begin(); kv != mDigiSlewCorr.end(); kv++ ) {
1073  delete kv->second;
1074  kv->second = nullptr;
1075  }
1076  mDigiSlewCorr.clear();
1077 
1078  mPulserPeakTot.clear();
1079  mPulserTimeDiff.clear();
1080 
1081  mJumpingPulsers.clear();
1082 
1083  return kStOk;
1084 }
1085 
1086 
1087 //_____________________________________________________________
1088 Int_t
1090 {
1091  if( mDoQA ) {
1092  LOG_INFO << "Finish() - writing *.etofCalib.root ..." << endm;
1093  setHistFileName();
1094  writeHistograms();
1095  }
1096 
1097  return kStOk;
1098 }
1099 
1100 //_____________________________________________________________
1101 Int_t
1103 {
1104  LOG_DEBUG << "StETofCalibMaker::Make(): starting ..." << endm;
1105 
1106  mEvent = ( StEvent* ) GetInputDS( "StEvent" );
1107  //mEvent = NULL; //don't check for StEvent for genDst.C testing. PW
1108 
1109  if ( mEvent ) {
1110  LOG_DEBUG << "Make(): running on StEvent" << endm;
1111 
1112  StETofCollection* etofCollection = mEvent->etofCollection();
1113 
1114  if( !etofCollection ) { //additional check for empty StEvents structures produced by other Makers. Needed for genDst.C
1115  LOG_WARN << "Make() - Found StEvent data structure, but no eTOF collection. Try MuDst processing instead" << endm;
1116  mMuDst = ( StMuDst* ) GetInputDS( "MuDst" );
1117 
1118  if( mMuDst ) {
1119  LOG_DEBUG << "Make() - running on MuDsts" << endm;
1120 
1121  processMuDst();
1122 
1123  return kStOk;
1124  }
1125  }
1126 
1127  processStEvent();
1128 
1129  return kStOk;
1130  }
1131  else {
1132  LOG_DEBUG << "Make(): no StEvent found" << endm;
1133 
1134  mMuDst = ( StMuDst* ) GetInputDS( "MuDst" );
1135 
1136  if( mMuDst ) {
1137  LOG_DEBUG << "Make(): running on MuDsts" << endm;
1138 
1139  processMuDst();
1140 
1141  return kStOk;
1142  }
1143  else {
1144  LOG_WARN << "Make() - no StMuDst or StEvent" << endm;
1145  return kStOk;
1146  }
1147  }
1148 
1149 }
1150 
1151 
1152 //_____________________________________________________________
1153 void
1155 {
1156  StETofCollection* etofCollection = mEvent->etofCollection();
1157 
1158  if( !etofCollection ) {
1159  LOG_WARN << "processStEvent() - no etof collection" << endm;
1160  return;
1161  }
1162 
1163  if( !etofCollection->digisPresent() ) {
1164  LOG_WARN << "processStEvent() - no digis present" << endm;
1165  return;
1166  }
1167 
1168  if( !etofCollection->etofHeader() ) {
1169  LOG_WARN << "processStEvent() - no header" << endm;
1170  return;
1171  }
1172 
1173  //---------------------------------
1174 
1175  StETofHeader* etofHeader = etofCollection->etofHeader();
1176  StSPtrVecETofDigi& etofDigis = etofCollection->etofDigis();
1177 
1178  /*if( mDoQA ){
1179  LOG_INFO << "filling missmatch histograms now" << endm;
1180  TClass* headerClass = etofHeader->IsA();
1181  if( headerClass->GetClassVersion() > 1 ){
1182  LOG_INFO << "getting missmatch vector" << endm;
1183  std::vector< Bool_t > vMissmatchVec = etofHeader->missMatchFlagVec(); //lookup error?
1184  int iGet4Id = 0;
1185  for( auto iMissMatchFlag : vMissmatchVec ){
1186  int iCounter = iGet4Id % 16; //probalby wrong!
1187  mHistograms[ "ETOF_QA_daqMissmatches_get4" ]->Fill(iGet4Id, iMissMatchFlag);
1188  mHistograms[ "ETOF_QA_daqMissmatches_counter" ]->Fill(iCounter, iMissMatchFlag);
1189  }
1190  }
1191  }*/
1192 
1193  size_t nDigis = etofDigis.size();
1194  if( mDebug ) {
1195  LOG_INFO << "processStEvent() - # fired eTOF digis : " << nDigis << endm;
1196  }
1197 
1198  mTriggerTime = triggerTime( etofHeader );
1199  mResetTime = fmod( resetTime( etofHeader ), eTofConst::bTofClockCycle );
1200 
1201  std::map< unsigned int, std::vector< unsigned int > > pulserCandMap;
1202 
1204  for( size_t i=0; i<nDigis; i++ ) {
1205  // LOG_INFO << "Digi array" << endm;
1206  StETofDigi* aDigi = etofDigis[ i ];
1207 
1208  if( !aDigi ) {
1209  LOG_WARN << "No digi found" << endm;
1210  continue;
1211  }
1212  // LOG_INFO << "Digi reset" << endm;
1214  resetToRaw( aDigi );
1215 
1216  // LOG_INFO << "Mapping reset" << endm;
1219  applyMapping( aDigi );
1220 
1222  if( mRunYear != 2018 ) {
1223  flagPulserDigis( aDigi, i, pulserCandMap );
1224  }
1225  }
1226 
1227  // LOG_INFO << "pulsers" << endm;
1228  if( mDebug ) {
1229  LOG_INFO << "size of pulserCandMap: " << pulserCandMap.size() << endm;
1230  }
1231  calculatePulserOffsets( pulserCandMap );
1232 
1233  // collect status bit information and fill good event flag for 2020+ data
1234  TClass* headerClass = etofHeader->IsA();
1235  if( headerClass->GetClassVersion() > 2 ){
1236 
1237  std::vector<bool> goodEventFlagVec;
1238  std::vector<bool> hasPulsersVec;
1239 
1240  //drag along pulser information
1241  for( unsigned int iCounter = 0; iCounter < 108; iCounter++){
1242  if ( !(mNPulsersCounter.count(iCounter) ) ){
1243  hasPulsersVec.push_back(false);
1244  }else{
1245  hasPulsersVec.push_back(mNPulsersCounter[iCounter] == 2);
1246  }
1247  }
1248  if (hasPulsersVec.size() == 108){
1249  //etofHeader->setHasPulsersVec(hasPulsersVec); // not working but not of relevance at the moment
1250  }
1251 
1252  //fill good event flag into header
1253  for( unsigned int iGet4 = 0; iGet4 < 1728; iGet4++){
1254  goodEventFlagVec.push_back(!etofHeader->missMatchFlagVec().at(iGet4));
1255  }
1256 
1257  if (goodEventFlagVec.size() == 1728){
1258  etofHeader->setGoodEventFlagVec(goodEventFlagVec);
1259  }
1260  }
1261 
1262 
1264  StructStuckFwDigi current = { -1, -1., -1. };
1265  StructStuckFwDigi prev = { -1, -1., -1. };
1266  int nDuplicates = 0;
1267 
1268  for( size_t i=0; i<nDigis; i++ ) {
1269  StETofDigi* aDigi = etofDigis[ i ];
1270 
1271  if( !aDigi ) {
1272  LOG_WARN << "No digi found" << endm;
1273  continue;
1274  }
1275 
1276  current.geomId = aDigi->sector() * 100000 + aDigi->zPlane() * 10000 + aDigi->counter() * 1000 + aDigi->strip() * 10 + aDigi->side();
1277  current.tot = aDigi->rawTot();
1278  current.time = aDigi->rawTime();
1279 
1280  // ignore digis that were sent in bulk from the same channel with exactly the same tot and time due to stuck firmware
1281  auto it = std::find( mStuckFwDigi.begin(), mStuckFwDigi.end(), current );
1282  if( it != mStuckFwDigi.end() ) {
1283  if( mDebug ) {
1284  LOG_INFO << "digi from stuck firmware (s" << aDigi->sector() << " m" << aDigi->zPlane() << " c" << aDigi->counter() << ") --> ignore" << endm;
1285  }
1286 
1287  nDuplicates++;
1288  continue;
1289  }
1290  else if( current == prev ) {
1291  mStuckFwDigi.push_back( current );
1292  resetToRaw( mMuDst->etofDigi( i-1 ) );
1293 
1294  nDuplicates++;
1295  continue;
1296  }
1297  else {
1298  prev = current;
1299  }
1300 
1301 
1304  applyCalibration( aDigi, etofHeader );
1305  }
1306 
1307  if( mDoQA && nDuplicates > 0 ) {
1308  LOG_INFO << "*** # duplicate digis from stuck firmware: " << nDuplicates << endm;
1309  for( const auto& v : mStuckFwDigi ) {
1310  LOG_INFO << "stuck channel with geomId: " << v.geomId << " " << v.tot << " " << v.time << endm;
1311  }
1312  }
1313  mStuckFwDigi.clear();
1314 }
1315 
1316 
1317 //_____________________________________________________________
1318 void
1320 {
1321  LOG_DEBUG << "processMuDst(): starting ..." << endm;
1322 
1323 
1324  if( !mMuDst->etofArray( muETofDigi ) ) {
1325  LOG_WARN << "processMuDst() - no digi array" << endm;
1326  return;
1327  }
1328 
1329  if( !mMuDst->numberOfETofDigi() ) {
1330  LOG_WARN << "processMuDst() - no digis present" << endm;
1331  return;
1332  }
1333 
1334  if( !mMuDst->etofArray( muETofHeader ) ) {
1335  LOG_WARN << "processMuDst() - no digi array" << endm;
1336  return;
1337  }
1338 
1339  StMuETofHeader* etofHeader = mMuDst->etofHeader();
1340  mNPulsersCounter.clear();
1341 
1342  //---------------------------------
1343 
1344 /* if( mDoQA ){
1345  LOG_INFO << "filling missmatch histograms now" << endm;
1346  TClass* headerClass = etofHeader->IsA();
1347  if( headerClass->GetClassVersion() > 1 ){
1348  LOG_INFO << "getting missmatch vector" << endm;
1349  std::vector< Bool_t > vMissmatchVec = etofHeader->missMatchFlagVec(); //lookup error?
1350  int iGet4Id = 0;
1351  for( auto iMissMatchFlag : vMissmatchVec ){
1352  int iCounter = iGet4Id % 16; //probalby wrong!
1353  mHistograms[ "ETOF_QA_daqMissmatches_get4" ]->Fill(iGet4Id, iMissMatchFlag);
1354  mHistograms[ "ETOF_QA_daqMissmatches_counter" ]->Fill(iCounter, iMissMatchFlag);
1355  }
1356  }
1357  }*/
1358 
1359  size_t nDigis = mMuDst->numberOfETofDigi();
1360  //LOG_INFO << "processMuDst() - # fired eTOF digis : " << nDigis << endm;
1361 
1362  mTriggerTime = triggerTime( ( StETofHeader* ) etofHeader );
1363  mResetTime = fmod( resetTime( ( StETofHeader* ) etofHeader ), eTofConst::bTofClockCycle );
1364  std::map< unsigned int, std::vector< unsigned int >> pulserCandMap;
1365 
1366 
1368  for( size_t i=0; i<nDigis; i++ ) {
1369  //LOG_INFO << "accessing etof digis: "<< i <<"/"<< nDigis << endm;
1370  StMuETofDigi* aDigi = mMuDst->etofDigi( i );
1371 
1372  if( !aDigi ) {
1373  LOG_WARN << "No digi found" << endm;
1374  continue;
1375  }
1377  //LOG_INFO << "resetting digi "<< i <<"/"<< nDigis << endm;
1378  resetToRaw( aDigi );
1379 
1380 
1383  //LOG_INFO << "mapping digi: "<< i <<"/"<< nDigis << endm;
1384  applyMapping( aDigi );
1385 
1386 
1387 
1389  //LOG_INFO << "pulser digi flagging: "<< i <<"/"<< nDigis << endm;
1390  if( mRunYear != 2018 ) {
1391  flagPulserDigis( aDigi, i, pulserCandMap );
1392  }
1393  }
1394 
1395 
1396  //LOG_INFO << "size of pulserCandMap: " << pulserCandMap.size() << endm;
1397 
1398  calculatePulserOffsets( pulserCandMap );
1399 
1400 
1401  // collect status bit information and fill good event flag for 2020+ data
1402  TClass* headerClass = etofHeader->IsA();
1403  if( headerClass->GetClassVersion() > 2 ){
1404 
1405  std::vector<bool> goodEventFlagVec;
1406  std::vector<bool> hasPulsersVec;//
1407 
1408  //drag along pulser information
1409  for( unsigned int iCounter = 0; iCounter < 108; iCounter++){
1410  if ( !(mNPulsersCounter.count(iCounter) ) ){
1411  hasPulsersVec.push_back(false);
1412  }else{
1413  hasPulsersVec.push_back(mNPulsersCounter[iCounter] == 2);
1414  }
1415  }
1416 
1417  if (hasPulsersVec.size() == 108){
1418  etofHeader->setHasPulsersVec(hasPulsersVec);
1419  }
1420 
1421  //fill good event flag into header
1422  for( unsigned int iGet4 = 0; iGet4 < 1728; iGet4++){
1423  goodEventFlagVec.push_back(!etofHeader->missMatchFlagVec().at(iGet4));
1424  }
1425 
1426  if (goodEventFlagVec.size() == 1728){
1427  etofHeader->setGoodEventFlagVec(goodEventFlagVec);
1428  }
1429  }
1430 
1432  StructStuckFwDigi current = { -1, -1., -1. };
1433  StructStuckFwDigi prev = { -1, -1., -1. };
1434  int nDuplicates = 0;
1435 
1436  for( size_t i=0; i<nDigis; i++ ) {
1437  StMuETofDigi* aDigi = mMuDst->etofDigi( i );
1438 
1439  if( !aDigi ) {
1440  LOG_WARN << "No digi found" << endm;
1441  continue;
1442  }
1443 
1444  current.geomId = aDigi->sector() * 100000 + aDigi->zPlane() * 10000 + aDigi->counter() * 1000 + aDigi->strip() * 10 + aDigi->side();
1445  current.tot = aDigi->rawTot();
1446  current.time = aDigi->rawTime();
1447 
1448  // ignore digis that were sent in bulk from the same channel with exactly the same tot and time due to stuck firmware
1449  auto it = std::find( mStuckFwDigi.begin(), mStuckFwDigi.end(), current );
1450  if( it != mStuckFwDigi.end() ) {
1451  if( mDebug ) {
1452  LOG_INFO << "digi from stuck firmware (s" << aDigi->sector() << " m" << aDigi->zPlane() << " c" << aDigi->counter() << ") --> ignore" << endm;
1453  }
1454 
1455  nDuplicates++;
1456  continue;
1457  }
1458  else if( current == prev ) {
1459  mStuckFwDigi.push_back( current );
1460  resetToRaw( mMuDst->etofDigi( i-1 ) );
1461 
1462  nDuplicates++;
1463  continue;
1464  }
1465  else {
1466  prev = current;
1467  }
1468 
1471  applyCalibration( aDigi, etofHeader );
1472  }
1473 
1474  if( mDoQA && nDuplicates > 0 ) {
1475  LOG_INFO << "*** # duplicate digis from stuck firmware: " << nDuplicates << endm;
1476  for( const auto& v : mStuckFwDigi ) {
1477  LOG_INFO << "stuck channel with geomId: " << v.geomId << " " << v.tot << " " << v.time << endm;
1478  }
1479  }
1480  mStuckFwDigi.clear();
1481 }
1482 //_____________________________________________________________
1483 
1484 
1485 //_____________________________________________________________
1486 bool
1487 StETofCalibMaker::isFileExisting( const std::string fileName )
1488 {
1489  std::ifstream infile( fileName );
1490  return infile.good();
1491 }
1492 
1493 
1494 //_____________________________________________________________
1498 void
1499 StETofCalibMaker::resetToRaw( StETofDigi* aDigi )
1500 {
1501  if( !aDigi ) return;
1502 
1503  aDigi->setGeoAddress( 0, 0, 0, 0, 0 );
1504  aDigi->setCalibTime( 0. );
1505  aDigi->setCalibTot( -1. );
1506 
1507  aDigi->setAssociatedHit( nullptr );
1508 }
1509 
1510 
1511 //_____________________________________________________________
1515 void
1516 StETofCalibMaker::applyMapping( StETofDigi* aDigi )
1517 {
1518  std::vector< unsigned int > geomVec;
1519 
1520  unsigned int rocId = aDigi->rocId();
1521  unsigned int get4Id = aDigi->get4Id();
1522  unsigned int elChan = aDigi->elChan();
1523 
1524  mHwMap->mapToGeom( rocId, get4Id, elChan, geomVec );
1525 
1526  if( geomVec.size() < 5 ) {
1527  if( mDebug ) {
1528  LOG_ERROR << "geometry vector has wrong size !!! --> skip digi" << endm;
1529  }
1530  return;
1531  }
1532 
1533  unsigned int sector = geomVec.at( 0 );
1534  unsigned int zplane = geomVec.at( 1 );
1535  unsigned int counter = geomVec.at( 2 );
1536  unsigned int strip = geomVec.at( 3 );
1537  unsigned int side = geomVec.at( 4 );
1538 
1539  if( mDebug && ( sector == 0 || zplane == 0 || counter == 0 || strip == 0 || side == 0 ) ) {
1540  LOG_ERROR << "geometry vector has entries equal to zero !!! --> skip digi" << endm;
1541  }
1542 
1543  aDigi->setGeoAddress( sector, zplane, counter, strip, side );
1544 
1545  if( mDebug ) {
1546  // print out the new information
1547  LOG_DEBUG << "sector, zplane, counter, strip, side: " << aDigi->sector() << ", ";
1548  LOG_DEBUG << aDigi->zPlane() << ", " << aDigi->counter() << ", ";
1549  LOG_DEBUG << aDigi->strip() << ", " << aDigi->side() << endm;
1550 
1551  LOG_DEBUG << "continuous module number: " << mHwMap->module( aDigi->sector(), aDigi->zPlane() ) << endm;
1552  }
1553 }
1554 
1555 
1556 //_____________________________________________________________
1560 void
1561 StETofCalibMaker::flagPulserDigis( StETofDigi* aDigi, unsigned int index, std::map< unsigned int, std::vector< unsigned int > >& pulserDigiMap )
1562 {
1563  bool isPulserCand = false;
1564 
1565  unsigned int key = aDigi->sector() * 1000 + aDigi->zPlane() * 100 + aDigi->counter() * 10 + aDigi->side();
1566 
1567 
1568  // pulser channel
1569  if( ( aDigi->strip() == 1 && aDigi->side() == 1 ) || ( aDigi->strip() == 32 && aDigi->side() == 2 ) ) {
1570  float timeToTrigger = aDigi->rawTime() - mTriggerTime;
1571 
1572 
1573  float totToPeak = aDigi->rawTot() - mPulserPeakTot.at( key );
1574  float totToHalfPeak = aDigi->rawTot() - mPulserPeakTot.at( key ) * 0.5;
1575 
1576 
1577  if( timeToTrigger > mPulserWindow.at( aDigi->rocId() ).first && timeToTrigger < mPulserWindow.at( aDigi->rocId() ).second ) {
1578  if( fabs( totToPeak ) < 25 || fabs( totToHalfPeak ) < 10 ) {
1579  isPulserCand = true;
1580  }
1581  }
1582  }
1583 
1584 
1585  if( isPulserCand ) {
1586  pulserDigiMap[ key ].push_back( index );
1587  }
1588 
1589 }
1590 
1591 
1592 //_____________________________________________________________
1598 void
1599 StETofCalibMaker::calculatePulserOffsets( std::map< unsigned int, std::vector< unsigned int > >& pulserDigiMap )
1600 {
1601  if( mDebug ) {
1602  for( auto it=pulserDigiMap.begin(); it!=pulserDigiMap.end(); it++ ) {
1603  LOG_DEBUG << "channel: " << it->first << " nCandidates: " << it->second.size() << endm;
1604  }
1605  }
1606 
1607  if( mReferencePulserIndex == 0 ) {
1608  if( mDebug ) {
1609  LOG_INFO << "reference pulser index is 0 --> pulser correction is turned off" << endm;
1610  }
1611  return;
1612  }
1613 
1614  if( mDebug ) {
1615  LOG_INFO << "reference pulser index: " << mReferencePulserIndex << endm;
1616  }
1617 
1618  std::map< int, double > pulserTimes;
1619  mNPulsersCounter.clear();
1620  mPulserPresent.clear(); //clear map of present pulsers in each event
1621 
1622  // loop over all candidates to find real pulser, save time in pulserTimes map
1623  for( auto it=pulserDigiMap.begin(); it!=pulserDigiMap.end(); it++ ) {
1624  if( it->second.size() == 0 ) {
1625  continue;
1626  }
1627  int sideIndex = it->first;
1628 
1629  double bestDiff = 100000;
1630  int candIndex = -1;
1631 
1632  for( size_t j=0; j<it->second.size(); j++ ) {
1633  double pulserTime = 0.;
1634  double pulserTot = 0.;
1635  if( mMuDst ) {
1636  pulserTime = mMuDst->etofDigi( it->second.at( j ) )->rawTime();
1637  pulserTot = mMuDst->etofDigi( it->second.at( j ) )->rawTot();
1638  } else if( mEvent ) {
1639  pulserTime = ( mEvent->etofCollection()->etofDigis() )[ it->second.at( j ) ]->rawTime();
1640  pulserTot = ( mEvent->etofCollection()->etofDigis() )[ it->second.at( j ) ]->rawTot();
1641  }
1642  double timeToTrigger = pulserTime - mTriggerTime;
1643  double totToPeak = pulserTot - mPulserPeakTot.at( sideIndex );
1644 
1645  if( mDebug && it->second.size() > 1 ) {
1646  LOG_INFO << it->second.size() << " pulsers @ " << sideIndex << " : timeToTrigger: " << timeToTrigger << " tot: " << pulserTot << endm;
1647  }
1648 
1649  // find "best fitting digi", remove other digis (likely misidentified noise)
1650  double currentDiff = fabs( timeToTrigger - mPulserPeakTime ) * 0.1 + fabs( totToPeak ); // might need better criterion? Normalisation to widths? PW
1651  if( currentDiff < bestDiff ) {
1652  bestDiff = currentDiff;
1653  candIndex = j;
1654  }
1655  }
1656 
1657  if( mDebug && it->second.size() > 1 ) {
1658  LOG_INFO << " --> selected CAND-INDEX: " << candIndex << endm;
1659  }
1660 
1661  double pulserTime = 0.;
1662 
1663  if( mMuDst ) {
1664  pulserTime = mMuDst->etofDigi( it->second.at( candIndex ) )->rawTime();
1665 
1666  // set calibTot to -999. to exclude it from being calibrated in the next step --> pulser will not be used to build hits
1667  mMuDst->etofDigi( it->second.at( candIndex ) )->setCalibTot( -999. );
1668  } else if( mEvent ) {
1669  pulserTime = ( mEvent->etofCollection()->etofDigis() )[ it->second.at( candIndex ) ]->rawTime();
1670 
1671  // set calibTot to -999. to exclude it from being calibrated in the next step --> pulser will not be used to build hits
1672  mEvent->etofCollection()->etofDigis() [ it->second.at( candIndex ) ]->setCalibTot( -999. );
1673  }
1674 
1675  pulserTimes[ sideIndex ] = pulserTime;
1676 
1677  int sector = sideIndex / 1000;
1678  int plane = ( sideIndex % 1000) / 100;
1679  int counter = ( sideIndex % 100) / 10;
1680  int key = 9 * ( sector - 13 ) + 3 * ( plane - 1 ) + ( counter - 1 );
1681  if( mNPulsersCounter.count( key ) ){
1682  mNPulsersCounter[key]++;
1683  }else{
1684  mNPulsersCounter[key] = 1;
1685  }
1686 
1687  mPulserPresent[ sideIndex ] = true;
1688  }
1689 
1690  double referenceTime = 0.;
1691 
1692  // update reference time (for QA)
1693  if( pulserTimes.count( mReferencePulserIndex ) ) {
1694  referenceTime = pulserTimes.at( mReferencePulserIndex ); //only updated for QA?? needed to remove smeared pulsers
1695  if( mDoQA ) {
1696  if( mDebug ) {
1697  LOG_INFO << "preliminary reference time:" << referenceTime << endm;
1698  }
1699  }
1700  }
1701 
1702 
1703  // deal with the pulser times --> tweek pulser times based on time differences inside/outside a Gbtx
1704  if( mUsePulserGbtxDiff ) {
1705  for( int gbtxId = 0; gbtxId < eTofConst::nModules * eTofConst::nSides; gbtxId++ ) {
1706  int sector = eTofConst::sectorStart + gbtxId / ( eTofConst::nPlanes * eTofConst::nSides );
1707  int zPlane = eTofConst::zPlaneStart + ( gbtxId % ( eTofConst::nPlanes * eTofConst::nSides ) ) / eTofConst::nSides;
1708  int side = eTofConst::sideStart + gbtxId % eTofConst::nSides;
1709 
1710  int partialKey = sector * 1000 + zPlane * 100 + side;
1711 
1712  vector< double > times( eTofConst::nCounters );
1713 
1714  double average = 0.;
1715  int nAv = 0;
1716 
1717  for( int counter = 1; counter <= eTofConst::nCounters; counter++ ) {
1718  int key = partialKey + 10 * counter;
1719 
1720  if( pulserTimes.count( key ) ) {
1721  if( mDoQA ) {// fill if all relevant pulsers are available.
1722  if( pulserTimes.count( partialKey + 10 ) ) {
1723  mHistograms.at( "pulserDigiTimeToCounter1" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, pulserTimes.at( partialKey + 10 ) - pulserTimes.at( key ) );
1724  }
1725  if( pulserTimes.count( partialKey + 20 ) ) {
1726  mHistograms.at( "pulserDigiTimeToCounter2" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, pulserTimes.at( partialKey + 20 ) - pulserTimes.at( key ) );
1727  }
1728  if( pulserTimes.count( partialKey + 30 ) ) {
1729  mHistograms.at( "pulserDigiTimeToCounter3" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, pulserTimes.at( partialKey + 30 ) - pulserTimes.at( key ) );
1730  }
1731  }
1732 
1733  times.at( counter - 1 ) = pulserTimes.at( key ) + mPulserTimeDiffGbtx.at( key ); //substract pulser time difference from database table
1734 
1735  bool isNonSmearedPulser = false;
1736  if( referenceTime != 0 ) {
1737  double dist = times.at( counter - 1 ) - referenceTime; //distance to reference pulser
1738  double redDist = mHistograms.at( "pulserDigiTimeDiff_GbtxCorrProfMod" )->GetBinContent( gbtxId * eTofConst::nCounters + counter ); // average distance to next clock edge for this pulser
1739 
1740  double modDist = fmod( dist - redDist, eTofConst::coarseClockCycle ); //Distance to "normal" offset. full clock cycle distances are thrown out? Why?
1741  modDist = modDist - eTofConst::coarseClockCycle * round( modDist / eTofConst::coarseClockCycle ); //substract a clock cycle if modDist > 0.5 coarseClockCycle
1742  //=> -0.5*coarseClockCycle < modDist < 0.5*coarseClockCycle => Distance to closest clock cycle
1743 
1744  if( redDist == 0 || fabs( modDist ) > 0.5 ) { //> 0.5ns + n*ClockCycle away from reference pulser. Hard cut? If the first pulser is off, all following will be neglected!
1745  if( redDist != 0 ) LOG_INFO << "too far away --> smeared pulser: " << key << "(" << gbtxId << "-" << counter << ")" << endm;
1746  redDist = dist; //empty in the beginning, Example distance to reference pulser
1747  }
1748  else {
1749  redDist += modDist; //adds always up?
1750  isNonSmearedPulser = true;
1751  }
1752 
1753  mHistograms.at( "pulserDigiTimeDiff_GbtxCorrProf" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, dist );
1754  mHistograms.at( "pulserDigiTimeDiff_GbtxCorrProfMod" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, redDist ); // TProfile! => Average!
1755 
1756  if( mDoQA ) {
1757  mHistograms.at( "pulserDigiTimeDiff_GbtxCorr" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, dist ); //Pulser offset on GBTX from database substracted
1758  mHistograms.at( "pulserDigiTimeDiff" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, pulserTimes.at( key ) - referenceTime ); //Pulser offset on GBTX not substracted
1759  }
1760  }
1761 
1762  // only use non-smeared pulsers for the
1763  if( isNonSmearedPulser ) {
1764  average += times.at( counter - 1 );
1765  nAv++;
1766  }
1767  else {
1768  times.at( counter - 1 ) = 0.;
1769  }
1770 
1771  }
1772  }
1773 
1774 
1775  if( nAv == eTofConst::nCounters ) { //all pulser present, check for single pulser jumps by comparing to average of the other two.
1776  double diff12 = fabs( times.at( 0 ) - times.at( 1 ) );
1777  double diff13 = fabs( times.at( 0 ) - times.at( 2 ) );
1778  double diff23 = fabs( times.at( 1 ) - times.at( 2 ) );
1779 
1780  if( diff12 > 0.2 && diff13 > 0.2 && diff23 < 0.2 ) {
1781  average -= times.at( 0 );
1782  nAv--;
1783 
1784  if( fabs( times.at( 0 ) - ( average / nAv ) + eTofConst::coarseClockCycle ) < 0.2 ) {
1785  mJumpingPulsers[ partialKey + 10 ] = 1;
1786 
1787  times.at( 0 ) += eTofConst::coarseClockCycle;
1788  average += times.at( 0 );
1789  nAv++;
1790  }
1791  else if( fabs( times.at( 0 ) - ( average / nAv ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1792  mJumpingPulsers[ partialKey + 10 ] = -1;
1793 
1794  times.at( 0 ) -= eTofConst::coarseClockCycle;
1795  average += times.at( 0 );
1796  nAv++;
1797  }
1798 
1799  if( mDoQA ) {
1800  mHistograms.at( "1_off" )->Fill( gbtxId * eTofConst::nCounters + 1.5, times.at( 0 ) - ( average / nAv ) );
1801  }
1802  }
1803  else if( diff12 > 0.2 && diff13 < 0.2 && diff23 > 0.2 ) {
1804  average -= times.at( 1 );
1805  nAv--;
1806 
1807  if( fabs( times.at( 1 ) - ( average / nAv ) + eTofConst::coarseClockCycle ) < 0.2 ) {
1808  mJumpingPulsers[ partialKey + 20 ] = 1;
1809 
1810  times.at( 1 ) += eTofConst::coarseClockCycle;
1811  average += times.at( 1 );
1812  nAv++;
1813  }
1814  else if( fabs( times.at( 1 ) - ( average / nAv ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1815  mJumpingPulsers[ partialKey + 20 ] = -1;
1816 
1817  times.at( 1 ) -= eTofConst::coarseClockCycle;
1818  average += times.at( 1 );
1819  nAv++;
1820  }
1821 
1822  if( mDoQA ) {
1823  mHistograms.at( "2_off" )->Fill( gbtxId * eTofConst::nCounters + 1.5, times.at( 1 ) - ( average / nAv ) );
1824  }
1825  }
1826  else if( diff12 < 0.2 && diff13 > 0.2 && diff23 > 0.2 ) {
1827  average -= times.at( 2 );
1828  nAv--;
1829 
1830  if( fabs( times.at( 2 ) - ( average / nAv ) + eTofConst::coarseClockCycle ) < 0.2 ) {
1831  mJumpingPulsers[ partialKey + 30 ] = 1;
1832 
1833  times.at( 2 ) += eTofConst::coarseClockCycle;
1834  average += times.at( 2 );
1835  nAv++;
1836  }
1837  else if( fabs( times.at( 2 ) - ( average / nAv ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1838  mJumpingPulsers[ partialKey + 30 ] = -1;
1839 
1840  times.at( 2 ) -= eTofConst::coarseClockCycle;
1841  average += times.at( 2 );
1842  nAv++;
1843  }
1844 
1845  if( mDoQA ) {
1846  mHistograms.at( "3_off" )->Fill( gbtxId * eTofConst::nCounters + 1.5, times.at( 2 ) - ( average / nAv ) );
1847  }
1848  }
1849  }
1850 
1851  if( nAv == eTofConst::nCounters - 1 ) {
1852  // if there are two pulsers, restore missing pulser from average of the other two
1853  if( times.at( 0 ) > 0 && times.at( 1 ) > 0 && fabs( fabs( times.at( 0 ) - times.at( 1 ) ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1854  if( mJumpingPulsers.count( partialKey + 10 ) ) {
1855  //LOG_INFO << gbtxId << " ### case 1 (1) ### " << endm;
1856  times.at( 0 ) += mJumpingPulsers.at( partialKey + 10 ) * eTofConst::coarseClockCycle;
1857  average += mJumpingPulsers.at( partialKey + 10 ) * eTofConst::coarseClockCycle;
1858  }
1859  else if( mJumpingPulsers.count( partialKey + 20 ) ) {
1860  //LOG_INFO << gbtxId << " ### case 1 (2) ### " << endm;
1861  times.at( 1 ) += mJumpingPulsers.at( partialKey + 20 ) * eTofConst::coarseClockCycle;
1862  average += mJumpingPulsers.at( partialKey + 20 ) * eTofConst::coarseClockCycle;
1863  }
1864  else {
1865  //LOG_INFO << gbtxId << " ### case 1 (3) ### " << endm;
1866  if( times.at( 0 ) < times.at( 1 ) ) {
1867  times.at( 0 ) += eTofConst::coarseClockCycle;
1868  average += eTofConst::coarseClockCycle;
1869  }
1870  else {
1871  times.at( 1 ) += eTofConst::coarseClockCycle;
1872  average += eTofConst::coarseClockCycle;
1873  }
1874  }
1875  }
1876  else if( times.at( 0 ) && times.at( 2 ) > 0 && fabs( fabs( times.at( 0 ) - times.at( 2 ) ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1877  if( mJumpingPulsers.count( partialKey + 10 ) ) {
1878  //LOG_INFO << gbtxId << " ### case 2 (1) ### " << endm;
1879  times.at( 0 ) += mJumpingPulsers.at( partialKey + 10 ) * eTofConst::coarseClockCycle;
1880  average += mJumpingPulsers.at( partialKey + 10 ) * eTofConst::coarseClockCycle;
1881  }
1882  else if( mJumpingPulsers.count( partialKey + 30 ) ) {
1883  //LOG_INFO << gbtxId << " ### case 2 (2) ### " << endm;
1884  times.at( 2 ) += mJumpingPulsers.at( partialKey + 30 ) * eTofConst::coarseClockCycle;
1885  average += mJumpingPulsers.at( partialKey + 30 ) * eTofConst::coarseClockCycle;
1886  }
1887  else {
1888  //LOG_INFO << gbtxId << " ### case 2 (3) ### " << endm;
1889  if( times.at( 0 ) < times.at( 2 ) ) {
1890  times.at( 0 ) += eTofConst::coarseClockCycle;
1891  average += eTofConst::coarseClockCycle;
1892  }
1893  else {
1894  times.at( 2 ) += eTofConst::coarseClockCycle;
1895  average += eTofConst::coarseClockCycle;
1896  }
1897  }
1898  }
1899  else if( times.at( 1 ) > 0 && times.at( 2 ) > 0 && fabs( fabs( times.at( 1 ) - times.at( 2 ) ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1900  if( mJumpingPulsers.count( partialKey + 20 ) ) {
1901  //LOG_INFO << gbtxId << " ### case 3 (1) ### " << endm;
1902  times.at( 1 ) += mJumpingPulsers.at( partialKey + 20 ) * eTofConst::coarseClockCycle;
1903  average += mJumpingPulsers.at( partialKey + 20 ) * eTofConst::coarseClockCycle;
1904  }
1905  else if( mJumpingPulsers.count( partialKey + 30 ) ) {
1906  //LOG_INFO << gbtxId << " ### case 3 (2) ### " << endm;
1907  times.at( 2 ) += mJumpingPulsers.at( partialKey + 30 ) * eTofConst::coarseClockCycle;
1908  average += mJumpingPulsers.at( partialKey + 30 ) * eTofConst::coarseClockCycle;
1909  }
1910  else {
1911  //LOG_INFO << gbtxId << " ### case 3 (3) ### " << endm;
1912  if( times.at( 1 ) < times.at( 2 ) ) {
1913  times.at( 1 ) += eTofConst::coarseClockCycle;
1914  average += eTofConst::coarseClockCycle;
1915  }
1916  else {
1917  times.at( 2 ) += eTofConst::coarseClockCycle;
1918  average += eTofConst::coarseClockCycle;
1919  }
1920  }
1921  }
1922  }
1923 
1924 
1925  if( nAv >= 2 ) {
1926  average /= nAv;
1927  }
1928 
1929  if( mDoQA && referenceTime != 0 ) {
1930  mHistograms.at( "pulserDigiTimeDiff_perGbtx" )->Fill( gbtxId * eTofConst::nCounters + 1.5, average - referenceTime );
1931  }
1932 
1933  for( int counter = eTofConst::counterStart; counter <= eTofConst::counterStop; counter++ ) {
1934  double diffToAv = 0.;
1935 
1936  if( times.at( counter - eTofConst::counterStart ) != 0. ) {
1937  diffToAv = times.at( counter - eTofConst::counterStart ) - average;
1938 
1939  if( fabs( diffToAv ) > 0.2 ) diffToAv = 0.; //removing didn't work
1940 
1941  if( mDoQA ) {
1942  mHistograms.at( "pulserDigiTimeDiff_toAverage" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, diffToAv );
1943  }
1944  }
1945 
1946  if( average != 0. ) {//only allow counter pulsers that are now close to average
1947  //times.at( counter - 1 ) = pulserTimes.at( key ) + mPulserTimeDiffGbtx.at( key )
1948  pulserTimes[ partialKey + 10 * counter ] = average + diffToAv - mPulserTimeDiffGbtx.at( partialKey + 10 * counter ); //restores original pulser times INCLUDING GBTX offset ?!?!
1949  //pulserTimes[ partialKey + 10 * counter ] = average;
1950  }
1951  else {
1952  if( pulserTimes.count( partialKey + 10 * counter ) ) {
1953  pulserTimes.erase( partialKey + 10 * counter );
1954  }
1955  }
1956  }
1957  }
1958  }
1959 
1960 
1961  // calculate difference to the reference
1962  referenceTime = 0.;
1963  if( pulserTimes.count( mReferencePulserIndex ) ) {
1964  if( mPulserTimeDiff.count( mReferencePulserIndex ) ){
1965  referenceTime = pulserTimes.at( mReferencePulserIndex ) - mPulserTimeDiff[ mReferencePulserIndex ];
1966  //LOG_INFO << "time of reference pulser updated: " << referenceTime << " with reference correction "<< mPulserTimeDiff[ mReferencePulserIndex ] << endm;
1967  }else{
1968  referenceTime = pulserTimes.at( mReferencePulserIndex );
1969  //LOG_INFO << "time of reference pulser updated: " << referenceTime << endm;
1970  }
1971  }
1972 
1973 
1974 
1975  if( referenceTime != 0 ) {
1976  int iLockThreshold = 0;
1977  mHistograms[ "pulserDigiTimeDiff_RefCorr" ]->Reset("ICESM");
1978 
1979  for( auto& kv : pulserTimes ) {
1980  // check if new pulser time difference seems reasonable ( only allow jumps in multiple of the coarse clock tick ) to avoid smeared pulsers
1981  if( mPulserTimeDiff.count( kv.first ) ) {//pulser time difference default available previous events
1982  //double modDist = fmod( mPulserTimeDiff.at( kv.first ) - ( kv.second - referenceTime ), eTofConst::coarseClockCycle );
1983  //modDist = modDist - eTofConst::coarseClockCycle * round( modDist / eTofConst::coarseClockCycle );
1984 
1985  double modDist = mPulserTimeDiff.at( kv.first ) - ( kv.second - referenceTime ); //test PW
1986  //modDist = modDist - eTofConst::coarseClockCycle * round( modDist / eTofConst::coarseClockCycle );
1987 
1988 
1989  if( fabs( modDist ) > 0.2 ) {
1990  mUnlockPulserState[ kv.first ]++;
1991 
1992  // LOG_INFO << " pulser time " << kv.first << " seems unreasonable (" << kv.second - referenceTime << ")";
1993  // LOG_INFO << " compared to previous stored value (" << mPulserTimeDiff.at( kv.first ) << ")" << endm;
1994 
1995  // only unlock pulser state if 10 consecutive events have a modDist larger then the threshold
1996  if( mUnlockPulserState.at( kv.first ) < iLockThreshold ) {
1997  // LOG_INFO << " --> ignore for now and move on" << endm;
1998  continue; //move on, don't update pulser times!
1999  }
2000  else{
2001  // LOG_INFO << " --> pulser state has been unlocked" << endm;
2002  }
2003 
2004  //fill 2d Hist here with GBTX and counter number
2005  mHistograms[ "pulserDigiTimeDiff_RefCorr" ]->Fill( kv.second - referenceTime );
2006 
2007  }
2008  else{
2009  if( mUnlockPulserState.count( kv.first ) ) {
2010  LOG_INFO << " --> new event doesn't have offset for pulser " << kv.first << " --> remove the entry" << endm;
2011  mUnlockPulserState.erase( kv.first );
2012  }
2013  }
2014  }
2015  //pulser time differece was set here!
2016 
2017 
2018  if( mDoQA ) {
2019  int sector = kv.first / 1000;
2020  int zPlane = ( kv.first % 1000 ) / 100;
2021  int counter = ( kv.first % 100 ) / 10;
2022  int side = kv.first % 10;
2023 
2024  int gbtxId = ( sector - eTofConst::sectorStart ) * ( eTofConst::nPlanes * eTofConst::nSides )
2025  + ( zPlane - eTofConst::zPlaneStart ) * eTofConst::nSides
2026  + ( side - eTofConst::sideStart );
2027 
2028  if( mPulserTimeDiff.count( kv.first ) ) {
2029  mHistograms.at( "pulserDigiTimeDiff_fullCorr" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, mPulserTimeDiff.at( kv.first ) );
2030  }
2031 
2032  mHistograms[ "pulserDigiPresence" ]->Fill(gbtxId * eTofConst::nCounters + counter - 0.5);
2033  }
2034  }
2035 
2036  //LOG_INFO << "Check " << referenceTime << endm;
2037  if( mDoQA ) {
2038  mHistograms[ "pulserDigiPresence" ]->Fill( -1 ); //use as event counter
2039  mHistograms[ "pulserDigiTimeDiff_CorrAgreement" ]->Fill( mHistograms[ "pulserDigiTimeDiff_RefCorr" ]->GetMaximum() );
2040  mHistograms[ "pulserDigiTimeDiff_CorrCommonJump" ]->Fill( mHistograms[ "pulserDigiTimeDiff_RefCorr" ]->GetEntries() );
2041  }
2042  if( ( mHistograms[ "pulserDigiTimeDiff_RefCorr" ]->GetEntries() > 150 ) && ( mHistograms[ "pulserDigiTimeDiff_RefCorr" ]->GetMaximum() > 15 ) ){
2043 
2044  //LOG_INFO << "Check: found " << mHistograms[ "pulserDigiTimeDiff_RefCorr" ]->GetMaximum() <<" agreeing time shifts. Shifting reference times " << endm;
2045  int iMaximumBin = mHistograms[ "pulserDigiTimeDiff_RefCorr" ]->GetMaximumBin();
2046  double dRefCorrAverage = 0;
2047  int iNRefCorrAgree = 0;
2048 
2049  for( auto& kv : pulserTimes ) { //build average of all agreeing pulsers
2050  if ( mHistograms[ "pulserDigiTimeDiff_RefCorr" ]->FindBin(kv.second - referenceTime) == iMaximumBin ){
2051  dRefCorrAverage += kv.second - referenceTime;
2052  iNRefCorrAgree++;
2053  }
2054  }
2055  dRefCorrAverage /= iNRefCorrAgree;
2056  referenceTime -= dRefCorrAverage;
2057  //pulserTimes.at( mReferencePulserIndex ) -= dRefCorrAverage;
2058  mPulserTimeDiff[ mReferencePulserIndex ] -= dRefCorrAverage;
2059 
2060  }else{
2061  //LOG_INFO << "Check: found only " << mHistograms[ "pulserDigiTimeDiff_RefCorr" ]->GetMaximum() <<" agreeing time shifts. Shifting pulser times " << endm;
2062  for( auto& kv : pulserTimes ) {
2063  if ( ! mPulserTimeDiff[ kv.first ] ) {
2064  mPulserTimeDiff[ kv.first ] = kv.second - referenceTime;
2065  continue;
2066  }
2067  if( mUnlockPulserState.count( kv.first ) ){
2068  if( mUnlockPulserState.at( kv.first ) > iLockThreshold ){// check if pulser is locked
2069  mPulserTimeDiff[ kv.first ] = kv.second - referenceTime;
2070  }
2071  }
2072  }
2073  }
2074  }
2075 }
2076 
2077 
2078 //_____________________________________________________________
2082 void
2083 StETofCalibMaker::applyCalibration( StETofDigi* aDigi, StETofHeader* etofHeader )
2084 {
2085  int key = aDigi->sector() * 100000 + aDigi->zPlane() * 10000 + aDigi->counter() * 1000 + aDigi->strip() * 10 + aDigi->side();
2086  if( !mStatus.count( key) || mStatus.at( key ) != 1 ) {
2087  if( mDebug ) {
2088  LOG_DEBUG << "status of channel with key " << key << " was not ok ---> skip calibrating this digi" << endm;
2089  }
2090  return;
2091  }
2092 
2093  // ignore digis flaged as pulsers ( calibTot = -999. )
2094  if( fabs( aDigi->calibTot() + 999. ) < 1.e-5 ) {
2095  if( mDebug ) {
2096  LOG_DEBUG << "digi flaged as pulser --> skip" << endm;
2097  }
2098  return;
2099  }
2100 
2101  float timeToTrigger = aDigi->rawTime() - mTriggerTime;
2102 
2103  // check if digi is inside the timing window and only calibrate those, do nothing digis outside the window ( calibTime = 0, calibTot = -1 )
2104  if( timeToTrigger > mTimingWindow.at( aDigi->rocId() ).first &&
2105  timeToTrigger < mTimingWindow.at( aDigi->rocId() ).second )
2106  {
2107 
2108  if( mStrictPulserHandling ){
2109  int PulserKey = aDigi->sector() * 1000 + aDigi->zPlane() * 100 + aDigi->side() + 10 * aDigi->counter();
2110  if( !mPulserPresent.count( PulserKey ) ) {
2111  if( mDebug ) {
2112  LOG_DEBUG << "no pulser in the same event for this counter --> digi skipped due to strict pulser handling" << endm;
2113  }
2114  return;
2115  }
2116  }
2117 
2118  double calibTot = aDigi->rawTot() * mGet4TotBinWidthNs * calibTotFactor( aDigi );
2119 
2120  aDigi->setCalibTot( calibTot );
2121 
2122  double calibTime = aDigi->rawTime() - mResetTime
2123  - resetTimeCorr()
2124  - calibTimeOffset( aDigi )
2125  - slewingTimeOffset( aDigi )
2126  - applyPulserOffset( aDigi );
2127 
2128  aDigi->setCalibTime( calibTime );
2129 
2130  if( mDebug ) {
2131  // print out the new information
2132  LOG_DEBUG << "raw Time, ToT: " << aDigi->rawTime() << ", " << aDigi->rawTot() << endm;
2133  LOG_DEBUG << "calibrated Time, ToT: " << aDigi->calibTime() << ", " << aDigi->calibTot() << endm;
2134  }
2135 
2136  }
2137  else{
2138  if( mDebug ) {
2139  LOG_DEBUG << "digi is outside the timing window (time to trigger = " << timeToTrigger << ") --> skip" << endm;
2140  }
2141  }
2142 }
2143 //_____________________________________________________________
2144 
2145 
2146 //_____________________________________________________________
2147 void
2148 StETofCalibMaker::resetToRaw( StMuETofDigi* aDigi )
2149 {
2150  if( !aDigi ) return;
2151 
2152  aDigi->setGeoAddress( 0, 0, 0, 0, 0 );
2153  aDigi->setCalibTime( 0. );
2154  aDigi->setCalibTot( -1. );
2155 
2156  aDigi->setAssociatedHitId( -1 );
2157 }
2158 
2159 
2160 //_____________________________________________________________
2161 void
2162 StETofCalibMaker::applyMapping( StMuETofDigi* aDigi )
2163 {
2164  applyMapping( ( StETofDigi* ) aDigi );
2165 }
2166 
2167 
2168 //_____________________________________________________________
2169 void
2170 StETofCalibMaker::flagPulserDigis( StMuETofDigi* aDigi, unsigned int index, std::map< unsigned int, std::vector< unsigned int > >& pulserDigiMap )
2171 {
2172  flagPulserDigis( ( StETofDigi* ) aDigi, index, pulserDigiMap );
2173 }
2174 
2175 
2176 //_____________________________________________________________
2177 void
2178 StETofCalibMaker::applyCalibration( StMuETofDigi* aDigi, StMuETofHeader* etofHeader )
2179 {
2180  applyCalibration( ( StETofDigi* ) aDigi, ( StETofHeader* ) etofHeader );
2181 }
2182 //_____________________________________________________________
2183 
2184 
2185 //_____________________________________________________________
2189 double
2190 StETofCalibMaker::calibTotFactor( StETofDigi* aDigi )
2191 {
2192  unsigned int key = aDigi->sector() * 100 + aDigi->zPlane() * 10 + aDigi->counter();
2193  unsigned int bin = aDigi->strip() + 32 * ( aDigi->side() - 1 );
2194 
2195  if( mDigiTotCorr.count( key ) ) {
2196  float binContent = mDigiTotCorr.at( key )->GetBinContent( bin );
2197 
2198  if( fabs( binContent ) > 1e-5 ) {
2199  if( mDebug ) {
2200  LOG_DEBUG << "calibTotFactor: histogram with key " << key << " at bin " << bin << " -> return bin content: " << binContent << endm;
2201  }
2202  return (1.0/binContent); //invert here to get to fixed mean value!
2203  }
2204  else {
2205  if( mDebug ) {
2206  LOG_WARN << "calibTotFactor: histogram with key " << key << " at bin " << bin << " has content of 0 -> return 1" << endm;
2207  }
2208  return 1.;
2209  }
2210  }
2211  else {
2212  if( mDebug ) {
2213  LOG_WARN << "calibTotFactor: required histogram with key " << key << " doesn't exist -> return 1" << endm;
2214  }
2215  return 1.;
2216  }
2217 }
2218 
2219 
2220 //_____________________________________________________________
2224 double
2225 StETofCalibMaker::calibTimeOffset( StETofDigi* aDigi )
2226 {
2227  unsigned int key = aDigi->sector() * 100 + aDigi->zPlane() * 10 + aDigi->counter();
2228  unsigned int bin = aDigi->strip() + 32 * ( aDigi->side() - 1 );
2229 
2230  if( mDigiTimeCorr.count( key ) ) {
2231  float binContent = mDigiTimeCorr.at( key )->GetBinContent( bin );
2232  if( mDebug ) {
2233  LOG_DEBUG << "calibTimeOffset: histogram with key " << key << " at bin " << bin << " -> return bin content: " << binContent << endm;
2234  }
2235  return binContent;
2236  }
2237  else {
2238  if( mDebug ) {
2239  LOG_WARN << "calibTimeOffset: required histogram with key " << key << " doesn't exist -> return 0" << endm;
2240  }
2241  return 0.;
2242  }
2243 }
2244 
2245 
2246 //_____________________________________________________________
2250 double
2251 StETofCalibMaker::slewingTimeOffset( StETofDigi* aDigi )
2252 {
2253  unsigned int key = aDigi->sector() * 100000 + aDigi->zPlane() * 10000 + aDigi->counter() * 1000 + aDigi->strip() * 10 + aDigi->side();
2254 
2255  if( mDigiSlewCorr.count( key ) ) {
2256 
2257  unsigned int totBin = mDigiSlewCorr.at( key )->FindBin( aDigi->rawTot() ); //adjusted. PW
2258 mDebug = true;
2259  if( mDigiSlewCorr.at( key )->GetBinEntries( totBin ) <= mMinDigisPerSlewBin && totBin < etofSlewing::nTotBins ) {
2260  if( mDebug ) {
2261  LOG_DEBUG << "slewingTimeOffset: insufficient statistics for slewing calibration in channel " << key << " at tot bin " << totBin << " --> return 0" << endm;
2262  }
2263  return 0.;
2264  }
2265 
2266  float val = mDigiSlewCorr.at( key )->Interpolate( aDigi->rawTot() ); //adjusted. PW
2267  if( mDebug ) {
2268  LOG_DEBUG << "slewingTimeOffset: histogram with key " << key << " with calib TOT of " << aDigi->calibTot() << " --> interpolated correction: " << val << endm;
2269  }
2270  return val;
2271  }
2272  else {
2273  if( mDebug ) {
2274  LOG_DEBUG << "slewingTimeOffset: required histogram with key " << key << " doesn't exist -> return 0" << endm;
2275  }
2276  return 0.;
2277  }
2278 }
2279 
2280 
2281 //_____________________________________________________________
2285 double
2286 StETofCalibMaker::applyPulserOffset( StETofDigi* aDigi )
2287 {
2288  int key = aDigi->sector() * 1000 + aDigi->zPlane() * 100 + aDigi->counter() * 10 + aDigi->side();
2289 
2290  if( !mPulserTimeDiff.count( key )) {
2291  return 0.;
2292  }
2293 
2294  return mPulserTimeDiff.at( key );
2295 }
2296 
2297 
2298 //_____________________________________________________________
2302 double
2303 StETofCalibMaker::triggerTime( StETofHeader* header )
2304 {
2305  // initialize trigger time with the one from the header in case the map of trigger time stamps per AFCK is empty
2306  double triggerTime = header->trgGdpbFullTime();
2307 
2308  // count the occurance of a given trigger time stamp in the GdbpTs map of the eTOF header
2309  std::map< uint64_t, short > countsGdpbTs;
2310  for( const auto& kv : header->rocGdpbTs() ) {
2311  if( mDebug ) {
2312  LOG_DEBUG << "triggerTime (" << std::hex << "Ox" << kv.first << std::dec << ") " << kv.second * eTofConst::coarseClockCycle * 1.e-9 << endm;
2313  }
2314  ++countsGdpbTs[ kv.second ];
2315  }
2316 
2317  // combine adjacent trigger times to get the number of right trigger time stamps without outliers
2318  // take the trigger Ts that occured most often in the combined counting map
2319  short maxCount = 0;
2320  short accCount = 0;
2321  uint64_t mostProbableTriggerTs = 0;
2322 
2323  for( auto it = countsGdpbTs.begin(); it != countsGdpbTs.end(); it++ ) {
2324  auto next = std::next( it, 1 );
2325  auto prev = std::prev( it, 1 );
2326 
2327  short countTs = it->second;
2328 
2329  if( next != countsGdpbTs.end() && ( next->first - it->first ) == 1 ) {
2330  countTs += next->second;
2331  }
2332  if( accCount > 0 && ( it->first - prev->first ) == 1 ) {
2333  countTs += prev->second;
2334  }
2335 
2336  if( countTs >= accCount ) {
2337  accCount = countTs;
2338 
2339  if( it->second > maxCount ) {
2340  maxCount = it->second;
2341  mostProbableTriggerTs = it->first;
2342  }
2343  }
2344  }
2345 
2346  if( mostProbableTriggerTs > 0 ) {
2347  triggerTime = mostProbableTriggerTs * eTofConst::coarseClockCycle;
2348  }
2349 
2350  if( mDebug ) {
2351  LOG_DEBUG << "trigger TS: " << mostProbableTriggerTs << " --> trigger time (ns): " << triggerTime << endm;
2352  }
2353 
2354  return triggerTime;
2355 }
2356 
2357 
2358 //_____________________________________________________________
2362 double
2363 StETofCalibMaker::resetTime( StETofHeader* header )
2364 {
2365  // count the occurance of a given reset time stamp in the StarTs map of the eTOF header
2366  std::map< uint64_t, short > countsStarTsRaw;
2367  for( const auto& kv : header->rocStarTs() ) {
2368  if( mDebug ) {
2369  LOG_DEBUG << "resetTime (" << std::hex << "Ox" << kv.first << std::dec << ") " << kv.second * eTofConst::coarseClockCycle * 1.e-9 << endm;
2370  }
2371 
2372  // in Run18 only one of the AFCKs was giving the correct reset time: 0x18e6
2373  if( mRunYear == 2018 && kv.first != 0x18e6 ) continue;
2374 
2375  if( kv.second != 0 ) {
2376  ++countsStarTsRaw[ kv.second ];
2377  }
2378  }
2379 
2380 
2381  // combine adjacent reset time values with the earlier one
2382  std::map< uint64_t, short > countsStarTs;
2383  for( auto it = countsStarTsRaw.begin(); it != countsStarTsRaw.end(); it++ ) {
2384  auto next = std::next( it, 1 );
2385 
2386  short countTs = it->second;
2387 
2388  if( next != countsStarTsRaw.end() && ( next->first - it->first ) == 1 ) {
2389  countTs += next->second;
2390  }
2391 
2392  countsStarTs[ it->first ] = countTs;
2393  }
2394 
2395 
2396 
2397 
2398  if( mDoQA ) {
2399  if( countsStarTs.size() == 0 ) {
2400  LOG_INFO << "all AFCKs report a reset time value 0" << endm;
2401  }
2402 
2403  for( const auto& kv : countsStarTs ) {
2404  LOG_DEBUG << "resetTime cand:" << kv.first << " (" << kv.second << " times)" << endm;
2405  mHistograms.at( "resetTimeCand_times" )->Fill( kv.second );
2406  }
2407 
2408  for( const auto& kv : header->rocStarTs() ) {
2409  unsigned int sector;
2410  mHwMap->mapToSector( kv.first, sector );
2411 
2412  LOG_DEBUG << "resetTime(" << sector << "): " << kv.second << endm;
2413 
2414  std::string histName = "resetTimeDifferenceToSector" + to_string( sector );
2415  for( const auto& jv : header->rocStarTs() ) {
2416  unsigned int sec;
2417  mHwMap->mapToSector( jv.first, sec );
2418 
2419  mHistograms.at( histName )->Fill( sec, ( int64_t ) jv.second - ( int64_t ) kv.second );
2420  }
2421  }
2422  }
2423 
2424 
2425  while( countsStarTs.size() > 0 ) {
2426  auto it = std::max_element( countsStarTs.begin(), countsStarTs.end(),
2427  []( const pair< uint64_t, short >& p1, const pair< uint64_t, short >& p2 ) {
2428  return p1.second < p2.second; } );
2429 
2430  double resetTime = it->first * eTofConst::coarseClockCycle;
2431 
2432 
2433  // only update reset time if it is at least two clock ticks away from the old reset time to avoid jitter
2434  if( abs( mResetTs - ( int64_t ) it->first ) < 2 ) {
2435  resetTime = mResetTs * eTofConst::coarseClockCycle;
2436  }
2437  else {
2438  LOG_INFO << "change in reset TS: " << mResetTs << " --> " << it->first << endm;
2439  mResetTs = it->first;
2440  }
2441 
2442 
2443  // Run19: trigger - reset time should be on the order of a few second up to 120 minutes (7.2*10^12 ns), i.e. max. run length
2444  // Run20: difference can be negative due to eTOF DAQ restarts at the beginning of runs while eTOF is put to "BUSY" in run control
2445  if( mTriggerTime - resetTime < 7.2e12 ) {
2446  if( mDebug ) {
2447  LOG_DEBUG << "reset time (ns): " << resetTime << " --> difference to trigger time in seconds: " << ( mTriggerTime - resetTime ) * 1.e-9 << endm;
2448  }
2449  LOG_DEBUG << "--> picked reset TS:" << mResetTs << endm;
2450 
2451  if( mDoQA ) {
2452  mHistograms.at( "resetTimeCand_picked" )->Fill( it->second );
2453 
2454  auto rawIt = std::max_element( countsStarTsRaw.begin(), countsStarTsRaw.end(),
2455  []( const pair< uint64_t, short >& p1, const pair< uint64_t, short >& p2 ) {
2456  return p1.second < p2.second; } );
2457 
2458  mHistograms.at( "resetTimeCand_compare" )->Fill( ( int64_t ) mResetTs - ( int64_t ) rawIt->first );
2459  mHistograms.at( "resetTimeCand_value" )->Fill( mResetTs % ( int ) eTofConst::bTofClockCycle );
2460  }
2461 
2462  return resetTime;
2463  }
2464  else {
2465  countsStarTs.erase( it );
2466  }
2467  }
2468 
2469  return 0.;
2470 }
2471 
2472 
2473 //_____________________________________________________________
2474 unsigned int
2475 StETofCalibMaker::channelToKey( const unsigned int channelId ) {
2476  unsigned int sector = ( channelId / eTofConst::nChannelsPerSector ) + eTofConst::sectorStart;
2477  unsigned int zPlane = ( ( channelId % eTofConst::nChannelsPerSector ) / eTofConst::nChannelsPerModule ) + eTofConst::zPlaneStart;
2478  unsigned int counter = ( ( channelId % eTofConst::nChannelsPerModule ) / eTofConst::nChannelsPerCounter ) + eTofConst::counterStart;
2479  unsigned int strip = ( ( channelId % eTofConst::nChannelsPerCounter ) / eTofConst::nSides ) + eTofConst::stripStart;
2480  unsigned int side = ( channelId % eTofConst::nSides ) + eTofConst::sideStart;
2481 
2482  return sector * 100000 + zPlane * 10000 + counter * 1000 + strip * 10 + side;
2483 }
2484 
2485 
2486 //_____________________________________________________________
2487 unsigned int
2488 StETofCalibMaker::detectorToKey( const unsigned int detectorId ) {
2489  unsigned int sector = ( detectorId / eTofConst::nCountersPerSector ) + eTofConst::sectorStart;
2490  unsigned int zPlane = ( ( detectorId % eTofConst::nCountersPerSector ) / eTofConst::nCounters ) + eTofConst::zPlaneStart;
2491  unsigned int counter = ( detectorId % eTofConst::nCounters ) + eTofConst::counterStart;
2492 
2493  return sector * 100 + zPlane * 10 + counter;
2494 }
2495 
2496 
2497 //_____________________________________________________________
2498 unsigned int
2499 StETofCalibMaker::sideToKey( const unsigned int sideId ) {
2500  unsigned int sector = ( sideId / ( eTofConst::nCountersPerSector * eTofConst::nSides ) ) + eTofConst::sectorStart;
2501  unsigned int zPlane = ( ( sideId % ( eTofConst::nCountersPerSector * eTofConst::nSides ) ) / ( eTofConst::nCounters * eTofConst::nSides ) ) + eTofConst::zPlaneStart;
2502  unsigned int counter = ( ( sideId % ( eTofConst::nCounters * eTofConst::nSides ) ) / eTofConst::nSides ) + eTofConst::counterStart;
2503  unsigned int side = ( sideId % eTofConst::nSides ) + eTofConst::sideStart;
2504 
2505  return sector * 1000 + zPlane * 100 + counter * 10 + side;
2506 }
2507 
2508 
2509 
2510 //_____________________________________________________________
2511 // setHistFileName uses the string argument from the chain being run to set
2512 // the name of the output histogram file.
2513 void
2514 StETofCalibMaker::setHistFileName()
2515 {
2516  string extension = ".etofCalib.root";
2517 
2518  if( GetChainOpt()->GetFileOut() != nullptr ) {
2519  TString outFile = GetChainOpt()->GetFileOut();
2520 
2521  mHistFileName = ( string ) outFile;
2522 
2523  // get rid of .root
2524  size_t lastindex = mHistFileName.find_last_of( "." );
2525  mHistFileName = mHistFileName.substr( 0, lastindex );
2526 
2527  // get rid of .MuDst or .event if necessary
2528  lastindex = mHistFileName.find_last_of( "." );
2529  mHistFileName = mHistFileName.substr( 0, lastindex );
2530 
2531  // get rid of directories
2532  lastindex = mHistFileName.find_last_of( "/" );
2533  mHistFileName = mHistFileName.substr( lastindex + 1, mHistFileName.length() );
2534 
2535  mHistFileName = mHistFileName + extension;
2536  } else {
2537  LOG_ERROR << "Cannot set the output filename for histograms" << endm;
2538  }
2539 }
2540 
2541 //_____________________________________________________________
2542 void
2543 StETofCalibMaker::bookHistograms()
2544 {
2545  LOG_INFO << "bookHistograms() ... " << endm;
2546 
2547  mHistograms[ "pulserDigiTimeDiff_GbtxCorrProf" ] = new TProfile( "pulserDigiTimeDiff_GbtxCorrProf", "time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216, "s" );
2548  mHistograms[ "pulserDigiTimeDiff_GbtxCorrProfMod" ] = new TProfile( "pulserDigiTimeDiff_GbtxCorrProfMod", "time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216, "s" );
2549  mHistograms[ "pulserDigiTimeDiff_fullCorr" ] = new TH2F( "pulserDigiTimeDiff_fullCorr", "time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2550  mHistograms[ "pulserDigiTimeDiff_RefCorr" ] = new TH1F("pulserDigiTimeDiff_RefCorr", "time difference of pulsers to reference; #Delta T (ns)", 45, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ));
2551 
2552  if( mDoQA ) {
2553  mHistograms[ "pulserDigiTimeDiff" ] = new TH2F( "pulserDigiTimeDiff", "time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2554  mHistograms[ "pulserDigiTimeToCounter1" ] = new TH2F( "pulserDigiTimeToCounter1", "time difference of pulsers to counter 1;pulser channel;#Delta T (ns)", 216, 0, 216, 2*360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2555  mHistograms[ "pulserDigiTimeToCounter2" ] = new TH2F( "pulserDigiTimeToCounter2", "time difference of pulsers to counter 2;pulser channel;#Delta T (ns)", 216, 0, 216, 2*360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2556  mHistograms[ "pulserDigiTimeToCounter3" ] = new TH2F( "pulserDigiTimeToCounter3", "time difference of pulsers to counter 3;pulser channel;#Delta T (ns)", 216, 0, 216, 2*360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2557 
2558  mHistograms[ "pulserDigiTimeDiff_GbtxCorr" ] = new TH2F( "pulserDigiTimeDiff_GbtxCorr", "time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2559  mHistograms[ "pulserDigiTimeDiff_perGbtx" ] = new TH2F( "pulserDigiTimeDiff_perGbtx", "average time difference of pulsers in one Gbtx to reference;pulser channel;#Delta T (ns)", 216/3, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2560  mHistograms[ "pulserDigiTimeDiff_toAverage" ] = new TH2F( "pulserDigiTimeDiff_toAverage", "time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216, 4*360, -359.5 * ( 6.25 / 112 ), 360.5 * ( 6.25 / 112 ) );
2561 
2562  mHistograms[ "1_off" ] = new TH2F( "1_off", "average time difference of pulsers in one Gbtx to reference;pulser channel;#Delta T (ns)", 216/3, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2563  mHistograms[ "2_off" ] = new TH2F( "2_off", "average time difference of pulsers in one Gbtx to reference;pulser channel;#Delta T (ns)", 216/3, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2564  mHistograms[ "3_off" ] = new TH2F( "3_off", "average time difference of pulsers in one Gbtx to reference;pulser channel;#Delta T (ns)", 216/3, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2565 
2566  mHistograms[ "pulserDigiTimeDiff_CorrAgreement" ] = new TH1F("pulserDigiTimeDiff_CorrAgreement", "Number of pulsers agreeing on a common shift between events; #pulsers", 218, -0.5, 217.5);
2567  mHistograms[ "pulserDigiTimeDiff_CorrCommonJump" ] = new TH1F("pulserDigiTimeDiff_CorrCommonJump", "Number of pulsers jumping at the same time between events; #pulsers", 218, -0.5, 217.5);
2568  mHistograms[ "pulserDigiPresence" ] = new TH1F( "pulserDigiPresence", "pulser presence (number of events at ( -1 );pulser channel", 218, -1.5, 216.5);
2569 
2570  for( int i=0; i<12; i++ ) {
2571  std::string histName = "resetTimeDifferenceToSector" + to_string( i + 13 );
2572  mHistograms[ histName ] = new TH2F( Form( "resetTimeDifferenceToSector%d", i + 13 ), Form("reset time difference to sector %d;sector;#DeltaT (clock ticks)", i + 13 ), 12, 12.5, 24.4, 5, -2.5, 2.5 );
2573  }
2574  mHistograms[ "ETOF_QA_daqMissmatches_get4" ] = new TProfile( "ETOF_QA_daqMissmatches_get4", "missmatch percentage for each get4; get4 []; missmatch percentage", 1728, 0.5, 1728.5 );
2575  mHistograms[ "ETOF_QA_daqMissmatches_counter" ] = new TProfile( "ETOF_QA_daqMissmatches_counter", "missmatch percentage for each counter; counter []; missmatch percentage", 108, 0.5, 108.5 );
2576  mHistograms[ "resetTimeCand_times" ] = new TH1F( "resetTimeCand_times", "sectors agreeing on reset time candidates;# sectors with common candidate;# events", 12, 0.5, 12.5 );
2577  mHistograms[ "resetTimeCand_picked" ] = new TH1F( "resetTimeCand_picked", "sectors agreeing on picked reset time;# sectors with picked reset time;# events", 12, 0.5, 12.5 );
2578  mHistograms[ "resetTimeCand_compare" ] = new TH1F( "resetTimeCand_compare", "difference between old and new way;#DeltaT (clock ticks);# events", 5, -2.5, 2.5 );
2579  mHistograms[ "resetTimeCand_value" ] = new TH1F( "resetTimeCand_value", "picked reset time value;clock ticks;# events", ( int ) eTofConst::bTofClockCycle, 0.5, 0.5 + eTofConst::bTofClockCycle );
2580  }
2581 
2582  for ( auto& kv : mHistograms ) {
2583  kv.second->SetDirectory( 0 );
2584  }
2585 }
2586 
2587 //_____________________________________________________________
2588 void
2589 StETofCalibMaker::writeHistograms()
2590 {
2591  if( mHistFileName != "" ) {
2592  LOG_INFO << "writing histograms to: " << mHistFileName.c_str() << endm;
2593 
2594  TFile histFile( mHistFileName.c_str(), "RECREATE", "etofCalib" );
2595  histFile.cd();
2596 
2597  for( int i=0; i<12; i++ ) {
2598  std::string histName = "resetTimeDifferenceToSector" + to_string( i + 13 );
2599  mHistograms.at( histName )->Scale( 12. / mHistograms.at( histName )->GetEntries() );
2600  }
2601 
2602  for ( const auto& kv : mHistograms ) {
2603  if( kv.second->GetEntries() > 0 ) kv.second->Write();
2604  }
2605 
2606  histFile.Close();
2607  }
2608  else {
2609  LOG_INFO << "histogram file name is empty string --> cannot write histograms" << endm;
2610  }
2611 }
unsigned int sector() const
Sector.
Definition: StMuETofDigi.h:212
double calibTime() const
calibrated time
Definition: StETofDigi.h:206
unsigned int zPlane() const
ZPlane.
Definition: StETofDigi.h:213
static StMuETofHeader * etofHeader()
returns pointer to the StMuETofHeader
Definition: StMuDst.h:431
unsigned int side() const
Side.
Definition: StMuETofDigi.h:217
unsigned int strip() const
Strip.
Definition: StMuETofDigi.h:215
double rawTot() const
Getter for uncalibrated Tot.
Definition: StMuETofDigi.h:208
unsigned int sector() const
Sector.
Definition: StETofDigi.h:212
unsigned int zPlane() const
ZPlane.
Definition: StMuETofDigi.h:213
unsigned int strip() const
Strip.
Definition: StETofDigi.h:215
double rawTime() const
Raw Time.
Definition: StETofDigi.h:205
double rawTot() const
Getter for uncalibrated Tot.
Definition: StETofDigi.h:208
static TClonesArray * etofArray(int type)
returns pointer to the n-th TClonesArray from the etof arrays // FS
Definition: StMuDst.h:289
static StMuETofDigi * etofDigi(int i)
returns pointer to the i-th StMuEtofDigi
Definition: StMuDst.h:427
unsigned int elChan() const
electronic Channel.
Definition: StETofDigi.h:218
unsigned int get4Id() const
get4Id.
Definition: StETofDigi.h:219
unsigned int counter() const
Counter.
Definition: StETofDigi.h:214
void setGeoAddress(const unsigned int iSector, const unsigned int iZPlane, const unsigned int iCounter, const unsigned int iChannel, const unsigned int iSide)
double calibTot() const
Getter for calibrated Tot.
Definition: StETofDigi.h:210
unsigned int rocId() const
RocId.
Definition: StETofDigi.h:220
std::vector< bool > missMatchFlagVec() const
Flag for each Get4 TDC to mark if it is available in this event.
unsigned int side() const
Side.
Definition: StETofDigi.h:217
StETofCalibMaker(const char *name="etofCalib")
default constructor
unsigned int counter() const
Counter.
Definition: StMuETofDigi.h:214
void setGeoAddress(const unsigned int iSector, const unsigned int iZPlane, const unsigned int iCounter, const unsigned int iChannel, const unsigned int iSide)
Definition: StETofDigi.cxx:152
std::vector< bool > missMatchFlagVec() const
Flag for each Get4 TDC to mark if it is available in this event.
Definition: Stypes.h:41
virtual TDataSet * Find(const char *path) const
Definition: TDataSet.cxx:362
double rawTime() const
Raw Time.
Definition: StMuETofDigi.h:205