56 #include "StChain/StChainOpt.h"
60 #include "StETofCollection.h"
61 #include "StETofHeader.h"
62 #include "StETofDigi.h"
64 #include "StMuDSTMaker/COMMON/StMuDst.h"
65 #include "StMuDSTMaker/COMMON/StMuETofDigi.h"
66 #include "StMuDSTMaker/COMMON/StMuETofHeader.h"
68 #include "StETofCalibMaker.h"
69 #include "StETofUtil/StETofHardwareMap.h"
70 #include "StETofUtil/StETofConstants.h"
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"
84 namespace etofSlewing {
85 const unsigned int nTotBins = 30;
88 const float conversionFactor = 1. / 100.;
98 mFileNameCalibParam(
"" ),
99 mFileNameElectronicsMap(
"" ),
100 mFileNameStatusMap(
"" ),
101 mFileNameTimingWindow(
"" ),
102 mFileNameSignalVelocity(
"" ),
103 mFileNameCalibHistograms(
"" ),
104 mFileNameResetTimeCorr(
"" ),
105 mFileNamePulserTotPeak(
"" ),
106 mFileNamePulserTimeDiffGbtx(
"" ),
108 mGet4TotBinWidthNs( 1. ),
109 mMinDigisPerSlewBin( 25 ),
110 mResetTimeCorr( 0. ),
114 mPulserPeakTime( 0. ),
115 mReferencePulserIndex( 0 ),
116 mStrictPulserHandling( 0 ),
117 mUsePulserGbtxDiff( true ),
123 LOG_DEBUG <<
"StETofCalibMaker::ctor" << endm;
126 mTimingWindow.clear();
127 mPulserWindow.clear();
128 mSignalVelocity.clear();
129 mDigiTotCorr.clear();
130 mDigiSlewCorr.clear();
132 mPulserPeakTot.clear();
133 mPulserTimeDiff.clear();
134 mPulserTimeDiffGbtx.clear();
135 mNPulsersCounter.clear();
136 mNStatusBitsCounter.clear();
137 mPulserPresent.clear();
139 mJumpingPulsers.clear();
141 mUnlockPulserState.clear();
143 mStuckFwDigi.clear();
150 StETofCalibMaker::~StETofCalibMaker()
158 StETofCalibMaker::Init()
160 LOG_INFO <<
"StETofCalibMaker::Init()" << endm;
170 StETofCalibMaker::InitRun( Int_t runnumber )
172 mRunYear = ( runnumber + 727000 ) / 1000000 + 1999;
173 LOG_INFO <<
"runnumber: " << runnumber <<
" --> year: " << mRunYear << endm;
176 std::ifstream paramFile;
192 if( mFileNameElectronicsMap.empty() ) {
193 LOG_INFO <<
"etofElectronicsMap: no filename provided --> load database table" << endm;
195 dbDataSet = GetDataBase(
"Geometry/etof/etofElectronicsMap" );
197 St_etofElectronicsMap* etofElectronicsMap =
static_cast< St_etofElectronicsMap*
> ( dbDataSet->
Find(
"etofElectronicsMap" ) );
198 if( !etofElectronicsMap ) {
199 LOG_ERROR <<
"unable to get the electronics map from the database" << endm;
206 LOG_INFO <<
"etofElectronicsMap: filename provided --> use parameter file: " << mFileNameElectronicsMap.c_str() << endm;
216 if( mFileNameStatusMap.empty() ) {
217 LOG_INFO <<
"etofStatusMap: no filename provided --> load database table" << endm;
219 dbDataSet = GetDataBase(
"Calibrations/etof/etofStatusMap" );
221 St_etofStatusMap* etofStatusMap =
static_cast< St_etofStatusMap*
> ( dbDataSet->
Find(
"etofStatusMap" ) );
222 if( !etofStatusMap ) {
223 LOG_ERROR <<
"unable to get the status map from the database" << endm;
227 etofStatusMap_st* statusMapTable = etofStatusMap->GetTable();
229 for(
size_t i=0; i<eTofConst::nChannelsInSystem; i++ ) {
230 mStatus[ channelToKey( i ) ] = (int) statusMapTable->status[ i ];
234 LOG_INFO <<
"etofStatusMap: filename provided --> use parameter file: " << mFileNameStatusMap.c_str() << endm;
236 paramFile.open( mFileNameStatusMap.c_str() );
238 if( !paramFile.is_open() ) {
239 LOG_ERROR <<
"unable to get the 'etofStatusMap' parameters from file --> file does not exist" << endm;
243 std::vector< int > status;
245 while( paramFile >> temp ) {
246 status.push_back( temp );
251 if( status.size() != eTofConst::nChannelsInSystem ) {
252 LOG_ERROR <<
"parameter file for 'etofStatusMap' has not the right amount of entries: ";
253 LOG_ERROR << status.size() <<
" instead of " << eTofConst::nChannelsInSystem <<
" !!!!" << endm;
257 for(
size_t i=0; i<eTofConst::nChannelsInSystem; i++ ) {
258 mStatus[ channelToKey( i ) ] = status[ i ];
262 for(
const auto& kv : mStatus ) {
263 LOG_DEBUG <<
"channel key: " << kv.first <<
" --> status = " << kv.second << endm;
269 mTimingWindow.clear();
270 mPulserWindow.clear();
272 if( mFileNameTimingWindow.empty() ) {
273 LOG_INFO <<
"etofTimingWindow: no filename provided --> load database table" << endm;
275 dbDataSet = GetDataBase(
"Calibrations/etof/etofTimingWindow" );
277 St_etofTimingWindow* etofTimingWindow =
static_cast< St_etofTimingWindow*
> ( dbDataSet->
Find(
"etofTimingWindow" ) );
278 if( !etofTimingWindow ) {
279 LOG_ERROR <<
"unable to get the timing window from the database" << endm;
283 etofTimingWindow_st* timingWindowTable = etofTimingWindow->GetTable();
285 for(
size_t i=0; i<12 ; i++ ) {
286 int key = timingWindowTable->afckAddress[ i ];
288 mTimingWindow[ key ] = std::make_pair( timingWindowTable->timingMin[ i ], timingWindowTable->timingMax[ i ] );
289 mPulserWindow[ key ] = std::make_pair( timingWindowTable->pulserMin[ i ], timingWindowTable->pulserMax[ i ] );
291 mPulserPeakTime = timingWindowTable->pulserPeak[ i ];
296 LOG_INFO <<
"etofTimingWindow: filename provided --> use parameter file: " << mFileNameTimingWindow.c_str() << endm;
298 paramFile.open( mFileNameTimingWindow.c_str() );
300 if( !paramFile.is_open() ) {
301 LOG_ERROR <<
"unable to get the 'etofTimingWindow' parameters from file --> file does not exist" << endm;
306 std::vector< float > times;
310 while( std::getline( paramFile, temp ) ) {
311 std::istringstream istring( temp );
314 istring >>std::hex >> address >> std::dec;
316 while( istring >> time ) {
317 times.push_back( time );
320 if( times.size() != 6 ) {
321 LOG_ERROR <<
"parameter file for 'etofTimingWindow' has not the right amount of entries: ";
322 LOG_ERROR <<
"times for address (0x" << std::hex << address << std::dec <<
") has " << times.size() <<
" instead of " << 6 <<
" !!!!" << endm;
327 mTimingWindow[ address ] = std::make_pair( times.at( 0 ), times.at( 1 ) );
328 mPulserWindow[ address ] = std::make_pair( times.at( 3 ), times.at( 4 ) );
330 mPulserPeakTime = times.at( 5 );
335 if( mTimingWindow.size() > 12 ) {
336 LOG_ERROR <<
" too many entries in mTimingWindow map ...." << endm;
339 if( mPulserWindow.size() > 12 ) {
340 LOG_ERROR <<
" too many entries in mPulserWindow map ...." << endm;
345 for(
const auto& kv : mTimingWindow ) {
346 LOG_DEBUG <<
"AFCK address: 0x" << std::hex << kv.first << std::dec <<
" --> timing window from " << kv.second.first <<
" to " << kv.second.second <<
" ns" << endm;
348 for(
const auto& kv : mPulserWindow ) {
349 LOG_DEBUG <<
"AFCK address: 0x" << std::hex << kv.first << std::dec <<
" --> pulser window from " << kv.second.first <<
" to " << kv.second.second <<
" ns" << endm;
351 LOG_DEBUG <<
"pulser time peak at " << mPulserPeakTime <<
" ns" << endm;
355 if( mFileNameCalibParam.empty() ) {
356 LOG_INFO <<
"etofCalibParam: no filename provided --> load database table" << endm;
358 dbDataSet = GetDataBase(
"Calibrations/etof/etofCalibParam" );
360 St_etofCalibParam* etofCalibParam =
static_cast< St_etofCalibParam*
> ( dbDataSet->
Find(
"etofCalibParam" ) );
361 if( !etofCalibParam ) {
362 LOG_ERROR <<
"unable to get the calibration params from the database" << endm;
366 etofCalibParam_st* calibParamTable = etofCalibParam->GetTable();
368 mGet4TotBinWidthNs = calibParamTable->get4TotBinWidthNs;
369 mMinDigisPerSlewBin = calibParamTable->minDigisInSlewBin;
372 if( mReferencePulserIndex == 0 ) {
373 mReferencePulserIndex = calibParamTable->referencePulserIndex;
376 LOG_INFO <<
"--- reference pulser index is set manually ---" << endm;
380 LOG_INFO <<
"etofCalibParam: filename provided --> use parameter file: " << mFileNameCalibParam.c_str() << endm;
382 paramFile.open( mFileNameCalibParam.c_str() );
384 if( !paramFile.is_open() ) {
385 LOG_ERROR <<
"unable to get the 'etofCalibParam' parameters from file --> file does not exist" << endm;
389 std::vector< float > param;
391 while( paramFile >> temp ) {
392 param.push_back( temp );
397 if( param.size() != 3 ) {
398 LOG_ERROR <<
"parameter file for 'etofCalibParam' has not the right amount of entries: ";
399 LOG_ERROR << param.size() <<
" instead of 3 !!!!" << endm;
403 if( param.at( 0 ) > 0. ) {
404 mGet4TotBinWidthNs = param.at( 0 );
406 if( param.at( 1 ) > 0 ) {
407 mMinDigisPerSlewBin = param.at( 1 );
411 if( param.at( 2 ) > 0 && mReferencePulserIndex == 0 ) {
412 mReferencePulserIndex = param.at( 2 );
415 LOG_INFO <<
"--- reference pulser index is set manually ---" << endm;
419 LOG_INFO <<
" Get4 TOT bin width to ns conversion factor: " << mGet4TotBinWidthNs << endm;
420 LOG_INFO <<
" minimal number of digis required in slewing bin: " << mMinDigisPerSlewBin << endm;
421 LOG_INFO <<
" reference pulser index: " << mReferencePulserIndex << endm;
427 mSignalVelocity.clear();
429 if( mFileNameSignalVelocity.empty() ) {
430 LOG_INFO <<
"etofSignalVelocity: no filename provided --> load database table" << endm;
432 dbDataSet = GetDataBase(
"Calibrations/etof/etofSignalVelocity" );
434 St_etofSignalVelocity* etofSignalVelocity =
static_cast< St_etofSignalVelocity*
> ( dbDataSet->
Find(
"etofSignalVelocity" ) );
435 if( !etofSignalVelocity ) {
436 LOG_ERROR <<
"unable to get the signal velocity from the database" << endm;
440 etofSignalVelocity_st* velocityTable = etofSignalVelocity->GetTable();
442 for(
size_t i=0; i<eTofConst::nCountersInSystem; i++ ) {
443 if( velocityTable->signalVelocity[ i ] > 0 ) {
444 mSignalVelocity[ detectorToKey( i ) ] = velocityTable->signalVelocity[ i ];
449 LOG_INFO <<
"etofSignalVelocity: filename provided --> use parameter file: " << mFileNameSignalVelocity.c_str() << endm;
451 paramFile.open( mFileNameSignalVelocity.c_str() );
453 if( !paramFile.is_open() ) {
454 LOG_ERROR <<
"unable to get the 'etofSignalVelocity' parameters from file --> file does not exist" << endm;
458 std::vector< float > velocity;
460 while( paramFile >> temp ) {
461 velocity.push_back( temp );
466 if( velocity.size() != eTofConst::nCountersInSystem ) {
467 LOG_ERROR <<
"parameter file for 'etofSignalVelocity' has not the right amount of entries: ";
468 LOG_ERROR << velocity.size() <<
" instead of " << eTofConst::nCountersInSystem <<
" !!!!" << endm;
472 for(
size_t i=0; i<eTofConst::nCountersInSystem; i++ ) {
473 if( velocity.at( i ) > 0 ) {
474 mSignalVelocity[ detectorToKey( i ) ] = velocity.at( i );
479 for(
const auto& kv : mSignalVelocity ) {
480 LOG_DEBUG <<
"counter key: " << kv.first <<
" --> signal velocity = " << kv.second <<
" cm / ns" << endm;
486 mDigiTotCorr.clear();
487 mDigiTimeCorr.clear();
488 mDigiSlewCorr.clear();
490 if( mFileNameCalibHistograms.empty() ) {
494 LOG_INFO <<
"etofDigiTotCorr: no filename provided --> load database table" << endm;
496 dbDataSet = GetDataBase(
"Calibrations/etof/etofDigiTotCorr" );
498 St_etofDigiTotCorr* etofDigiTotCorr =
static_cast< St_etofDigiTotCorr*
> ( dbDataSet->
Find(
"etofDigiTotCorr" ) );
499 if( !etofDigiTotCorr ) {
500 LOG_ERROR <<
"unable to get the digi tot correction from the database" << endm;
504 etofDigiTotCorr_st* digiTotCorrTable = etofDigiTotCorr->GetTable();
506 for(
size_t i=0; i<eTofConst::nChannelsInSystem; i++ ) {
507 unsigned int key = channelToKey( i );
508 unsigned int detector = key / 1000;
510 if( mDigiTotCorr.count( detector ) == 0 ) {
511 TString name = Form(
"digiTotCorr_%d", detector );
512 mDigiTotCorr[ detector ] =
new TH1F( name, name, 2 * eTofConst::nStrips, 0, 2 * eTofConst::nStrips );
515 unsigned int strip = ( key % 1000 ) / 10;
516 unsigned int side = key % 10;
519 LOG_DEBUG << i <<
" " << detector <<
" " << strip <<
" " << side <<
" " << digiTotCorrTable->totCorr[ i ] << endm;
522 mDigiTotCorr.at( detector )->SetBinContent( strip + eTofConst::nStrips * ( side - 1 ) , digiTotCorrTable->totCorr[ i ] );
525 for(
auto& kv : mDigiTotCorr ) {
526 kv.second->SetDirectory( 0 );
533 LOG_INFO <<
"etofDigiTimeCorr: no filename provided --> load database table" << endm;
536 dbDataSet = GetDataBase(
"Calibrations/etof/etofDigiTimeCorr" );
538 St_etofDigiTimeCorr* etofDigiTimeCorr =
static_cast< St_etofDigiTimeCorr*
> ( dbDataSet->
Find(
"etofDigiTimeCorr" ) );
539 if( !etofDigiTimeCorr ) {
540 LOG_ERROR <<
"unable to get the digi time correction from the database" << endm;
544 etofDigiTimeCorr_st* digiTimeCorrTable = etofDigiTimeCorr->GetTable();
546 for(
size_t i=0; i<eTofConst::nChannelsInSystem; i++ ) {
547 unsigned int key = channelToKey( i );
548 unsigned int detector = key / 1000;
550 if( mDigiTimeCorr.count( detector ) == 0 ) {
551 TString name = Form(
"digiTimeCorr_%d", detector );
552 mDigiTimeCorr[ detector ] =
new TH1F( name, name, 2 * eTofConst::nStrips, 0, 2 * eTofConst::nStrips );
555 unsigned int strip = ( key % 1000 ) / 10;
556 unsigned int side = key % 10;
559 LOG_DEBUG << i <<
" " << detector <<
" " << strip <<
" " << side <<
" " << digiTimeCorrTable->timeCorr[ i ] << endm;
562 mDigiTimeCorr.at( detector )->SetBinContent( strip + eTofConst::nStrips * ( side - 1 ) , digiTimeCorrTable->timeCorr[ i ] );
565 for(
auto& kv : mDigiTimeCorr ) {
566 kv.second->SetDirectory( 0 );
573 LOG_INFO <<
"etofDigiSlewCorr: no filename provided --> load database table" << endm;
575 dbDataSet = GetDataBase(
"Calibrations/etof/etofDigiSlewCorr" );
577 St_etofDigiSlewCorr* etofDigiSlewCorr =
static_cast< St_etofDigiSlewCorr*
> ( dbDataSet->
Find(
"etofDigiSlewCorr" ) );
578 if( !etofDigiSlewCorr ) {
579 LOG_ERROR <<
"unable to get the digi slewing correction from the database" << endm;
583 etofDigiSlewCorr_st* digiSlewCorrTable = etofDigiSlewCorr->GetTable();
585 for(
size_t i=0; i<eTofConst::nChannelsInSystem; i++ ) {
586 unsigned int key = channelToKey( i );
589 std::vector< float > binEdges;
590 std::vector< float > corr;
592 binEdges.push_back( 0. );
593 for(
size_t j=0; j<etofSlewing::nTotBins; j++ ) {
594 int tableIndex = i * etofSlewing::nTotBins + j;
596 if( i != digiSlewCorrTable->channelId[ tableIndex ] ) {
597 LOG_ERROR <<
"something went wrong in reading the database table for eTOF slewing corrections" << endm;
601 binEdges.push_back( digiSlewCorrTable->upperTotEdge[ tableIndex ] * etofSlewing::conversionFactor );
602 corr.push_back( digiSlewCorrTable->corr[ tableIndex ] * etofSlewing::conversionFactor );
606 if( mDigiSlewCorr.count( key ) == 0 ) {
607 TString name = Form(
"digiSlewCorr_%d", key );
608 mDigiSlewCorr[ key ] =
new TProfile( name, name, etofSlewing::nTotBins, binEdges.data() );
612 for(
size_t j=0; j<etofSlewing::nTotBins; j++ ) {
613 float totCenter = mDigiSlewCorr.at( key )->GetBinCenter( j+1 );
614 mDigiSlewCorr.at( key )->Fill( totCenter, corr.at( j ), mMinDigisPerSlewBin + 1 );
618 for(
auto& kv : mDigiSlewCorr ) {
619 kv.second->SetDirectory( 0 );
623 LOG_INFO <<
"etof calibration histograms: filename provided --> use parameter file: " << mFileNameCalibHistograms.c_str() << endm;
625 if( !isFileExisting( mFileNameCalibHistograms ) ) {
626 LOG_ERROR <<
"unable to get the 'etofDigiTotCorr', 'etofDigiTimeCorr', 'etofDigiSlewCorr' parameters from file --> file does not exist" << endm;
629 TFile* histFile =
new TFile( mFileNameCalibHistograms.c_str(),
"READ" );
632 LOG_ERROR <<
"No calibration file found: " << mFileNameCalibHistograms.c_str() << endm;
633 LOG_INFO <<
"setting all parameters to default" << endm;
634 }
else if( histFile->IsZombie() ){
635 LOG_ERROR <<
"Zombie calibration file: " << mFileNameCalibHistograms.c_str() << endm;
636 LOG_INFO <<
"stopping execution" << endm;
640 TFile* histOffsetFile =
new TFile( mFileNameOffsetHistograms.c_str(),
"READ" );
642 if( !histOffsetFile ) {
643 LOG_INFO <<
"No offset file found: " << mFileNameOffsetHistograms.c_str() << endm;
644 LOG_INFO <<
"setting all parameters to default" << endm;
645 }
else if( histOffsetFile->IsZombie() ) {
646 LOG_ERROR <<
"Zombie offset file: " << mFileNameOffsetHistograms.c_str() << endm;
647 LOG_INFO <<
"stopping execution" << endm;
650 LOG_INFO <<
"Successfully opened RunOffset file "<< mFileNameOffsetHistograms.c_str() << endm;
655 TString hPosOffsetName = Form(
"calib_Run%d_PosOffsets_pfx", runnumber );
656 TString hT0OffsetName = Form(
"calib_Run%d_T0Offsets_pfx", runnumber );
657 TProfile* hPosOffsetProfile =
nullptr;
658 TProfile* hT0OffsetProfile =
nullptr;
660 if( histOffsetFile && !(histOffsetFile->IsZombie()) ){
661 LOG_INFO <<
"Getting run offset histograms Pos: "<< hPosOffsetName <<
" T0: "<< hT0OffsetName << endm;
662 hPosOffsetProfile = ( TProfile* ) histOffsetFile->Get( hPosOffsetName );
663 hT0OffsetProfile = ( TProfile* ) histOffsetFile->Get( hT0OffsetName );
666 for(
unsigned int sector = eTofConst::sectorStart; sector <= eTofConst::sectorStop; sector++ ) {
668 if( mRunYear == 2018 && sector != 18 )
continue;
670 for(
unsigned int zPlane = eTofConst::zPlaneStart; zPlane <= eTofConst::zPlaneStop; zPlane++ ) {
671 for(
unsigned int counter = eTofConst::counterStart; counter <= eTofConst::counterStop; counter++ ) {
673 unsigned int key = sector * 100 + zPlane * 10 + counter;
675 LOG_DEBUG <<
"detectorId key: " << sector <<
" " << zPlane <<
" " << counter << endm;
684 if( mDigiTotCorr.count( key ) == 0 ) {
685 TString name = Form(
"digiTotCorr_%d", key );
686 mDigiTotCorr[ key ] =
new TH1F( name, name, 2 * eTofConst::nStrips, 0, 2 * eTofConst::nStrips );
689 hname = Form(
"calib_Sector%02d_ZPlane%d_Det%d_TotDigi_pfx", sector, zPlane, counter );
690 hProfile = ( TProfile* ) histFile->Get( hname );
693 for(
size_t i=1; i<=2 * eTofConst::nStrips; i++ ) {
694 if( hProfile->GetBinContent( i ) != 0 ) {
695 mDigiTotCorr.at( key )->SetBinContent( i , hProfile->GetBinContent( i ) );
698 mDigiTotCorr.at( key )->SetBinContent( i , 1. );
702 if( mDigiTotCorr.at( key )->GetBinContent( i ) < 0.05 || mDigiTotCorr.at( key )->GetBinContent( i ) > 10 ) {
703 mDigiTotCorr.at( key )->SetBinContent( i , 1. );
708 if( isFileExisting( mFileNameCalibHistograms ) ) {
709 LOG_WARN <<
"unable to find histogram: " << hname << endm;
712 for(
size_t i=1; i<=2 * eTofConst::nStrips; i++ ) {
713 mDigiTotCorr.at( key )->SetBinContent( i , 1. );
717 LOG_DEBUG <<
"histogram " << mDigiTotCorr.at( key )->GetName() <<
" filled" << endm;
723 if( mDigiTimeCorr.count( key ) == 0 ) {
724 TString name = Form(
"digiTimeCorr_%d", key );
725 mDigiTimeCorr[ key ] =
new TH1F( name, name, 2 * eTofConst::nStrips, 0, 2 * eTofConst::nStrips );
729 hname = Form(
"calib_Sector%02d_ZPlane%d_Det%d_Pos_pfx", sector, zPlane, counter );
730 hProfile = ( TProfile* ) histFile->Get( hname );
732 double dRunOffset = 0;
733 if (hPosOffsetProfile) {
734 int mCounterBin = hPosOffsetProfile->FindBin( 9*(sector -13) + 3 * (zPlane - 1) + counter );
735 dRunOffset = hPosOffsetProfile->GetBinContent( mCounterBin );
736 LOG_DEBUG <<
"setting run position offset to "<< dRunOffset<<
" for counter "<< ( 9*(sector -13) + 3 * (zPlane - 1) + counter ) << endm;
738 LOG_INFO <<
"position offset histogram "<<hPosOffsetName <<
" not found" << endm;
742 for(
size_t i=1; i<= eTofConst::nStrips; i++ ) {
743 mDigiTimeCorr.at( key )->AddBinContent( i , -1. * ( hProfile->GetBinContent( i ) + dRunOffset ) / mSignalVelocity.at( key ) );
744 mDigiTimeCorr.at( key )->AddBinContent( eTofConst::nStrips + i , ( hProfile->GetBinContent( i ) + dRunOffset ) / mSignalVelocity.at( key ) );
748 if( isFileExisting( mFileNameCalibHistograms ) ) {
749 LOG_WARN <<
"unable to find histogram: " << hname << endm;
754 hname = Form(
"calib_Sector%02d_ZPlane%d_Det%d_T0corr", sector, zPlane, counter );
755 hProfile = ( TProfile* ) histFile->Get( hname );
758 if (hT0OffsetProfile) {
759 int mCounterBin = hT0OffsetProfile->FindBin( 9*(sector -13) + 3 * (zPlane - 1) + counter );
760 dRunOffset = hT0OffsetProfile->GetBinContent( mCounterBin );
762 LOG_DEBUG <<
"setting run time offset to "<< dRunOffset<<
" for counter "<< ( 9*(sector -13) + 3 * (zPlane - 1) + counter ) << endm;
764 if( hProfile && hProfile->GetNbinsX() == 1 ) {
765 LOG_DEBUG <<
"T0 histogram with 1 bin: " << key << endm;
766 for(
size_t i=1; i<=2 * eTofConst::nStrips; i++ ) {
767 mDigiTimeCorr.at( key )->AddBinContent( i , hProfile->GetBinContent( 1 ) + dRunOffset );
770 else if( hProfile && hProfile->GetNbinsX() == eTofConst::nStrips ) {
771 LOG_DEBUG <<
"T0 histogram with 32 bins: " << key << endm;
772 for(
size_t i=1; i<= eTofConst::nStrips; i++ ) {
773 mDigiTimeCorr.at( key )->AddBinContent( i , hProfile->GetBinContent( i ) + dRunOffset );
774 mDigiTimeCorr.at( key )->AddBinContent( i + eTofConst::nStrips , hProfile->GetBinContent( i ) + dRunOffset );
778 if( isFileExisting( mFileNameCalibHistograms ) ) {
779 LOG_WARN <<
"unable to find histogram: " << hname << endm;
783 LOG_DEBUG <<
"histogram " << mDigiTimeCorr.at( key )->GetName() <<
" filled" << endm;
789 for(
unsigned int chan = 0; chan < eTofConst::nChannelsPerCounter; chan++ ) {
790 unsigned int strip = ( chan % 32 ) + 1;
791 unsigned int side = ( chan / 32 ) + 1;
792 unsigned int key = sector * 100000 + zPlane * 10000 + counter * 1000 + strip * 10 + side;
794 LOG_DEBUG <<
"channelId key: " << sector <<
" " << zPlane <<
" " << counter <<
" " << strip <<
" " << side << endm;
796 hname = Form(
"calib_Sector%02d_ZPlane%d_Det%d_Chan%d_Walk_pfx", sector, zPlane, counter, chan );
797 hProfile = ( TProfile* ) histFile->Get( hname );
801 LOG_DEBUG <<
"unable to find histogram: " << hname <<
"--> check detector-wise" << endm;
802 hname = Form(
"calib_Sector%02d_ZPlane%d_Det%d_AvWalk_pfx", sector, zPlane, counter );
803 hProfile = ( TProfile* ) histFile->Get( hname );
807 unsigned int nbins = hProfile->GetNbinsX();
809 if( mDigiSlewCorr.count( key ) == 0 ) {
811 std::vector< float > bins( nbins + 1 );
813 for(
size_t i=0; i<nbins; i++ ) {
814 bins.at( i ) = hProfile->GetXaxis()->GetBinLowEdge( i+1 );
816 bins.at( nbins ) = hProfile->GetXaxis()->GetBinUpEdge( nbins );
818 TString name = Form(
"digiSlewCorr_%d", key );
819 mDigiSlewCorr[ key ] =
new TProfile( name, name, nbins, bins.data() );
822 for(
size_t iTotBin=1; iTotBin<=nbins; iTotBin++ ) {
823 mDigiSlewCorr.at( key )->Fill( hProfile->GetBinCenter( iTotBin ), hProfile->GetBinContent( iTotBin ), hProfile->GetBinEntries( iTotBin ) );
827 LOG_DEBUG <<
"unable to find histogram: " << hname << endm;
835 for(
auto& kv : mDigiTotCorr ) {
836 kv.second->SetDirectory( 0 );
838 for(
auto& kv : mDigiTimeCorr ) {
839 kv.second->SetDirectory( 0 );
841 for(
auto& kv : mDigiSlewCorr ) {
842 kv.second->SetDirectory( 0 );
855 if( mFileNameResetTimeCorr.empty() ) {
856 LOG_INFO <<
"etofResetTimeCorr: no filename provided --> load database table" << endm;
858 dbDataSet = GetDataBase(
"Calibrations/etof/etofResetTimeCorr" );
860 St_etofResetTimeCorr* etofResetTimeCorr =
static_cast< St_etofResetTimeCorr*
> ( dbDataSet->
Find(
"etofResetTimeCorr" ) );
861 if( !etofResetTimeCorr ) {
862 LOG_ERROR <<
"unable to get the reset time correction from the database" << endm;
866 etofResetTimeCorr_st* resetTimeCorrTable = etofResetTimeCorr->GetTable();
868 mResetTimeCorr = resetTimeCorrTable->resetTimeOffset;
871 LOG_INFO <<
"etofResetTimeCorr: filename provided --> use parameter file: " << mFileNameResetTimeCorr.c_str() << endm;
873 paramFile.open( mFileNameResetTimeCorr.c_str() );
875 if( !paramFile.is_open() ) {
876 LOG_ERROR <<
"unable to get the 'etofResetTimeCorr' parameters from file --> file does not exist" << endm;
880 std::map< unsigned int, float > param;
883 while( paramFile >> temp ) {
885 param[ (
unsigned int ) temp ] = temp2;
890 for(
const auto& kv : param ) {
891 LOG_DEBUG <<
"run: " << kv.first <<
" --> reset time corr = " << kv.second <<
" ns" << endm;
894 if( param.count( runnumber ) ) {
895 mResetTimeCorr = param.at( runnumber );
899 LOG_INFO <<
"additional reset time offset correction: " << mResetTimeCorr <<
" ns" << endm;
904 mPulserPeakTot.clear();
906 if( mFileNamePulserTotPeak.empty() ) {
907 LOG_INFO <<
"etofPulserPeakTot: no filename provided --> load database table" << endm;
909 dbDataSet = GetDataBase(
"Calibrations/etof/etofPulserTotPeak" );
911 St_etofPulserTotPeak* etofPulserTotPeak =
static_cast< St_etofPulserTotPeak*
> ( dbDataSet->
Find(
"etofPulserTotPeak" ) );
912 if( !etofPulserTotPeak ) {
913 LOG_ERROR <<
"unable to get the pulser tot peak parameters from the database" << endm;
917 etofPulserTotPeak_st* pulserTotTable = etofPulserTotPeak->GetTable();
919 for(
size_t i=0; i<eTofConst::nCountersInSystem * 2; i++ ) {
920 if( pulserTotTable->pulserTot[ i ] > 0 ) {
921 mPulserPeakTot[ sideToKey( i ) ] = ( int ) pulserTotTable->pulserTot[ i ];
926 LOG_INFO <<
"etofPulserPeakTot: filename provided --> use parameter file: " << mFileNamePulserTotPeak.c_str() << endm;
928 paramFile.open( mFileNamePulserTotPeak.c_str() );
930 if( !paramFile.is_open() ) {
931 LOG_ERROR <<
"unable to get the 'etofPulserTotPeak' parameters from file --> file does not exist" << endm;
935 std::vector< float > param;
937 while( paramFile >> temp ) {
938 param.push_back( temp );
943 if( param.size() != eTofConst::nCountersInSystem * 2 ) {
944 LOG_ERROR <<
"parameter file for 'etofPulserTotPeak' has not the right amount of entries: ";
945 LOG_ERROR << param.size() <<
" instead of " << eTofConst::nCountersInSystem * 2 <<
" !!!!" << endm;
949 for(
size_t i=0; i<eTofConst::nCountersInSystem * 2; i++ ) {
950 if( param.at( i ) > 0 ) {
951 mPulserPeakTot[ sideToKey( i ) ] = param.at( i );
956 for(
const auto& kv : mPulserPeakTot ) {
957 LOG_DEBUG <<
"side key: " << kv.first <<
" --> pulser peak tot = " << kv.second <<
" (bin)" << endm;
963 mPulserTimeDiffGbtx.clear();
965 if( mFileNamePulserTimeDiffGbtx.empty() ) {
966 LOG_INFO <<
"etofPulserTimeDiffGbtx: no filename provided --> load database table" << endm;
968 dbDataSet = GetDataBase(
"Calibrations/etof/etofPulserTimeDiffGbtx" );
970 St_etofPulserTimeDiffGbtx* etofPulserTimeDiffGbtx =
static_cast< St_etofPulserTimeDiffGbtx*
> ( dbDataSet->
Find(
"etofPulserTimeDiffGbtx" ) );
971 if( !etofPulserTimeDiffGbtx ) {
972 LOG_ERROR <<
"unable to get the pulser time difference parameters from the database" << endm;
976 etofPulserTimeDiffGbtx_st* pulserTimeTable = etofPulserTimeDiffGbtx->GetTable();
978 for(
size_t i=0; i<eTofConst::nModules * eTofConst::nSides; i++ ) {
979 unsigned int sector = eTofConst::sectorStart + i / ( eTofConst::nPlanes * eTofConst::nSides );
980 unsigned int zPlane = eTofConst::zPlaneStart + ( i % ( eTofConst::nPlanes * eTofConst::nSides ) ) / eTofConst::nSides;
981 unsigned int side = eTofConst::counterStart + i % eTofConst::nSides;
983 mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 1 * 10 + side ] = 0.;
984 mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 2 * 10 + side ] = pulserTimeTable->pulserTimeDiffGbtx[ i * 2 + 0 ];
985 mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 3 * 10 + side ] = pulserTimeTable->pulserTimeDiffGbtx[ i * 2 + 1 ];
989 LOG_INFO <<
"etofPulserTimeDiff: filename provided --> use parameter file: " << mFileNamePulserTimeDiffGbtx.c_str() << endm;
991 paramFile.open( mFileNamePulserTimeDiffGbtx.c_str() );
993 if( !paramFile.is_open() ) {
994 LOG_ERROR <<
"unable to get the 'etotPulserTimeDiffGbtc' parameters from file --> file does not exist" << endm;
998 std::vector< float > param;
1000 while( paramFile >> temp ) {
1001 param.push_back( temp );
1006 if( param.size() != eTofConst::nModules * eTofConst::nSides * 2 ) {
1007 LOG_ERROR <<
"parameter file for 'etofPulserTimeDiffGbtx' has not the right amount of entries: ";
1008 LOG_ERROR << param.size() <<
" instead of " << eTofConst::nModules * eTofConst::nSides <<
" !!!!" << endm;
1012 for(
size_t i=0; i<eTofConst::nModules * eTofConst::nSides; i++ ) {
1013 unsigned int sector = eTofConst::sectorStart + i / ( eTofConst::nPlanes * eTofConst::nSides );
1014 unsigned int zPlane = eTofConst::zPlaneStart + ( i % ( eTofConst::nPlanes * eTofConst::nSides ) ) / eTofConst::nSides;
1015 unsigned int side = eTofConst::counterStart + i % eTofConst::nSides;
1017 mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 1 * 10 + side ] = 0.;
1018 mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 2 * 10 + side ] = param.at( i * 2 + 0 );
1019 mPulserTimeDiffGbtx[ sector * 1000 + zPlane * 100 + 3 * 10 + side ] = param.at( i * 2 + 1 );
1025 bool allZero =
true;
1027 for(
const auto& kv : mPulserTimeDiffGbtx ) {
1028 if( fabs( kv.second ) > 1e-4 ) allZero =
false;
1031 LOG_INFO <<
"side key: " << kv.first <<
" --> pulser time diff inside Gbtx ( to counter 1 ) = " << kv.second << endm;
1036 mUsePulserGbtxDiff =
false;
1037 LOG_INFO <<
"the use of pulser relations inside a Gbtx is turned off" << endm;
1041 mHistograms.at(
"pulserDigiTimeDiff_GbtxCorrProfMod" )->Reset();
1052 StETofCalibMaker::FinishRun( Int_t runnumber )
1060 for(
auto kv = mDigiTotCorr.begin(); kv != mDigiTotCorr.end(); kv++ ) {
1062 kv->second =
nullptr;
1064 mDigiTotCorr.clear();
1066 for(
auto kv = mDigiTimeCorr.begin(); kv != mDigiTimeCorr.end(); kv++ ) {
1068 kv->second =
nullptr;
1070 mDigiTimeCorr.clear();
1072 for(
auto kv = mDigiSlewCorr.begin(); kv != mDigiSlewCorr.end(); kv++ ) {
1074 kv->second =
nullptr;
1076 mDigiSlewCorr.clear();
1078 mPulserPeakTot.clear();
1079 mPulserTimeDiff.clear();
1081 mJumpingPulsers.clear();
1092 LOG_INFO <<
"Finish() - writing *.etofCalib.root ..." << endm;
1104 LOG_DEBUG <<
"StETofCalibMaker::Make(): starting ..." << endm;
1106 mEvent = (
StEvent* ) GetInputDS(
"StEvent" );
1110 LOG_DEBUG <<
"Make(): running on StEvent" << endm;
1114 if( !etofCollection ) {
1115 LOG_WARN <<
"Make() - Found StEvent data structure, but no eTOF collection. Try MuDst processing instead" << endm;
1116 mMuDst = (
StMuDst* ) GetInputDS(
"MuDst" );
1119 LOG_DEBUG <<
"Make() - running on MuDsts" << endm;
1132 LOG_DEBUG <<
"Make(): no StEvent found" << endm;
1134 mMuDst = (
StMuDst* ) GetInputDS(
"MuDst" );
1137 LOG_DEBUG <<
"Make(): running on MuDsts" << endm;
1144 LOG_WARN <<
"Make() - no StMuDst or StEvent" << endm;
1158 if( !etofCollection ) {
1159 LOG_WARN <<
"processStEvent() - no etof collection" << endm;
1163 if( !etofCollection->digisPresent() ) {
1164 LOG_WARN <<
"processStEvent() - no digis present" << endm;
1168 if( !etofCollection->etofHeader() ) {
1169 LOG_WARN <<
"processStEvent() - no header" << endm;
1175 StETofHeader* etofHeader = etofCollection->etofHeader();
1176 StSPtrVecETofDigi& etofDigis = etofCollection->etofDigis();
1193 size_t nDigis = etofDigis.size();
1195 LOG_INFO <<
"processStEvent() - # fired eTOF digis : " << nDigis << endm;
1198 mTriggerTime = triggerTime( etofHeader );
1199 mResetTime = fmod( resetTime( etofHeader ), eTofConst::bTofClockCycle );
1201 std::map< unsigned int, std::vector< unsigned int > > pulserCandMap;
1204 for(
size_t i=0; i<nDigis; i++ ) {
1209 LOG_WARN <<
"No digi found" << endm;
1214 resetToRaw( aDigi );
1219 applyMapping( aDigi );
1222 if( mRunYear != 2018 ) {
1223 flagPulserDigis( aDigi, i, pulserCandMap );
1229 LOG_INFO <<
"size of pulserCandMap: " << pulserCandMap.size() << endm;
1231 calculatePulserOffsets( pulserCandMap );
1234 TClass* headerClass = etofHeader->IsA();
1235 if( headerClass->GetClassVersion() > 2 ){
1237 std::vector<bool> goodEventFlagVec;
1238 std::vector<bool> hasPulsersVec;
1241 for(
unsigned int iCounter = 0; iCounter < 108; iCounter++){
1242 if ( !(mNPulsersCounter.count(iCounter) ) ){
1243 hasPulsersVec.push_back(
false);
1245 hasPulsersVec.push_back(mNPulsersCounter[iCounter] == 2);
1248 if (hasPulsersVec.size() == 108){
1253 for(
unsigned int iGet4 = 0; iGet4 < 1728; iGet4++){
1257 if (goodEventFlagVec.size() == 1728){
1258 etofHeader->setGoodEventFlagVec(goodEventFlagVec);
1266 int nDuplicates = 0;
1268 for(
size_t i=0; i<nDigis; i++ ) {
1272 LOG_WARN <<
"No digi found" << endm;
1277 current.tot = aDigi->
rawTot();
1278 current.time = aDigi->
rawTime();
1281 auto it = std::find( mStuckFwDigi.begin(), mStuckFwDigi.end(), current );
1282 if( it != mStuckFwDigi.end() ) {
1284 LOG_INFO <<
"digi from stuck firmware (s" << aDigi->
sector() <<
" m" << aDigi->
zPlane() <<
" c" << aDigi->
counter() <<
") --> ignore" << endm;
1290 else if( current == prev ) {
1291 mStuckFwDigi.push_back( current );
1292 resetToRaw( mMuDst->
etofDigi( i-1 ) );
1304 applyCalibration( aDigi, etofHeader );
1307 if( mDoQA && nDuplicates > 0 ) {
1308 LOG_INFO <<
"*** # duplicate digis from stuck firmware: " << nDuplicates << endm;
1309 for(
const auto& v : mStuckFwDigi ) {
1310 LOG_INFO <<
"stuck channel with geomId: " << v.geomId <<
" " << v.tot <<
" " << v.time << endm;
1313 mStuckFwDigi.clear();
1321 LOG_DEBUG <<
"processMuDst(): starting ..." << endm;
1324 if( !mMuDst->
etofArray( muETofDigi ) ) {
1325 LOG_WARN <<
"processMuDst() - no digi array" << endm;
1329 if( !mMuDst->numberOfETofDigi() ) {
1330 LOG_WARN <<
"processMuDst() - no digis present" << endm;
1334 if( !mMuDst->
etofArray( muETofHeader ) ) {
1335 LOG_WARN <<
"processMuDst() - no digi array" << endm;
1340 mNPulsersCounter.clear();
1359 size_t nDigis = mMuDst->numberOfETofDigi();
1362 mTriggerTime = triggerTime( (
StETofHeader* ) etofHeader );
1363 mResetTime = fmod( resetTime( (
StETofHeader* ) etofHeader ), eTofConst::bTofClockCycle );
1364 std::map< unsigned int, std::vector< unsigned int >> pulserCandMap;
1368 for(
size_t i=0; i<nDigis; i++ ) {
1373 LOG_WARN <<
"No digi found" << endm;
1378 resetToRaw( aDigi );
1384 applyMapping( aDigi );
1390 if( mRunYear != 2018 ) {
1391 flagPulserDigis( aDigi, i, pulserCandMap );
1398 calculatePulserOffsets( pulserCandMap );
1402 TClass* headerClass = etofHeader->IsA();
1403 if( headerClass->GetClassVersion() > 2 ){
1405 std::vector<bool> goodEventFlagVec;
1406 std::vector<bool> hasPulsersVec;
1409 for(
unsigned int iCounter = 0; iCounter < 108; iCounter++){
1410 if ( !(mNPulsersCounter.count(iCounter) ) ){
1411 hasPulsersVec.push_back(
false);
1413 hasPulsersVec.push_back(mNPulsersCounter[iCounter] == 2);
1417 if (hasPulsersVec.size() == 108){
1418 etofHeader->setHasPulsersVec(hasPulsersVec);
1422 for(
unsigned int iGet4 = 0; iGet4 < 1728; iGet4++){
1426 if (goodEventFlagVec.size() == 1728){
1427 etofHeader->setGoodEventFlagVec(goodEventFlagVec);
1434 int nDuplicates = 0;
1436 for(
size_t i=0; i<nDigis; i++ ) {
1440 LOG_WARN <<
"No digi found" << endm;
1445 current.tot = aDigi->
rawTot();
1446 current.time = aDigi->
rawTime();
1449 auto it = std::find( mStuckFwDigi.begin(), mStuckFwDigi.end(), current );
1450 if( it != mStuckFwDigi.end() ) {
1452 LOG_INFO <<
"digi from stuck firmware (s" << aDigi->
sector() <<
" m" << aDigi->
zPlane() <<
" c" << aDigi->
counter() <<
") --> ignore" << endm;
1458 else if( current == prev ) {
1459 mStuckFwDigi.push_back( current );
1460 resetToRaw( mMuDst->
etofDigi( i-1 ) );
1471 applyCalibration( aDigi, etofHeader );
1474 if( mDoQA && nDuplicates > 0 ) {
1475 LOG_INFO <<
"*** # duplicate digis from stuck firmware: " << nDuplicates << endm;
1476 for(
const auto& v : mStuckFwDigi ) {
1477 LOG_INFO <<
"stuck channel with geomId: " << v.geomId <<
" " << v.tot <<
" " << v.time << endm;
1480 mStuckFwDigi.clear();
1487 StETofCalibMaker::isFileExisting(
const std::string fileName )
1489 std::ifstream infile( fileName );
1490 return infile.good();
1499 StETofCalibMaker::resetToRaw(
StETofDigi* aDigi )
1501 if( !aDigi )
return;
1504 aDigi->setCalibTime( 0. );
1505 aDigi->setCalibTot( -1. );
1507 aDigi->setAssociatedHit(
nullptr );
1516 StETofCalibMaker::applyMapping(
StETofDigi* aDigi )
1518 std::vector< unsigned int > geomVec;
1520 unsigned int rocId = aDigi->
rocId();
1521 unsigned int get4Id = aDigi->
get4Id();
1522 unsigned int elChan = aDigi->
elChan();
1524 mHwMap->mapToGeom( rocId, get4Id, elChan, geomVec );
1526 if( geomVec.size() < 5 ) {
1528 LOG_ERROR <<
"geometry vector has wrong size !!! --> skip digi" << endm;
1533 unsigned int sector = geomVec.at( 0 );
1534 unsigned int zplane = geomVec.at( 1 );
1535 unsigned int counter = geomVec.at( 2 );
1536 unsigned int strip = geomVec.at( 3 );
1537 unsigned int side = geomVec.at( 4 );
1539 if( mDebug && ( sector == 0 || zplane == 0 || counter == 0 || strip == 0 || side == 0 ) ) {
1540 LOG_ERROR <<
"geometry vector has entries equal to zero !!! --> skip digi" << endm;
1543 aDigi->
setGeoAddress( sector, zplane, counter, strip, side );
1547 LOG_DEBUG <<
"sector, zplane, counter, strip, side: " << aDigi->
sector() <<
", ";
1548 LOG_DEBUG << aDigi->
zPlane() <<
", " << aDigi->
counter() <<
", ";
1549 LOG_DEBUG << aDigi->
strip() <<
", " << aDigi->
side() << endm;
1551 LOG_DEBUG <<
"continuous module number: " << mHwMap->module( aDigi->
sector(), aDigi->
zPlane() ) << endm;
1561 StETofCalibMaker::flagPulserDigis(
StETofDigi* aDigi,
unsigned int index, std::map<
unsigned int, std::vector< unsigned int > >& pulserDigiMap )
1563 bool isPulserCand =
false;
1569 if( ( aDigi->
strip() == 1 && aDigi->
side() == 1 ) || ( aDigi->
strip() == 32 && aDigi->
side() == 2 ) ) {
1570 float timeToTrigger = aDigi->
rawTime() - mTriggerTime;
1573 float totToPeak = aDigi->
rawTot() - mPulserPeakTot.at( key );
1574 float totToHalfPeak = aDigi->
rawTot() - mPulserPeakTot.at( key ) * 0.5;
1577 if( timeToTrigger > mPulserWindow.at( aDigi->
rocId() ).first && timeToTrigger < mPulserWindow.at( aDigi->
rocId() ).second ) {
1578 if( fabs( totToPeak ) < 25 || fabs( totToHalfPeak ) < 10 ) {
1579 isPulserCand =
true;
1585 if( isPulserCand ) {
1586 pulserDigiMap[ key ].push_back( index );
1599 StETofCalibMaker::calculatePulserOffsets( std::map<
unsigned int, std::vector< unsigned int > >& pulserDigiMap )
1602 for(
auto it=pulserDigiMap.begin(); it!=pulserDigiMap.end(); it++ ) {
1603 LOG_DEBUG <<
"channel: " << it->first <<
" nCandidates: " << it->second.size() << endm;
1607 if( mReferencePulserIndex == 0 ) {
1609 LOG_INFO <<
"reference pulser index is 0 --> pulser correction is turned off" << endm;
1615 LOG_INFO <<
"reference pulser index: " << mReferencePulserIndex << endm;
1618 std::map< int, double > pulserTimes;
1619 mNPulsersCounter.clear();
1620 mPulserPresent.clear();
1623 for(
auto it=pulserDigiMap.begin(); it!=pulserDigiMap.end(); it++ ) {
1624 if( it->second.size() == 0 ) {
1627 int sideIndex = it->first;
1629 double bestDiff = 100000;
1632 for(
size_t j=0; j<it->second.size(); j++ ) {
1633 double pulserTime = 0.;
1634 double pulserTot = 0.;
1636 pulserTime = mMuDst->
etofDigi( it->second.at( j ) )->rawTime();
1637 pulserTot = mMuDst->
etofDigi( it->second.at( j ) )->rawTot();
1638 }
else if( mEvent ) {
1639 pulserTime = ( mEvent->etofCollection()->etofDigis() )[ it->second.at( j ) ]->
rawTime();
1640 pulserTot = ( mEvent->etofCollection()->etofDigis() )[ it->second.at( j ) ]->rawTot();
1642 double timeToTrigger = pulserTime - mTriggerTime;
1643 double totToPeak = pulserTot - mPulserPeakTot.at( sideIndex );
1645 if( mDebug && it->second.size() > 1 ) {
1646 LOG_INFO << it->second.size() <<
" pulsers @ " << sideIndex <<
" : timeToTrigger: " << timeToTrigger <<
" tot: " << pulserTot << endm;
1650 double currentDiff = fabs( timeToTrigger - mPulserPeakTime ) * 0.1 + fabs( totToPeak );
1651 if( currentDiff < bestDiff ) {
1652 bestDiff = currentDiff;
1657 if( mDebug && it->second.size() > 1 ) {
1658 LOG_INFO <<
" --> selected CAND-INDEX: " << candIndex << endm;
1661 double pulserTime = 0.;
1664 pulserTime = mMuDst->
etofDigi( it->second.at( candIndex ) )->rawTime();
1667 mMuDst->
etofDigi( it->second.at( candIndex ) )->setCalibTot( -999. );
1668 }
else if( mEvent ) {
1669 pulserTime = ( mEvent->etofCollection()->etofDigis() )[ it->second.at( candIndex ) ]->
rawTime();
1672 mEvent->etofCollection()->etofDigis() [ it->second.at( candIndex ) ]->setCalibTot( -999. );
1675 pulserTimes[ sideIndex ] = pulserTime;
1677 int sector = sideIndex / 1000;
1678 int plane = ( sideIndex % 1000) / 100;
1679 int counter = ( sideIndex % 100) / 10;
1680 int key = 9 * ( sector - 13 ) + 3 * ( plane - 1 ) + ( counter - 1 );
1681 if( mNPulsersCounter.count( key ) ){
1682 mNPulsersCounter[key]++;
1684 mNPulsersCounter[key] = 1;
1687 mPulserPresent[ sideIndex ] =
true;
1690 double referenceTime = 0.;
1693 if( pulserTimes.count( mReferencePulserIndex ) ) {
1694 referenceTime = pulserTimes.at( mReferencePulserIndex );
1697 LOG_INFO <<
"preliminary reference time:" << referenceTime << endm;
1704 if( mUsePulserGbtxDiff ) {
1705 for(
int gbtxId = 0; gbtxId < eTofConst::nModules * eTofConst::nSides; gbtxId++ ) {
1706 int sector = eTofConst::sectorStart + gbtxId / ( eTofConst::nPlanes * eTofConst::nSides );
1707 int zPlane = eTofConst::zPlaneStart + ( gbtxId % ( eTofConst::nPlanes * eTofConst::nSides ) ) / eTofConst::nSides;
1708 int side = eTofConst::sideStart + gbtxId % eTofConst::nSides;
1710 int partialKey = sector * 1000 + zPlane * 100 + side;
1712 vector< double > times( eTofConst::nCounters );
1714 double average = 0.;
1717 for(
int counter = 1; counter <= eTofConst::nCounters; counter++ ) {
1718 int key = partialKey + 10 * counter;
1720 if( pulserTimes.count( key ) ) {
1722 if( pulserTimes.count( partialKey + 10 ) ) {
1723 mHistograms.at(
"pulserDigiTimeToCounter1" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, pulserTimes.at( partialKey + 10 ) - pulserTimes.at( key ) );
1725 if( pulserTimes.count( partialKey + 20 ) ) {
1726 mHistograms.at(
"pulserDigiTimeToCounter2" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, pulserTimes.at( partialKey + 20 ) - pulserTimes.at( key ) );
1728 if( pulserTimes.count( partialKey + 30 ) ) {
1729 mHistograms.at(
"pulserDigiTimeToCounter3" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, pulserTimes.at( partialKey + 30 ) - pulserTimes.at( key ) );
1733 times.at( counter - 1 ) = pulserTimes.at( key ) + mPulserTimeDiffGbtx.at( key );
1735 bool isNonSmearedPulser =
false;
1736 if( referenceTime != 0 ) {
1737 double dist = times.at( counter - 1 ) - referenceTime;
1738 double redDist = mHistograms.at(
"pulserDigiTimeDiff_GbtxCorrProfMod" )->GetBinContent( gbtxId * eTofConst::nCounters + counter );
1740 double modDist = fmod( dist - redDist, eTofConst::coarseClockCycle );
1741 modDist = modDist - eTofConst::coarseClockCycle * round( modDist / eTofConst::coarseClockCycle );
1744 if( redDist == 0 || fabs( modDist ) > 0.5 ) {
1745 if( redDist != 0 ) LOG_INFO <<
"too far away --> smeared pulser: " << key <<
"(" << gbtxId <<
"-" << counter <<
")" << endm;
1750 isNonSmearedPulser =
true;
1753 mHistograms.at(
"pulserDigiTimeDiff_GbtxCorrProf" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, dist );
1754 mHistograms.at(
"pulserDigiTimeDiff_GbtxCorrProfMod" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, redDist );
1757 mHistograms.at(
"pulserDigiTimeDiff_GbtxCorr" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, dist );
1758 mHistograms.at(
"pulserDigiTimeDiff" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, pulserTimes.at( key ) - referenceTime );
1763 if( isNonSmearedPulser ) {
1764 average += times.at( counter - 1 );
1768 times.at( counter - 1 ) = 0.;
1775 if( nAv == eTofConst::nCounters ) {
1776 double diff12 = fabs( times.at( 0 ) - times.at( 1 ) );
1777 double diff13 = fabs( times.at( 0 ) - times.at( 2 ) );
1778 double diff23 = fabs( times.at( 1 ) - times.at( 2 ) );
1780 if( diff12 > 0.2 && diff13 > 0.2 && diff23 < 0.2 ) {
1781 average -= times.at( 0 );
1784 if( fabs( times.at( 0 ) - ( average / nAv ) + eTofConst::coarseClockCycle ) < 0.2 ) {
1785 mJumpingPulsers[ partialKey + 10 ] = 1;
1787 times.at( 0 ) += eTofConst::coarseClockCycle;
1788 average += times.at( 0 );
1791 else if( fabs( times.at( 0 ) - ( average / nAv ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1792 mJumpingPulsers[ partialKey + 10 ] = -1;
1794 times.at( 0 ) -= eTofConst::coarseClockCycle;
1795 average += times.at( 0 );
1800 mHistograms.at(
"1_off" )->Fill( gbtxId * eTofConst::nCounters + 1.5, times.at( 0 ) - ( average / nAv ) );
1803 else if( diff12 > 0.2 && diff13 < 0.2 && diff23 > 0.2 ) {
1804 average -= times.at( 1 );
1807 if( fabs( times.at( 1 ) - ( average / nAv ) + eTofConst::coarseClockCycle ) < 0.2 ) {
1808 mJumpingPulsers[ partialKey + 20 ] = 1;
1810 times.at( 1 ) += eTofConst::coarseClockCycle;
1811 average += times.at( 1 );
1814 else if( fabs( times.at( 1 ) - ( average / nAv ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1815 mJumpingPulsers[ partialKey + 20 ] = -1;
1817 times.at( 1 ) -= eTofConst::coarseClockCycle;
1818 average += times.at( 1 );
1823 mHistograms.at(
"2_off" )->Fill( gbtxId * eTofConst::nCounters + 1.5, times.at( 1 ) - ( average / nAv ) );
1826 else if( diff12 < 0.2 && diff13 > 0.2 && diff23 > 0.2 ) {
1827 average -= times.at( 2 );
1830 if( fabs( times.at( 2 ) - ( average / nAv ) + eTofConst::coarseClockCycle ) < 0.2 ) {
1831 mJumpingPulsers[ partialKey + 30 ] = 1;
1833 times.at( 2 ) += eTofConst::coarseClockCycle;
1834 average += times.at( 2 );
1837 else if( fabs( times.at( 2 ) - ( average / nAv ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1838 mJumpingPulsers[ partialKey + 30 ] = -1;
1840 times.at( 2 ) -= eTofConst::coarseClockCycle;
1841 average += times.at( 2 );
1846 mHistograms.at(
"3_off" )->Fill( gbtxId * eTofConst::nCounters + 1.5, times.at( 2 ) - ( average / nAv ) );
1851 if( nAv == eTofConst::nCounters - 1 ) {
1853 if( times.at( 0 ) > 0 && times.at( 1 ) > 0 && fabs( fabs( times.at( 0 ) - times.at( 1 ) ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1854 if( mJumpingPulsers.count( partialKey + 10 ) ) {
1856 times.at( 0 ) += mJumpingPulsers.at( partialKey + 10 ) * eTofConst::coarseClockCycle;
1857 average += mJumpingPulsers.at( partialKey + 10 ) * eTofConst::coarseClockCycle;
1859 else if( mJumpingPulsers.count( partialKey + 20 ) ) {
1861 times.at( 1 ) += mJumpingPulsers.at( partialKey + 20 ) * eTofConst::coarseClockCycle;
1862 average += mJumpingPulsers.at( partialKey + 20 ) * eTofConst::coarseClockCycle;
1866 if( times.at( 0 ) < times.at( 1 ) ) {
1867 times.at( 0 ) += eTofConst::coarseClockCycle;
1868 average += eTofConst::coarseClockCycle;
1871 times.at( 1 ) += eTofConst::coarseClockCycle;
1872 average += eTofConst::coarseClockCycle;
1876 else if( times.at( 0 ) && times.at( 2 ) > 0 && fabs( fabs( times.at( 0 ) - times.at( 2 ) ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1877 if( mJumpingPulsers.count( partialKey + 10 ) ) {
1879 times.at( 0 ) += mJumpingPulsers.at( partialKey + 10 ) * eTofConst::coarseClockCycle;
1880 average += mJumpingPulsers.at( partialKey + 10 ) * eTofConst::coarseClockCycle;
1882 else if( mJumpingPulsers.count( partialKey + 30 ) ) {
1884 times.at( 2 ) += mJumpingPulsers.at( partialKey + 30 ) * eTofConst::coarseClockCycle;
1885 average += mJumpingPulsers.at( partialKey + 30 ) * eTofConst::coarseClockCycle;
1889 if( times.at( 0 ) < times.at( 2 ) ) {
1890 times.at( 0 ) += eTofConst::coarseClockCycle;
1891 average += eTofConst::coarseClockCycle;
1894 times.at( 2 ) += eTofConst::coarseClockCycle;
1895 average += eTofConst::coarseClockCycle;
1899 else if( times.at( 1 ) > 0 && times.at( 2 ) > 0 && fabs( fabs( times.at( 1 ) - times.at( 2 ) ) - eTofConst::coarseClockCycle ) < 0.2 ) {
1900 if( mJumpingPulsers.count( partialKey + 20 ) ) {
1902 times.at( 1 ) += mJumpingPulsers.at( partialKey + 20 ) * eTofConst::coarseClockCycle;
1903 average += mJumpingPulsers.at( partialKey + 20 ) * eTofConst::coarseClockCycle;
1905 else if( mJumpingPulsers.count( partialKey + 30 ) ) {
1907 times.at( 2 ) += mJumpingPulsers.at( partialKey + 30 ) * eTofConst::coarseClockCycle;
1908 average += mJumpingPulsers.at( partialKey + 30 ) * eTofConst::coarseClockCycle;
1912 if( times.at( 1 ) < times.at( 2 ) ) {
1913 times.at( 1 ) += eTofConst::coarseClockCycle;
1914 average += eTofConst::coarseClockCycle;
1917 times.at( 2 ) += eTofConst::coarseClockCycle;
1918 average += eTofConst::coarseClockCycle;
1929 if( mDoQA && referenceTime != 0 ) {
1930 mHistograms.at(
"pulserDigiTimeDiff_perGbtx" )->Fill( gbtxId * eTofConst::nCounters + 1.5, average - referenceTime );
1933 for(
int counter = eTofConst::counterStart; counter <= eTofConst::counterStop; counter++ ) {
1934 double diffToAv = 0.;
1936 if( times.at( counter - eTofConst::counterStart ) != 0. ) {
1937 diffToAv = times.at( counter - eTofConst::counterStart ) - average;
1939 if( fabs( diffToAv ) > 0.2 ) diffToAv = 0.;
1942 mHistograms.at(
"pulserDigiTimeDiff_toAverage" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, diffToAv );
1946 if( average != 0. ) {
1948 pulserTimes[ partialKey + 10 * counter ] = average + diffToAv - mPulserTimeDiffGbtx.at( partialKey + 10 * counter );
1952 if( pulserTimes.count( partialKey + 10 * counter ) ) {
1953 pulserTimes.erase( partialKey + 10 * counter );
1963 if( pulserTimes.count( mReferencePulserIndex ) ) {
1964 if( mPulserTimeDiff.count( mReferencePulserIndex ) ){
1965 referenceTime = pulserTimes.at( mReferencePulserIndex ) - mPulserTimeDiff[ mReferencePulserIndex ];
1968 referenceTime = pulserTimes.at( mReferencePulserIndex );
1975 if( referenceTime != 0 ) {
1976 int iLockThreshold = 0;
1977 mHistograms[
"pulserDigiTimeDiff_RefCorr" ]->Reset(
"ICESM");
1979 for(
auto& kv : pulserTimes ) {
1981 if( mPulserTimeDiff.count( kv.first ) ) {
1985 double modDist = mPulserTimeDiff.at( kv.first ) - ( kv.second - referenceTime );
1989 if( fabs( modDist ) > 0.2 ) {
1990 mUnlockPulserState[ kv.first ]++;
1996 if( mUnlockPulserState.at( kv.first ) < iLockThreshold ) {
2005 mHistograms[
"pulserDigiTimeDiff_RefCorr" ]->Fill( kv.second - referenceTime );
2009 if( mUnlockPulserState.count( kv.first ) ) {
2010 LOG_INFO <<
" --> new event doesn't have offset for pulser " << kv.first <<
" --> remove the entry" << endm;
2011 mUnlockPulserState.erase( kv.first );
2019 int sector = kv.first / 1000;
2020 int zPlane = ( kv.first % 1000 ) / 100;
2021 int counter = ( kv.first % 100 ) / 10;
2022 int side = kv.first % 10;
2024 int gbtxId = ( sector - eTofConst::sectorStart ) * ( eTofConst::nPlanes * eTofConst::nSides )
2025 + ( zPlane - eTofConst::zPlaneStart ) * eTofConst::nSides
2026 + ( side - eTofConst::sideStart );
2028 if( mPulserTimeDiff.count( kv.first ) ) {
2029 mHistograms.at(
"pulserDigiTimeDiff_fullCorr" )->Fill( gbtxId * eTofConst::nCounters + counter - 0.5, mPulserTimeDiff.at( kv.first ) );
2032 mHistograms[
"pulserDigiPresence" ]->Fill(gbtxId * eTofConst::nCounters + counter - 0.5);
2038 mHistograms[
"pulserDigiPresence" ]->Fill( -1 );
2039 mHistograms[
"pulserDigiTimeDiff_CorrAgreement" ]->Fill( mHistograms[
"pulserDigiTimeDiff_RefCorr" ]->GetMaximum() );
2040 mHistograms[
"pulserDigiTimeDiff_CorrCommonJump" ]->Fill( mHistograms[
"pulserDigiTimeDiff_RefCorr" ]->GetEntries() );
2042 if( ( mHistograms[
"pulserDigiTimeDiff_RefCorr" ]->GetEntries() > 150 ) && ( mHistograms[
"pulserDigiTimeDiff_RefCorr" ]->GetMaximum() > 15 ) ){
2045 int iMaximumBin = mHistograms[
"pulserDigiTimeDiff_RefCorr" ]->GetMaximumBin();
2046 double dRefCorrAverage = 0;
2047 int iNRefCorrAgree = 0;
2049 for(
auto& kv : pulserTimes ) {
2050 if ( mHistograms[
"pulserDigiTimeDiff_RefCorr" ]->FindBin(kv.second - referenceTime) == iMaximumBin ){
2051 dRefCorrAverage += kv.second - referenceTime;
2055 dRefCorrAverage /= iNRefCorrAgree;
2056 referenceTime -= dRefCorrAverage;
2058 mPulserTimeDiff[ mReferencePulserIndex ] -= dRefCorrAverage;
2062 for(
auto& kv : pulserTimes ) {
2063 if ( ! mPulserTimeDiff[ kv.first ] ) {
2064 mPulserTimeDiff[ kv.first ] = kv.second - referenceTime;
2067 if( mUnlockPulserState.count( kv.first ) ){
2068 if( mUnlockPulserState.at( kv.first ) > iLockThreshold ){
2069 mPulserTimeDiff[ kv.first ] = kv.second - referenceTime;
2086 if( !mStatus.count( key) || mStatus.at( key ) != 1 ) {
2088 LOG_DEBUG <<
"status of channel with key " << key <<
" was not ok ---> skip calibrating this digi" << endm;
2094 if( fabs( aDigi->
calibTot() + 999. ) < 1.e-5 ) {
2096 LOG_DEBUG <<
"digi flaged as pulser --> skip" << endm;
2101 float timeToTrigger = aDigi->
rawTime() - mTriggerTime;
2104 if( timeToTrigger > mTimingWindow.at( aDigi->
rocId() ).first &&
2105 timeToTrigger < mTimingWindow.at( aDigi->
rocId() ).second )
2108 if( mStrictPulserHandling ){
2110 if( !mPulserPresent.count( PulserKey ) ) {
2112 LOG_DEBUG <<
"no pulser in the same event for this counter --> digi skipped due to strict pulser handling" << endm;
2118 double calibTot = aDigi->
rawTot() * mGet4TotBinWidthNs * calibTotFactor( aDigi );
2120 aDigi->setCalibTot( calibTot );
2122 double calibTime = aDigi->
rawTime() - mResetTime
2124 - calibTimeOffset( aDigi )
2125 - slewingTimeOffset( aDigi )
2126 - applyPulserOffset( aDigi );
2128 aDigi->setCalibTime( calibTime );
2132 LOG_DEBUG <<
"raw Time, ToT: " << aDigi->
rawTime() <<
", " << aDigi->
rawTot() << endm;
2133 LOG_DEBUG <<
"calibrated Time, ToT: " << aDigi->
calibTime() <<
", " << aDigi->
calibTot() << endm;
2139 LOG_DEBUG <<
"digi is outside the timing window (time to trigger = " << timeToTrigger <<
") --> skip" << endm;
2150 if( !aDigi )
return;
2153 aDigi->setCalibTime( 0. );
2154 aDigi->setCalibTot( -1. );
2156 aDigi->setAssociatedHitId( -1 );
2170 StETofCalibMaker::flagPulserDigis(
StMuETofDigi* aDigi,
unsigned int index, std::map<
unsigned int, std::vector< unsigned int > >& pulserDigiMap )
2172 flagPulserDigis( (
StETofDigi* ) aDigi, index, pulserDigiMap );
2190 StETofCalibMaker::calibTotFactor(
StETofDigi* aDigi )
2193 unsigned int bin = aDigi->
strip() + 32 * ( aDigi->
side() - 1 );
2195 if( mDigiTotCorr.count( key ) ) {
2196 float binContent = mDigiTotCorr.at( key )->GetBinContent( bin );
2198 if( fabs( binContent ) > 1e-5 ) {
2200 LOG_DEBUG <<
"calibTotFactor: histogram with key " << key <<
" at bin " << bin <<
" -> return bin content: " << binContent << endm;
2202 return (1.0/binContent);
2206 LOG_WARN <<
"calibTotFactor: histogram with key " << key <<
" at bin " << bin <<
" has content of 0 -> return 1" << endm;
2213 LOG_WARN <<
"calibTotFactor: required histogram with key " << key <<
" doesn't exist -> return 1" << endm;
2225 StETofCalibMaker::calibTimeOffset(
StETofDigi* aDigi )
2228 unsigned int bin = aDigi->
strip() + 32 * ( aDigi->
side() - 1 );
2230 if( mDigiTimeCorr.count( key ) ) {
2231 float binContent = mDigiTimeCorr.at( key )->GetBinContent( bin );
2233 LOG_DEBUG <<
"calibTimeOffset: histogram with key " << key <<
" at bin " << bin <<
" -> return bin content: " << binContent << endm;
2239 LOG_WARN <<
"calibTimeOffset: required histogram with key " << key <<
" doesn't exist -> return 0" << endm;
2251 StETofCalibMaker::slewingTimeOffset(
StETofDigi* aDigi )
2255 if( mDigiSlewCorr.count( key ) ) {
2257 unsigned int totBin = mDigiSlewCorr.at( key )->FindBin( aDigi->
rawTot() );
2259 if( mDigiSlewCorr.at( key )->GetBinEntries( totBin ) <= mMinDigisPerSlewBin && totBin < etofSlewing::nTotBins ) {
2261 LOG_DEBUG <<
"slewingTimeOffset: insufficient statistics for slewing calibration in channel " << key <<
" at tot bin " << totBin <<
" --> return 0" << endm;
2266 float val = mDigiSlewCorr.at( key )->Interpolate( aDigi->
rawTot() );
2268 LOG_DEBUG <<
"slewingTimeOffset: histogram with key " << key <<
" with calib TOT of " << aDigi->
calibTot() <<
" --> interpolated correction: " << val << endm;
2274 LOG_DEBUG <<
"slewingTimeOffset: required histogram with key " << key <<
" doesn't exist -> return 0" << endm;
2286 StETofCalibMaker::applyPulserOffset(
StETofDigi* aDigi )
2290 if( !mPulserTimeDiff.count( key )) {
2294 return mPulserTimeDiff.at( key );
2306 double triggerTime = header->trgGdpbFullTime();
2309 std::map< uint64_t, short > countsGdpbTs;
2310 for(
const auto& kv : header->rocGdpbTs() ) {
2312 LOG_DEBUG <<
"triggerTime (" << std::hex <<
"Ox" << kv.first << std::dec <<
") " << kv.second * eTofConst::coarseClockCycle * 1.e-9 << endm;
2314 ++countsGdpbTs[ kv.second ];
2321 uint64_t mostProbableTriggerTs = 0;
2323 for(
auto it = countsGdpbTs.begin(); it != countsGdpbTs.end(); it++ ) {
2324 auto next = std::next( it, 1 );
2325 auto prev = std::prev( it, 1 );
2327 short countTs = it->second;
2329 if( next != countsGdpbTs.end() && ( next->first - it->first ) == 1 ) {
2330 countTs += next->second;
2332 if( accCount > 0 && ( it->first - prev->first ) == 1 ) {
2333 countTs += prev->second;
2336 if( countTs >= accCount ) {
2339 if( it->second > maxCount ) {
2340 maxCount = it->second;
2341 mostProbableTriggerTs = it->first;
2346 if( mostProbableTriggerTs > 0 ) {
2347 triggerTime = mostProbableTriggerTs * eTofConst::coarseClockCycle;
2351 LOG_DEBUG <<
"trigger TS: " << mostProbableTriggerTs <<
" --> trigger time (ns): " << triggerTime << endm;
2366 std::map< uint64_t, short > countsStarTsRaw;
2367 for(
const auto& kv : header->rocStarTs() ) {
2369 LOG_DEBUG <<
"resetTime (" << std::hex <<
"Ox" << kv.first << std::dec <<
") " << kv.second * eTofConst::coarseClockCycle * 1.e-9 << endm;
2373 if( mRunYear == 2018 && kv.first != 0x18e6 )
continue;
2375 if( kv.second != 0 ) {
2376 ++countsStarTsRaw[ kv.second ];
2382 std::map< uint64_t, short > countsStarTs;
2383 for(
auto it = countsStarTsRaw.begin(); it != countsStarTsRaw.end(); it++ ) {
2384 auto next = std::next( it, 1 );
2386 short countTs = it->second;
2388 if( next != countsStarTsRaw.end() && ( next->first - it->first ) == 1 ) {
2389 countTs += next->second;
2392 countsStarTs[ it->first ] = countTs;
2399 if( countsStarTs.size() == 0 ) {
2400 LOG_INFO <<
"all AFCKs report a reset time value 0" << endm;
2403 for(
const auto& kv : countsStarTs ) {
2404 LOG_DEBUG <<
"resetTime cand:" << kv.first <<
" (" << kv.second <<
" times)" << endm;
2405 mHistograms.at(
"resetTimeCand_times" )->Fill( kv.second );
2408 for(
const auto& kv : header->rocStarTs() ) {
2409 unsigned int sector;
2410 mHwMap->mapToSector( kv.first, sector );
2412 LOG_DEBUG <<
"resetTime(" << sector <<
"): " << kv.second << endm;
2414 std::string histName =
"resetTimeDifferenceToSector" + to_string( sector );
2415 for(
const auto& jv : header->rocStarTs() ) {
2417 mHwMap->mapToSector( jv.first, sec );
2419 mHistograms.at( histName )->Fill( sec, ( int64_t ) jv.second - ( int64_t ) kv.second );
2425 while( countsStarTs.size() > 0 ) {
2426 auto it = std::max_element( countsStarTs.begin(), countsStarTs.end(),
2427 [](
const pair< uint64_t, short >& p1,
const pair< uint64_t, short >& p2 ) {
2428 return p1.second < p2.second; } );
2430 double resetTime = it->first * eTofConst::coarseClockCycle;
2434 if( abs( mResetTs - ( int64_t ) it->first ) < 2 ) {
2435 resetTime = mResetTs * eTofConst::coarseClockCycle;
2438 LOG_INFO <<
"change in reset TS: " << mResetTs <<
" --> " << it->first << endm;
2439 mResetTs = it->first;
2445 if( mTriggerTime - resetTime < 7.2e12 ) {
2447 LOG_DEBUG <<
"reset time (ns): " << resetTime <<
" --> difference to trigger time in seconds: " << ( mTriggerTime - resetTime ) * 1.e-9 << endm;
2449 LOG_DEBUG <<
"--> picked reset TS:" << mResetTs << endm;
2452 mHistograms.at(
"resetTimeCand_picked" )->Fill( it->second );
2454 auto rawIt = std::max_element( countsStarTsRaw.begin(), countsStarTsRaw.end(),
2455 [](
const pair< uint64_t, short >& p1,
const pair< uint64_t, short >& p2 ) {
2456 return p1.second < p2.second; } );
2458 mHistograms.at(
"resetTimeCand_compare" )->Fill( ( int64_t ) mResetTs - ( int64_t ) rawIt->first );
2459 mHistograms.at(
"resetTimeCand_value" )->Fill( mResetTs % (
int ) eTofConst::bTofClockCycle );
2465 countsStarTs.erase( it );
2475 StETofCalibMaker::channelToKey(
const unsigned int channelId ) {
2476 unsigned int sector = ( channelId / eTofConst::nChannelsPerSector ) + eTofConst::sectorStart;
2477 unsigned int zPlane = ( ( channelId % eTofConst::nChannelsPerSector ) / eTofConst::nChannelsPerModule ) + eTofConst::zPlaneStart;
2478 unsigned int counter = ( ( channelId % eTofConst::nChannelsPerModule ) / eTofConst::nChannelsPerCounter ) + eTofConst::counterStart;
2479 unsigned int strip = ( ( channelId % eTofConst::nChannelsPerCounter ) / eTofConst::nSides ) + eTofConst::stripStart;
2480 unsigned int side = ( channelId % eTofConst::nSides ) + eTofConst::sideStart;
2482 return sector * 100000 + zPlane * 10000 + counter * 1000 + strip * 10 + side;
2488 StETofCalibMaker::detectorToKey(
const unsigned int detectorId ) {
2489 unsigned int sector = ( detectorId / eTofConst::nCountersPerSector ) + eTofConst::sectorStart;
2490 unsigned int zPlane = ( ( detectorId % eTofConst::nCountersPerSector ) / eTofConst::nCounters ) + eTofConst::zPlaneStart;
2491 unsigned int counter = ( detectorId % eTofConst::nCounters ) + eTofConst::counterStart;
2493 return sector * 100 + zPlane * 10 + counter;
2499 StETofCalibMaker::sideToKey(
const unsigned int sideId ) {
2500 unsigned int sector = ( sideId / ( eTofConst::nCountersPerSector * eTofConst::nSides ) ) + eTofConst::sectorStart;
2501 unsigned int zPlane = ( ( sideId % ( eTofConst::nCountersPerSector * eTofConst::nSides ) ) / ( eTofConst::nCounters * eTofConst::nSides ) ) + eTofConst::zPlaneStart;
2502 unsigned int counter = ( ( sideId % ( eTofConst::nCounters * eTofConst::nSides ) ) / eTofConst::nSides ) + eTofConst::counterStart;
2503 unsigned int side = ( sideId % eTofConst::nSides ) + eTofConst::sideStart;
2505 return sector * 1000 + zPlane * 100 + counter * 10 + side;
2514 StETofCalibMaker::setHistFileName()
2516 string extension =
".etofCalib.root";
2518 if( GetChainOpt()->GetFileOut() !=
nullptr ) {
2519 TString outFile = GetChainOpt()->GetFileOut();
2521 mHistFileName = ( string ) outFile;
2524 size_t lastindex = mHistFileName.find_last_of(
"." );
2525 mHistFileName = mHistFileName.substr( 0, lastindex );
2528 lastindex = mHistFileName.find_last_of(
"." );
2529 mHistFileName = mHistFileName.substr( 0, lastindex );
2532 lastindex = mHistFileName.find_last_of(
"/" );
2533 mHistFileName = mHistFileName.substr( lastindex + 1, mHistFileName.length() );
2535 mHistFileName = mHistFileName + extension;
2537 LOG_ERROR <<
"Cannot set the output filename for histograms" << endm;
2543 StETofCalibMaker::bookHistograms()
2545 LOG_INFO <<
"bookHistograms() ... " << endm;
2547 mHistograms[
"pulserDigiTimeDiff_GbtxCorrProf" ] =
new TProfile(
"pulserDigiTimeDiff_GbtxCorrProf",
"time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216,
"s" );
2548 mHistograms[
"pulserDigiTimeDiff_GbtxCorrProfMod" ] =
new TProfile(
"pulserDigiTimeDiff_GbtxCorrProfMod",
"time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216,
"s" );
2549 mHistograms[
"pulserDigiTimeDiff_fullCorr" ] =
new TH2F(
"pulserDigiTimeDiff_fullCorr",
"time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2550 mHistograms[
"pulserDigiTimeDiff_RefCorr" ] =
new TH1F(
"pulserDigiTimeDiff_RefCorr",
"time difference of pulsers to reference; #Delta T (ns)", 45, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ));
2553 mHistograms[
"pulserDigiTimeDiff" ] =
new TH2F(
"pulserDigiTimeDiff",
"time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2554 mHistograms[
"pulserDigiTimeToCounter1" ] =
new TH2F(
"pulserDigiTimeToCounter1",
"time difference of pulsers to counter 1;pulser channel;#Delta T (ns)", 216, 0, 216, 2*360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2555 mHistograms[
"pulserDigiTimeToCounter2" ] =
new TH2F(
"pulserDigiTimeToCounter2",
"time difference of pulsers to counter 2;pulser channel;#Delta T (ns)", 216, 0, 216, 2*360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2556 mHistograms[
"pulserDigiTimeToCounter3" ] =
new TH2F(
"pulserDigiTimeToCounter3",
"time difference of pulsers to counter 3;pulser channel;#Delta T (ns)", 216, 0, 216, 2*360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2558 mHistograms[
"pulserDigiTimeDiff_GbtxCorr" ] =
new TH2F(
"pulserDigiTimeDiff_GbtxCorr",
"time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2559 mHistograms[
"pulserDigiTimeDiff_perGbtx" ] =
new TH2F(
"pulserDigiTimeDiff_perGbtx",
"average time difference of pulsers in one Gbtx to reference;pulser channel;#Delta T (ns)", 216/3, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2560 mHistograms[
"pulserDigiTimeDiff_toAverage" ] =
new TH2F(
"pulserDigiTimeDiff_toAverage",
"time difference of pulsers to reference;pulser channel;#Delta T (ns)", 216, 0, 216, 4*360, -359.5 * ( 6.25 / 112 ), 360.5 * ( 6.25 / 112 ) );
2562 mHistograms[
"1_off" ] =
new TH2F(
"1_off",
"average time difference of pulsers in one Gbtx to reference;pulser channel;#Delta T (ns)", 216/3, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2563 mHistograms[
"2_off" ] =
new TH2F(
"2_off",
"average time difference of pulsers in one Gbtx to reference;pulser channel;#Delta T (ns)", 216/3, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2564 mHistograms[
"3_off" ] =
new TH2F(
"3_off",
"average time difference of pulsers in one Gbtx to reference;pulser channel;#Delta T (ns)", 216/3, 0, 216, 360, -179.5 * ( 6.25 / 112 ), 180.5 * ( 6.25 / 112 ) );
2566 mHistograms[
"pulserDigiTimeDiff_CorrAgreement" ] =
new TH1F(
"pulserDigiTimeDiff_CorrAgreement",
"Number of pulsers agreeing on a common shift between events; #pulsers", 218, -0.5, 217.5);
2567 mHistograms[
"pulserDigiTimeDiff_CorrCommonJump" ] =
new TH1F(
"pulserDigiTimeDiff_CorrCommonJump",
"Number of pulsers jumping at the same time between events; #pulsers", 218, -0.5, 217.5);
2568 mHistograms[
"pulserDigiPresence" ] =
new TH1F(
"pulserDigiPresence",
"pulser presence (number of events at ( -1 );pulser channel", 218, -1.5, 216.5);
2570 for(
int i=0; i<12; i++ ) {
2571 std::string histName =
"resetTimeDifferenceToSector" + to_string( i + 13 );
2572 mHistograms[ histName ] =
new TH2F( Form(
"resetTimeDifferenceToSector%d", i + 13 ), Form(
"reset time difference to sector %d;sector;#DeltaT (clock ticks)", i + 13 ), 12, 12.5, 24.4, 5, -2.5, 2.5 );
2574 mHistograms[
"ETOF_QA_daqMissmatches_get4" ] =
new TProfile(
"ETOF_QA_daqMissmatches_get4",
"missmatch percentage for each get4; get4 []; missmatch percentage", 1728, 0.5, 1728.5 );
2575 mHistograms[
"ETOF_QA_daqMissmatches_counter" ] =
new TProfile(
"ETOF_QA_daqMissmatches_counter",
"missmatch percentage for each counter; counter []; missmatch percentage", 108, 0.5, 108.5 );
2576 mHistograms[
"resetTimeCand_times" ] =
new TH1F(
"resetTimeCand_times",
"sectors agreeing on reset time candidates;# sectors with common candidate;# events", 12, 0.5, 12.5 );
2577 mHistograms[
"resetTimeCand_picked" ] =
new TH1F(
"resetTimeCand_picked",
"sectors agreeing on picked reset time;# sectors with picked reset time;# events", 12, 0.5, 12.5 );
2578 mHistograms[
"resetTimeCand_compare" ] =
new TH1F(
"resetTimeCand_compare",
"difference between old and new way;#DeltaT (clock ticks);# events", 5, -2.5, 2.5 );
2579 mHistograms[
"resetTimeCand_value" ] =
new TH1F(
"resetTimeCand_value",
"picked reset time value;clock ticks;# events", (
int ) eTofConst::bTofClockCycle, 0.5, 0.5 + eTofConst::bTofClockCycle );
2582 for (
auto& kv : mHistograms ) {
2583 kv.second->SetDirectory( 0 );
2589 StETofCalibMaker::writeHistograms()
2591 if( mHistFileName !=
"" ) {
2592 LOG_INFO <<
"writing histograms to: " << mHistFileName.c_str() << endm;
2594 TFile histFile( mHistFileName.c_str(),
"RECREATE",
"etofCalib" );
2597 for(
int i=0; i<12; i++ ) {
2598 std::string histName =
"resetTimeDifferenceToSector" + to_string( i + 13 );
2599 mHistograms.at( histName )->Scale( 12. / mHistograms.at( histName )->GetEntries() );
2602 for (
const auto& kv : mHistograms ) {
2603 if( kv.second->GetEntries() > 0 ) kv.second->Write();
2609 LOG_INFO <<
"histogram file name is empty string --> cannot write histograms" << endm;
unsigned int sector() const
Sector.
double calibTime() const
calibrated time
unsigned int zPlane() const
ZPlane.
static StMuETofHeader * etofHeader()
returns pointer to the StMuETofHeader
unsigned int side() const
Side.
unsigned int strip() const
Strip.
double rawTot() const
Getter for uncalibrated Tot.
unsigned int sector() const
Sector.
unsigned int zPlane() const
ZPlane.
unsigned int strip() const
Strip.
double rawTime() const
Raw Time.
double rawTot() const
Getter for uncalibrated Tot.
static TClonesArray * etofArray(int type)
returns pointer to the n-th TClonesArray from the etof arrays // FS
static StMuETofDigi * etofDigi(int i)
returns pointer to the i-th StMuEtofDigi
unsigned int elChan() const
electronic Channel.
unsigned int get4Id() const
get4Id.
unsigned int counter() const
Counter.
void setGeoAddress(const unsigned int iSector, const unsigned int iZPlane, const unsigned int iCounter, const unsigned int iChannel, const unsigned int iSide)
double calibTot() const
Getter for calibrated Tot.
unsigned int rocId() const
RocId.
unsigned int side() const
Side.
StETofCalibMaker(const char *name="etofCalib")
default constructor
unsigned int counter() const
Counter.
void setGeoAddress(const unsigned int iSector, const unsigned int iZPlane, const unsigned int iCounter, const unsigned int iChannel, const unsigned int iSide)
virtual TDataSet * Find(const char *path) const
double rawTime() const
Raw Time.