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