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