diff --git a/TileCalorimeter/TileExample/TileRecEx/share/jobOptions_TileCalibRec.py b/TileCalorimeter/TileExample/TileRecEx/share/jobOptions_TileCalibRec.py index ade24cf238cb37caf405969b1056217df229ea96..08125adf4081c1597be9a779f1119b72c33ae5be 100644 --- a/TileCalorimeter/TileExample/TileRecEx/share/jobOptions_TileCalibRec.py +++ b/TileCalorimeter/TileExample/TileRecEx/share/jobOptions_TileCalibRec.py @@ -468,6 +468,9 @@ if not 'doTileMF' in dir(): if not 'doTileOF1' in dir(): doTileOF1 = False +if not 'doTileWiener' in dir(): + doTileWiener = False + if not 'doTileFit' in dir(): doTileFit = not TileCompareMode and ReadDigits @@ -513,7 +516,7 @@ if not 'OfcFromCOOL' in dir(): else: OfcFromCOOL = False -if useRODReco or doTileOpt2 or doTileMF or doTileOF1 or doTileOptATLAS or doTileFitCool or TileCompareMode or not 'TileUseCOOL' in dir(): +if useRODReco or doTileOpt2 or doTileMF or doTileOF1 or doTileOptATLAS or doTileWiener or doTileFitCool or TileCompareMode or not 'TileUseCOOL' in dir(): TileUseCOOL = True TileUseCOOLOFC = not ReadPool or OfcFromCOOL @@ -835,12 +838,19 @@ else: from xAODEventInfoCnv.xAODEventInfoCreator import xAODMaker__EventInfoCnvAlg topSequence+=xAODMaker__EventInfoCnvAlg() +#============================================================ +#=== configure BunchCrossingTool +#============================================================ +from TrigBunchCrossingTool.BunchCrossingTool import BunchCrossingTool +ToolSvc += BunchCrossingTool("LHC" if globalflags.DataSource() == "data" else "MC") + #============================================================= #=== read ByteStream and reconstruct data #============================================================= tileRawChannelBuilderFitFilter = None tileRawChannelBuilderOpt2Filter = None tileRawChannelBuilderOptATLAS = None +tileRawChannelBuilderWienerFilter = None tileDigitsContainer = '' if not ReadPool: include( "ByteStreamCnvSvcBase/BSAddProvSvc_RDO_jobOptions.py" ) @@ -864,6 +874,7 @@ else: tileRawChannelBuilderFitFilter = theTileRawChannelGetter.TileRawChannelBuilderFitFilter() tileRawChannelBuilderOpt2Filter = theTileRawChannelGetter.TileRawChannelBuilderOpt2Filter() tileRawChannelBuilderOptATLAS = theTileRawChannelGetter.TileRawChannelBuilderOptATLAS() + tileRawChannelBuilderWienerFilter = theTileRawChannelGetter.TileRawChannelBuilderWienerFilter() if doRecoESD: topSequence.TileRChMaker.TileDigitsContainer="TileDigitsFlt" tileDigitsContainer = 'TileDigitsFlt' @@ -940,6 +951,20 @@ if doTileOF1: print ToolSvc.TileRawChannelBuilderOF1 +if doTileWiener and tileRawChannelBuilderWienerFilter: + if PhaseFromCOOL: + tileRawChannelBuilderWienerFilter.correctTime = False # do not need to correct time with best phase + + tileRawChannelBuilderWienerFilter.BestPhase = PhaseFromCOOL # Phase from COOL or assume phase=0 + + if TileMonoRun or TileRampRun: + if TileCompareMode or TileEmulateDSP: + tileRawChannelBuilderWienerFilter.EmulateDSP = True # use dsp emulation + tileRawChannelBuilderWienerFilter.UseDSPCorrection = not TileBiGainRun + tileRawChannelBuilderWienerFilter.OutputLevel = 1 + + print tileRawChannelBuilderWienerFilter + if (doEventDisplay or doCreatePool): # create TileHit from TileRawChannel and store it in TileHitVec from TileRecAlgs.TileHitFromRawChGetter import * @@ -1290,10 +1315,13 @@ if doTileMon: if doTileOpt2: theTileRawChannelMon.TileRawChannelContainer = "TileRawChannelOpt2" - + if doTileMF: theTileRawChannelMon.TileRawChannelContainer = "TileRawChannelMF" + if doTileWiener: + theTileRawChannelMon.TileRawChannelContainer = "TileRawChannelWiener" + if useRODReco: theTileRawChannelMon.TileRawChannelContainerDSP = "TileRawChannelCnt" diff --git a/TileCalorimeter/TileIdentifier/TileIdentifier/TileFragHash.h b/TileCalorimeter/TileIdentifier/TileIdentifier/TileFragHash.h index ea10c7de9ab169e3c82b4e6a57cc2cbef9606e82..a8a9d4940e203af6408102bf2e8ff73d32679ce7 100755 --- a/TileCalorimeter/TileIdentifier/TileIdentifier/TileFragHash.h +++ b/TileCalorimeter/TileIdentifier/TileIdentifier/TileFragHash.h @@ -32,7 +32,8 @@ class TileFragHash { /** initialize */ enum TYPE {Beam=255, Default=0, Digitizer=0, OptFilterDsp=1, OptFilterOffline=2, OptFilterDspCompressed=3, - ManyAmps=4, MF=5, FitFilter=6, FitFilterCool=7, FlatFilter=8 }; + ManyAmps=4, MF=5, FitFilter=6, FitFilterCool=7, FlatFilter=8, + WienerFilterOffline=9 }; void initialize(const TileHWID * tileHWID, TYPE type=Default ); diff --git a/TileCalorimeter/TileRec/TileRec/TileAANtuple.h b/TileCalorimeter/TileRec/TileRec/TileAANtuple.h index d632283c02d35aa9e8329d656ffc7b8823513d78..d9c4aaa4c4a480ffd6298dd6a103a1fedab63946 100755 --- a/TileCalorimeter/TileRec/TileRec/TileAANtuple.h +++ b/TileCalorimeter/TileRec/TileRec/TileAANtuple.h @@ -22,6 +22,7 @@ /// TileRawChannelContainerFit key of fit filtered RawChannels in TDS /// TileRawChannelContainerFitCool key of fit filtered RawChannels in TDS /// TileRawChannelOF1 key of OF1 filtered RawChannels in TDS +/// TileRawChannelWiener key of Wiener filtered RawChannels in TDS /// TileLaserObject string key of Laser Object in TDS /// CalibrateEnergy bool If calibration should be applied to energy /// CalibMode bool If data is in calibration mode (bi-gain input) @@ -300,6 +301,11 @@ class TileAANtuple : public AthAlgorithm { float m_chi2MF[N_ROS2][N_MODULES][N_CHANS] = {{{0}}}; float m_pedMF[N_ROS2][N_MODULES][N_CHANS] = {{{0}}}; + float m_eWiener[N_ROS2][N_MODULES][N_CHANS] = {{{0}}}; + float m_tWiener[N_ROS2][N_MODULES][N_CHANS] = {{{0}}}; + float m_pedWiener[N_ROS2][N_MODULES][N_CHANS] = {{{0}}}; + float m_chi2Wiener[N_ROS2][N_MODULES][N_CHANS] = {{{0}}}; + short m_ROD_GlobalCRC[N_ROS][N_MODULES] = {{0}}; short m_ROD_BCID[N_ROS][N_MODULES] = {{0}}; short m_ROD_DMUBCIDErr[N_ROS][N_MODULES][N_DMUS] = {{{0}}}; @@ -349,6 +355,7 @@ class TileAANtuple : public AthAlgorithm { SG::ReadHandleKey<TileRawChannelContainer> m_dspRawChannelContainerKey; SG::ReadHandleKey<TileRawChannelContainer> m_mfRawChannelContainerKey; SG::ReadHandleKey<TileRawChannelContainer> m_of1RawChannelContainerKey; + SG::ReadHandleKey<TileRawChannelContainer> m_wienerRawChannelContainerKey; SG::ReadHandleKey<TileLaserObject> m_laserObjectKey; SG::ReadHandleKey<TileRawChannelContainer> m_tileMuRcvRawChannelContainerKey; // TMDB SG::ReadHandleKey<TileDigitsContainer> m_tileMuRcvDigitsContainerKey; // TMDB diff --git a/TileCalorimeter/TileRec/share/TileDefaults_jobOptions.py b/TileCalorimeter/TileRec/share/TileDefaults_jobOptions.py index 84f446c1b7cc5dbd3125b450d9dfee274d4998c6..5bec251868ec7e0c9c9f339f1cf88680891441c2 100755 --- a/TileCalorimeter/TileRec/share/TileDefaults_jobOptions.py +++ b/TileCalorimeter/TileRec/share/TileDefaults_jobOptions.py @@ -11,6 +11,7 @@ if ('doSim' in dir()) and doSim: doTileQIE = False doTileOptATLAS = False doTileMF = False + doTileWiener = False doTileFit = False doTileFitCool = False TileCisRun = False @@ -61,6 +62,11 @@ if not 'doTileMF' in dir(): else: jobproperties.TileRecFlags.doTileMF = doTileMF +if not 'doTileWiener' in dir(): + doTileWiener = jobproperties.TileRecFlags.doTileWiener() +else: + jobproperties.TileRecFlags.doTileWiener = doTileWiener + if not 'doTileFit' in dir(): doTileFit = jobproperties.TileRecFlags.doTileFit() else: diff --git a/TileCalorimeter/TileRec/share/TileNtuple_jobOptions.py b/TileCalorimeter/TileRec/share/TileNtuple_jobOptions.py index ebca0914cdd328c666d5fec366cc812d5295e1ed..46861e8bd2ca891023fa8a9a99209aae74d1b69e 100755 --- a/TileCalorimeter/TileRec/share/TileNtuple_jobOptions.py +++ b/TileCalorimeter/TileRec/share/TileNtuple_jobOptions.py @@ -18,6 +18,7 @@ TileNtuple.TileRawChannelContainerOpt = "" TileNtuple.TileRawChannelContainerQIE = "" TileNtuple.TileRawChannelContainerOF1 = "" TileNtuple.TileRawChannelContainerMF = "" +TileNtuple.TileRawChannelContainerWiener = "" TileNtuple.TileRawChannelContainerDsp = "" TileNtuple.TileLaserObject = "" TileNtuple.TileMuRcvRawChannelContainer= "" @@ -87,6 +88,9 @@ else: if doTileMF: TileNtuple.TileRawChannelContainerMF = "TileRawChannelMF" + if doTileWiener: + TileNtuple.TileRawChannelContainerWiener = "TileRawChannelWiener" + if useTMDB: TileNtuple.TileMuRcvRawChannelContainer = "MuRcvRawChCnt" TileNtuple.TileMuRcvDigitsContainer = "MuRcvDigitsCnt" diff --git a/TileCalorimeter/TileRec/src/TileAANtuple.cxx b/TileCalorimeter/TileRec/src/TileAANtuple.cxx index cf6b2e57a65406dbe3a2080121b2fdcfb05d4df9..18098aa5efd00da150b270fb736bfc1ce2d14419 100755 --- a/TileCalorimeter/TileRec/src/TileAANtuple.cxx +++ b/TileCalorimeter/TileRec/src/TileAANtuple.cxx @@ -149,6 +149,7 @@ TileAANtuple::TileAANtuple(std::string name, ISvcLocator* pSvcLocator) declareProperty("TileRawChannelContainerOF1", m_of1RawChannelContainerKey = ""); // declareProperty("TileRawChannelContainerDsp", m_dspRawChannelContainerKey = ""); // declareProperty("TileRawChannelContainerMF", m_mfRawChannelContainerKey = ""); // + declareProperty("TileRawChannelContainerWiener", m_wienerRawChannelContainerKey = "");// declareProperty("TileMuRcvRawChannelContainer", m_tileMuRcvRawChannelContainerKey = "MuRcvRawChCnt");// TMDB declareProperty("TileMuRcvDigitsContainer", m_tileMuRcvDigitsContainerKey = "MuRcvDigitsCnt");// TMDB declareProperty("TileMuRcvContainer", m_tileMuRcvContainerKey = "TileMuRcvCnt");// TMDB @@ -230,6 +231,7 @@ StatusCode TileAANtuple::initialize() { ATH_CHECK( m_qieRawChannelContainerKey.initialize(SG::AllowEmpty) ); ATH_CHECK( m_dspRawChannelContainerKey.initialize(SG::AllowEmpty) ); ATH_CHECK( m_of1RawChannelContainerKey.initialize(SG::AllowEmpty) ); + ATH_CHECK( m_wienerRawChannelContainerKey.initialize(SG::AllowEmpty) ); ATH_CHECK( m_l2CntKey.initialize(m_compareMode) ); ATH_MSG_INFO( "initialization completed" ) ; @@ -338,14 +340,15 @@ StatusCode TileAANtuple::execute() { // store TileRawChannels // start from DSP channels - so we can find out what is the DSP units - empty &= storeRawChannels(ctx, m_dspRawChannelContainerKey, m_arrays->m_eDsp, m_arrays->m_tDsp, m_arrays->m_chi2Dsp, m_arrays->m_pedDsp, true ).isFailure(); - empty &= storeRawChannels(ctx, m_rawChannelContainerKey, m_arrays->m_ene, m_arrays->m_time, m_arrays->m_chi2, m_arrays->m_ped, false).isFailure(); - empty &= storeMFRawChannels(ctx, m_mfRawChannelContainerKey, m_arrays->m_eMF, m_arrays->m_tMF, m_arrays->m_chi2MF, m_arrays->m_pedMF, false).isFailure(); - empty &= storeRawChannels(ctx, m_fitRawChannelContainerKey, m_arrays->m_eFit, m_arrays->m_tFit, m_arrays->m_chi2Fit, m_arrays->m_pedFit, false).isFailure(); - empty &= storeRawChannels(ctx, m_fitcRawChannelContainerKey, m_arrays->m_eFitc, m_arrays->m_tFitc, m_arrays->m_chi2Fitc,m_arrays->m_pedFitc,false).isFailure(); - empty &= storeRawChannels(ctx, m_optRawChannelContainerKey, m_arrays->m_eOpt, m_arrays->m_tOpt, m_arrays->m_chi2Opt, m_arrays->m_pedOpt, false).isFailure(); - empty &= storeRawChannels(ctx, m_qieRawChannelContainerKey, m_arrays->m_eQIE, m_arrays->m_tQIE, m_arrays->m_chi2QIE, m_arrays->m_pedQIE, false).isFailure(); - empty &= storeRawChannels(ctx, m_of1RawChannelContainerKey, m_arrays->m_eOF1, m_arrays->m_tOF1, m_arrays->m_chi2OF1, m_arrays->m_pedOF1, false).isFailure(); + empty &= storeRawChannels(ctx, m_dspRawChannelContainerKey, m_arrays->m_eDsp, m_arrays->m_tDsp, m_arrays->m_chi2Dsp, m_arrays->m_pedDsp, true ).isFailure(); + empty &= storeRawChannels(ctx, m_rawChannelContainerKey, m_arrays->m_ene, m_arrays->m_time, m_arrays->m_chi2, m_arrays->m_ped, false).isFailure(); + empty &= storeMFRawChannels(ctx, m_mfRawChannelContainerKey, m_arrays->m_eMF, m_arrays->m_tMF, m_arrays->m_chi2MF, m_arrays->m_pedMF, false).isFailure(); + empty &= storeRawChannels(ctx, m_fitRawChannelContainerKey, m_arrays->m_eFit, m_arrays->m_tFit, m_arrays->m_chi2Fit, m_arrays->m_pedFit, false).isFailure(); + empty &= storeRawChannels(ctx, m_fitcRawChannelContainerKey, m_arrays->m_eFitc, m_arrays->m_tFitc, m_arrays->m_chi2Fitc, m_arrays->m_pedFitc, false).isFailure(); + empty &= storeRawChannels(ctx, m_optRawChannelContainerKey, m_arrays->m_eOpt, m_arrays->m_tOpt, m_arrays->m_chi2Opt, m_arrays->m_pedOpt, false).isFailure(); + empty &= storeRawChannels(ctx, m_qieRawChannelContainerKey, m_arrays->m_eQIE, m_arrays->m_tQIE, m_arrays->m_chi2QIE, m_arrays->m_pedQIE, false).isFailure(); + empty &= storeRawChannels(ctx, m_of1RawChannelContainerKey, m_arrays->m_eOF1, m_arrays->m_tOF1, m_arrays->m_chi2OF1, m_arrays->m_pedOF1, false).isFailure(); + empty &= storeRawChannels(ctx, m_wienerRawChannelContainerKey, m_arrays->m_eWiener, m_arrays->m_tWiener, m_arrays->m_chi2Wiener, m_arrays->m_pedWiener, false).isFailure(); // store TMDB data // @@ -1827,6 +1830,7 @@ void TileAANtuple::DIGI_addBranch(void) || !m_dspRawChannelContainerKey.empty() || !m_mfRawChannelContainerKey.empty() || !m_of1RawChannelContainerKey.empty() + || !m_wienerRawChannelContainerKey.empty() || !m_bsInput) ) { m_ntuplePtr->Branch(NAME2("gain",f_suf), m_arrays->m_gain[ir], NAME3("gain", f_suf,"[4][64][48]/S")); // short @@ -1847,6 +1851,7 @@ void TileAANtuple::DIGI_addBranch(void) || !m_qieRawChannelContainerKey.empty() || !m_of1RawChannelContainerKey.empty() || !m_dspRawChannelContainerKey.empty() + || !m_wienerRawChannelContainerKey.empty() || m_bsInput) { m_ntuplePtr->Branch(NAME2("gain",f_suf), m_arrays->m_gain[ir], NAME3("gain", f_suf,"[4][64][48]/S")); // short @@ -1911,31 +1916,38 @@ void TileAANtuple::DIGI_addBranch(void) } if (!m_qieRawChannelContainerKey.empty()) { - m_ntuplePtr->Branch(NAME2("eQIE",f_suf), m_arrays->m_eQIE[ir], NAME3("eQIE",f_suf,"[4][64][48]/F")); // float - m_ntuplePtr->Branch(NAME2("tQIE",f_suf), m_arrays->m_tQIE[ir], NAME3("tQIE",f_suf,"[4][64][48]/F")); // float - m_ntuplePtr->Branch(NAME2("pedQIE",f_suf), m_arrays->m_pedQIE[ir], NAME3("pedQIE",f_suf,"[4][64][48]/F")); // float - m_ntuplePtr->Branch(NAME2("chi2QIE",f_suf), m_arrays->m_chi2QIE[ir], NAME3("chi2QIE",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("eQIE",f_suf), m_arrays->m_eQIE[ir], NAME3("eQIE",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("tQIE",f_suf), m_arrays->m_tQIE[ir], NAME3("tQIE",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("pedQIE",f_suf), m_arrays->m_pedQIE[ir], NAME3("pedQIE",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("chi2QIE",f_suf), m_arrays->m_chi2QIE[ir], NAME3("chi2QIE",f_suf,"[4][64][48]/F")); // float } if (!m_of1RawChannelContainerKey.empty()) { - m_ntuplePtr->Branch(NAME2("eOF1",f_suf), m_arrays->m_eOF1[ir], NAME3("eOF1",f_suf,"[4][64][48]/F")); // float - m_ntuplePtr->Branch(NAME2("tOF1",f_suf), m_arrays->m_tOF1[ir], NAME3("tOF1",f_suf,"[4][64][48]/F")); // float - m_ntuplePtr->Branch(NAME2("pedOF1",f_suf), m_arrays->m_pedOF1[ir], NAME3("pedOF1",f_suf,"[4][64][48]/F")); // float - m_ntuplePtr->Branch(NAME2("chi2OF1",f_suf), m_arrays->m_chi2OF1[ir], NAME3("chi2OF1",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("eOF1",f_suf), m_arrays->m_eOF1[ir], NAME3("eOF1",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("tOF1",f_suf), m_arrays->m_tOF1[ir], NAME3("tOF1",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("pedOF1",f_suf), m_arrays->m_pedOF1[ir], NAME3("pedOF1",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("chi2OF1",f_suf), m_arrays->m_chi2OF1[ir], NAME3("chi2OF1",f_suf,"[4][64][48]/F")); // float } if (!m_dspRawChannelContainerKey.empty()) { - m_ntuplePtr->Branch(NAME2("eDsp",f_suf), m_arrays->m_eDsp[ir], NAME3("eDsp",f_suf,"[4][64][48]/F")); // float - m_ntuplePtr->Branch(NAME2("tDsp",f_suf), m_arrays->m_tDsp[ir], NAME3("tDsp",f_suf,"[4][64][48]/F")); // float - m_ntuplePtr->Branch(NAME2("pedDsp",f_suf), m_arrays->m_pedDsp[ir], NAME3("pedDsp",f_suf,"[4][64][48]/F")); // float - m_ntuplePtr->Branch(NAME2("chi2Dsp",f_suf), m_arrays->m_chi2Dsp[ir], NAME3("chi2Dsp",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("eDsp",f_suf), m_arrays->m_eDsp[ir], NAME3("eDsp",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("tDsp",f_suf), m_arrays->m_tDsp[ir], NAME3("tDsp",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("pedDsp",f_suf), m_arrays->m_pedDsp[ir], NAME3("pedDsp",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("chi2Dsp",f_suf), m_arrays->m_chi2Dsp[ir], NAME3("chi2Dsp",f_suf,"[4][64][48]/F")); // float + } + + if (!m_wienerRawChannelContainerKey.empty()) { + m_ntuplePtr->Branch(NAME2("eWiener",f_suf), m_arrays->m_eWiener[ir], NAME3("eWiener",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("tWiener",f_suf), m_arrays->m_tWiener[ir], NAME3("tWiener",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("pedWiener",f_suf), m_arrays->m_pedWiener[ir], NAME3("pedWiener",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("chi2Wiener",f_suf), m_arrays->m_chi2Wiener[ir], NAME3("chi2Wiener",f_suf,"[4][64][48]/F")); // float } if (!m_mfRawChannelContainerKey.empty()) { - m_ntuplePtr->Branch(NAME2("eMF",f_suf), m_arrays->m_eMF[ir], NAME3("eMF",f_suf,"[4][64][48][7]/F")); // float - m_ntuplePtr->Branch(NAME2("tMF",f_suf), m_arrays->m_tMF[ir], NAME3("tMF",f_suf,"[4][64][48][7]/F")); // float - m_ntuplePtr->Branch(NAME2("chi2MF",f_suf), m_arrays->m_chi2MF[ir], NAME3("chi2MF",f_suf,"[4][64][48]/F")); // float - m_ntuplePtr->Branch(NAME2("pedMF",f_suf), m_arrays->m_pedMF[ir], NAME3("pedMF",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("eMF",f_suf), m_arrays->m_eMF[ir], NAME3("eMF",f_suf,"[4][64][48][7]/F")); // float + m_ntuplePtr->Branch(NAME2("tMF",f_suf), m_arrays->m_tMF[ir], NAME3("tMF",f_suf,"[4][64][48][7]/F")); // float + m_ntuplePtr->Branch(NAME2("chi2MF",f_suf), m_arrays->m_chi2MF[ir], NAME3("chi2MF",f_suf,"[4][64][48]/F")); // float + m_ntuplePtr->Branch(NAME2("pedMF",f_suf), m_arrays->m_pedMF[ir], NAME3("pedMF",f_suf,"[4][64][48]/F")); // float } if (m_bsInput) { @@ -2066,6 +2078,13 @@ void TileAANtuple::DIGI_clearBranch(void) { CLEAR2(m_arrays->m_chi2MF, size); CLEAR2(m_arrays->m_pedMF, size); } + + if (!m_wienerRawChannelContainerKey.empty()) { + CLEAR2(m_arrays->m_eWiener, size); + CLEAR2(m_arrays->m_tWiener, size); + CLEAR2(m_arrays->m_pedWiener, size); + CLEAR2(m_arrays->m_chi2Wiener, size); + } if (m_bsInput) { CLEAR1(m_arrays->m_ROD_GlobalCRC); diff --git a/TileCalorimeter/TileRecUtils/CMakeLists.txt b/TileCalorimeter/TileRecUtils/CMakeLists.txt index 6d193e46c8baba970f12120d76401c3a32cf4c27..4ef20b1858d191e38c28cfc2d754b69a902b53bc 100644 --- a/TileCalorimeter/TileRecUtils/CMakeLists.txt +++ b/TileCalorimeter/TileRecUtils/CMakeLists.txt @@ -29,7 +29,8 @@ atlas_depends_on_subdirs( PUBLIC Event/xAOD/xAODEventInfo TileCalorimeter/TileCalib/TileCalibBlobObjs TileCalorimeter/TileDetDescr - Tools/PathResolver ) + Tools/PathResolver + Trigger/TrigAnalysis/TrigAnalysisInterfaces ) # External dependencies: find_package( Boost COMPONENTS filesystem thread system ) diff --git a/TileCalorimeter/TileRecUtils/TileRecUtils/TileRawChannelBuilderWienerFilter.h b/TileCalorimeter/TileRecUtils/TileRecUtils/TileRawChannelBuilderWienerFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..51f66af48bbb97751bbe7cb5549ea742e37f2f45 --- /dev/null +++ b/TileCalorimeter/TileRecUtils/TileRecUtils/TileRawChannelBuilderWienerFilter.h @@ -0,0 +1,111 @@ +/* + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TILERECUTILS_TILERAWCHANNELBUILDERWIENERFILTER_H +#define TILERECUTILS_TILERAWCHANNELBUILDERWIENERFILTER_H + +////////////////////////////////////////////////////////////////////// +// +// AUTHOR : G. Goncalves <ginaciog@cern.ch> +// CREATED: May 2019 +// +// TileRawChannelBuilderWienerFilter.h +// +// implementation of the Wiener Filter for +// energy/time reconstruction in TileCal +// +////////////////////////////////////////////////////////////////////// + +// Tile includes +#include "TileRecUtils/TileRawChannelBuilder.h" +#include "TileConditions/ITileCondToolOfc.h" +#include "TileConditions/TileCondToolTiming.h" +#include "TileConditions/TileCondToolNoiseSample.h" + +// Atlas includes +#include "TrigAnalysisInterfaces/IBunchCrossingTool.h" + +#include <vector> +#include <string> + +/** + * + * @class TileRawChannelBuilderWienerFilter + * @brief Reconstructs Tile digitized pulses (ie, computes amplitude, time and pedestal) as a linear combination of the samples + * + * This class implements an energy reconstruction method known as Wiener + * Filter. It takes as input the digital samples from the digitizer boards + * in the front-end electronics and outputs the amplitude, time and pedestal + * of the pulse. Full details and fundaments of the method can be found in + * the ATL-TILECAL-PROC-2019-002 note. Two different versions of the algorithms + * are currently used: Wiener Optimal (with uses seven sets of weights, optimized + * for the first three, the middle and the last three BCIDs in the train) and + * Wiener General (one set of weights for all BCIDs). + */ +class TileRawChannelBuilderWienerFilter: public TileRawChannelBuilder { + public: + + TileRawChannelBuilderWienerFilter(const std::string& type, const std::string& name, + const IInterface *parent); //!< Constructor + ~TileRawChannelBuilderWienerFilter(); //!< Destructor + + // virtual methods + virtual StatusCode initialize(); //!< Initialize method + //virtual StatusCode execute(); + virtual StatusCode finalize(); //!< Finalize method + + // Inherited from TileRawChannelBuilder + virtual TileRawChannel* rawChannel(const TileDigits* digits); + + /** + * AlgTool InterfaceID + */ + static const InterfaceID& interfaceID(); + + private: + + ToolHandle<TileCondToolNoiseSample> m_tileToolNoiseSample{this, + "TileCondToolNoiseSample", "TileCondToolNoiseSample", "Tile noise sample tool"}; + + ToolHandle<Trig::IBunchCrossingTool> m_bunchCrossingTool; + + //!< Applies OF algorithm + double filter(int ros, int drawer, int channel, int &gain, double &pedestal, double &litude, double &time); + int findMaxDigitPosition(); //!< Finds maximum digit position in the pulse + //!< Gets pedestal estimation for OF1 + float getPedestal(int ros, int drawer, int channel, int gain); + //!< Computes A,time,ped using OF. If iterations are required, the Iterator method is used + double compute(int ros, int drawer, int channel, int gain, double &pedestal, double &litude, double &time, double& phase); + //!< Gets the BCID index within the train + int getBCIDIndex(); + + int m_maxIterations; //!< maximum number of iteration to perform + int m_pedestalMode; //!< pedestal mode to use + bool m_confTB; //!< use testbeam configuration + double m_timeForConvergence; //!< minimum time difference to quit iteration procedure + bool m_isMC; //!< bool variable for MC: true=> MC; false=> data + bool m_minus1Iter; //!< bool variable for whether to apply -1 iteration (initial phase guess) + bool m_correctAmplitude; //!< If true, resulting amplitude is corrected when using weights for tau=0 without iteration + bool m_correctTimeNI; //!< If true, resulting time is corrected when using method without iteration + + bool m_bestPhase; // if true, use best phase from COOL DB in "fixed phase" mode (i.e., no iterations) + bool m_emulateDsp; // if true, emulate DSP reconstruction algorithm + int m_nSignal; //!< internal counters + int m_nNegative; //!< internal counters + int m_nCenter; //!< internal counters + int m_nConst; //!< internal counters + + int m_nSamples; //!< number of samples in the data + int m_t0SamplePosition; //!< position of peak sample = (m_nSamples-1)/2 + double m_maxTime; //!< max allowed time = 25*(m_nSamples-1)/2 + double m_minTime; //!< min allowed time = -25*(m_nSamples-1)/2 + + + std::vector<float> m_digits; + + static const float m_gfcWiener[8]; //!< Wiener General weights + static const float m_ofcWiener[7][8]; //!< Wiener Optimal weights +}; + +#endif diff --git a/TileCalorimeter/TileRecUtils/doc/packagedoc.h b/TileCalorimeter/TileRecUtils/doc/packagedoc.h index 023f44a613dc50272c04d1b5c0b8b18439318986..0433ebbe8d3eef4197e7c4e70787f79c5168e978 100644 --- a/TileCalorimeter/TileRecUtils/doc/packagedoc.h +++ b/TileCalorimeter/TileRecUtils/doc/packagedoc.h @@ -37,6 +37,8 @@ This package contains the algorithms and tools used for TileCal reconstruction - TileRawChannelBuilderManyAmps: used by TileDigitsMaker. TileRawChannels are built using the method developed by Frank Marritt, based on fitting the digital samples to to a series of pulses (for pileup handling). Default method for full ATLAS reconstruction. + - TileRawChannelBuilderWienerFilter: used by TileDigitsMaker. TileRawChannels are built using the Wiener Filtering method, based on the cross-correlation between the samples and the noise minimizing its contribution. + - TileCellBuilder: creates TileCells from TileRawChannels, which are stored in a container. By default uses Opt2 Optimal Filter, and on-th-fly cell masking. - TileCellFakeProb: scales down the energy of TileCells due to simulated failure of drawers diff --git a/TileCalorimeter/TileRecUtils/python/TileRawChannelGetter.py b/TileCalorimeter/TileRecUtils/python/TileRawChannelGetter.py index f2d0da7be47987d04d0b6117656036e3ce06d74c..53d60801804876979bc0bffa31087a59d99cab68 100644 --- a/TileCalorimeter/TileRecUtils/python/TileRawChannelGetter.py +++ b/TileCalorimeter/TileRecUtils/python/TileRawChannelGetter.py @@ -156,6 +156,7 @@ class TileRawChannelGetter ( Configured) : if (jobproperties.TileRecFlags.doTileMF() or (not jobproperties.TileRecFlags.OfcFromCOOL() and (jobproperties.TileRecFlags.doTileOF1() + or jobproperties.TileRecFlags.doTileWiener() or jobproperties.TileRecFlags.doTileOpt2() or jobproperties.TileRecFlags.doTileOptATLAS()))): @@ -167,6 +168,7 @@ class TileRawChannelGetter ( Configured) : if jobproperties.TileRecFlags.OfcFromCOOL(): if (jobproperties.TileRecFlags.doTileMF() + or jobproperties.TileRecFlags.doTileWiener() or jobproperties.TileRecFlags.doTileOpt2() or jobproperties.TileRecFlags.doTileOptATLAS()): @@ -510,6 +512,35 @@ class TileRawChannelGetter ( Configured) : theTileRawChannelMaker.TileRawChannelBuilder += [theTileRawChannelBuilderOptATLAS] self._TileRawChannelBuilderOptATLAS = theTileRawChannelBuilderOptATLAS + + if jobproperties.TileRecFlags.doTileWiener(): + try: + from TileRecUtils.TileRecUtilsConf import TileRawChannelBuilderWienerFilter + theTileRawChannelBuilderWienerFilter= TileRawChannelBuilderWienerFilter() + except Exception: + mlog.error("could not get handle to TileRawChannelBuilderWienerFilter Quit") + traceback.print_exc() + return False + + #TileRawChannelBuilderWienerFilter Options: + jobproperties.TileRecFlags.TileRawChannelContainer = "TileRawChannelWiener" + theTileRawChannelBuilderWienerFilter.TileRawChannelContainer = "TileRawChannelWiener" + theTileRawChannelBuilderWienerFilter.RunType = jobproperties.TileRecFlags.TileRunType() + theTileRawChannelBuilderWienerFilter.calibrateEnergy = jobproperties.TileRecFlags.calibrateEnergy() + theTileRawChannelBuilderWienerFilter.correctTime = jobproperties.TileRecFlags.correctTime() + theTileRawChannelBuilderWienerFilter.NoiseFilterTools = NoiseFilterTools + theTileRawChannelBuilderWienerFilter.BestPhase = False # no point to use best phase with interations + theTileRawChannelBuilderWienerFilter.MC = globalflags.DataSource()!='data' + theTileRawChannelBuilderWienerFilter.PedestalMode = 1 + theTileRawChannelBuilderWienerFilter.MaxIterations = 5 + theTileRawChannelBuilderWienerFilter.Minus1Iteration = True + theTileRawChannelBuilderWienerFilter.AmplitudeCorrection = False # don't need correction after iterations + theTileRawChannelBuilderWienerFilter.TimeCorrection = False # don't need correction after iterations + theTileRawChannelBuilderWienerFilter.DSPContainer = TileRawChannelContainerDSP + + mlog.info(" adding now TileRawChannelBuilderWienerFilter to the algorithm: %s", theTileRawChannelMaker.name()) + + theTileRawChannelMaker.TileRawChannelBuilder += [theTileRawChannelBuilderWienerFilter] # now add algorithm to topSequence @@ -547,6 +578,7 @@ class TileRawChannelGetter ( Configured) : jobproperties.TileRecFlags.doTileManyAmps = False jobproperties.TileRecFlags.doTileMF = False jobproperties.TileRecFlags.doTileOF1 = False + jobproperties.TileRecFlags.doTileWiener = False jobproperties.TileRecFlags.OfcFromCOOL = False jobproperties.TileRecFlags.print_JobProperties('tree&value') diff --git a/TileCalorimeter/TileRecUtils/python/TileRawChannelGetter_DigiHSTruth.py b/TileCalorimeter/TileRecUtils/python/TileRawChannelGetter_DigiHSTruth.py index d7db276ef175c247d012d7c3f09dc94e94cfb1b7..9d9b9cac4c1e1cea70f05cbc5c4ca52e0cebd2bc 100644 --- a/TileCalorimeter/TileRecUtils/python/TileRawChannelGetter_DigiHSTruth.py +++ b/TileCalorimeter/TileRecUtils/python/TileRawChannelGetter_DigiHSTruth.py @@ -158,6 +158,7 @@ class TileRawChannelGetter_DigiHSTruth ( Configured) : jobproperties.TileRecFlags.doTileManyAmps = False jobproperties.TileRecFlags.doTileMF = False jobproperties.TileRecFlags.doTileOF1 = False + jobproperties.TileRecFlags.doTileWiener = False jobproperties.TileRecFlags.OfcFromCOOL = False jobproperties.TileRecFlags.print_JobProperties('tree&value') diff --git a/TileCalorimeter/TileRecUtils/python/TileRecFlags.py b/TileCalorimeter/TileRecUtils/python/TileRecFlags.py index 8f085d31df1447807431bdca5b0b96f04e05eb69..670eda6eff58461bc2b06f9b7746186c90c236a8 100644 --- a/TileCalorimeter/TileRecUtils/python/TileRecFlags.py +++ b/TileCalorimeter/TileRecUtils/python/TileRecFlags.py @@ -96,6 +96,14 @@ class doTileMF(JobProperty): allowedTypes = ['bool'] StoredValue = False +# +class doTileWiener(JobProperty): + """ Use Wiener filter method for energy reconstruction + """ + statusOn = True + allowedTypes = ['bool'] + StoredValue = False + # class noiseFilter(JobProperty): """ appply one or another noise suppression algorithm at the level of raw channels @@ -266,6 +274,7 @@ list_jobproperties = [ doTileManyAmps, doTileMF, doTileOptATLAS, + doTileWiener, TileRawChannelContainer, TileDigitsContainer, TileRunType, diff --git a/TileCalorimeter/TileRecUtils/share/TileCellBuilder_test.py b/TileCalorimeter/TileRecUtils/share/TileCellBuilder_test.py index d4cdcdba668799c4eb3eef2da4f9eab20e7683de..d2ff6eafaaaffb4c49d0dc8def944f28f86bae74 100644 --- a/TileCalorimeter/TileRecUtils/share/TileCellBuilder_test.py +++ b/TileCalorimeter/TileRecUtils/share/TileCellBuilder_test.py @@ -261,6 +261,7 @@ class TileFragHash: FitFilter = 6 FitFilterCool = 7 FlatFilter = 8 + WienerFilterOffline = 9 from AthenaPython.PyAthenaComps import Alg, StatusCode diff --git a/TileCalorimeter/TileRecUtils/share/TileDQstatusAlg_test.py b/TileCalorimeter/TileRecUtils/share/TileDQstatusAlg_test.py index 68a1ebd256f93a37a1a6b5092f29db95d052fa5c..126fcf97625a33b8b1b5412936769c5bdfdc4b80 100644 --- a/TileCalorimeter/TileRecUtils/share/TileDQstatusAlg_test.py +++ b/TileCalorimeter/TileRecUtils/share/TileDQstatusAlg_test.py @@ -62,6 +62,7 @@ class TileFragHash: FitFilter = 6 FitFilterCool = 7 FlatFilter = 8 + WienerFilterOffline = 9 BEAM_TDC_FRAG = 0x000 diff --git a/TileCalorimeter/TileRecUtils/share/TileDQstatusTool_test.py b/TileCalorimeter/TileRecUtils/share/TileDQstatusTool_test.py index 4ad299625abf26fdaa8a899817497853efd22436..bf6681b6d63aad2f906f502e7cafdf75e44b88ef 100644 --- a/TileCalorimeter/TileRecUtils/share/TileDQstatusTool_test.py +++ b/TileCalorimeter/TileRecUtils/share/TileDQstatusTool_test.py @@ -62,6 +62,7 @@ class TileFragHash: FitFilter = 6 FitFilterCool = 7 FlatFilter = 8 + WienerFilterOffline = 9 BEAM_TDC_FRAG = 0x000 diff --git a/TileCalorimeter/TileRecUtils/share/TileRawChannelBuilder_test.py b/TileCalorimeter/TileRecUtils/share/TileRawChannelBuilder_test.py index b3e1e0fbf13e7ebb2b61fa52490377549a68c268..166f8a1b392a4ce50970e9479285b799bcba6a07 100644 --- a/TileCalorimeter/TileRecUtils/share/TileRawChannelBuilder_test.py +++ b/TileCalorimeter/TileRecUtils/share/TileRawChannelBuilder_test.py @@ -190,6 +190,7 @@ class TileFragHash: FitFilter = 6 FitFilterCool = 7 FlatFilter = 8 + WienerFilterOffline = 9 from AthenaPython.PyAthenaComps import Alg, StatusCode diff --git a/TileCalorimeter/TileRecUtils/src/TileRawChannelBuilderWienerFilter.cxx b/TileCalorimeter/TileRecUtils/src/TileRawChannelBuilderWienerFilter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..1d67b4a1b2021b4b0d644cb669dd9b71d0cc98e6 --- /dev/null +++ b/TileCalorimeter/TileRecUtils/src/TileRawChannelBuilderWienerFilter.cxx @@ -0,0 +1,442 @@ +/* + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +*/ + +#include "TileEvent/TileRawChannel.h" +#include "TileCalibBlobObjs/TileCalibUtils.h" + +// Gaudi includes +#include "GaudiKernel/Property.h" +#include "GeoModelInterfaces/IGeoModelSvc.h" + +// Atlas includes +#include "AthAllocators/DataPool.h" +#include "AthenaKernel/errorcheck.h" +#include "xAODEventInfo/EventInfo.h" + +// Tile includes +#include "TileRecUtils/TileRawChannelBuilderWienerFilter.h" +#include "TileConditions/TilePulseShapes.h" +#include "CLHEP/Matrix/Matrix.h" +#include "TileEvent/TileRawChannelContainer.h" +#include "TileEvent/TileDigitsContainer.h" +#include "TileEvent/TileDigits.h" +#include "CaloIdentifier/TileID.h" +#include "TileIdentifier/TileHWID.h" +#include "TileConditions/TileInfo.h" + +//using namespace std; +#include <algorithm> + +//interface stuff +static const InterfaceID IID_ITileRawChannelBuilderWienerFilter("TileRawChannelBuilderWienerFilter", 1, 0); + + +const InterfaceID& TileRawChannelBuilderWienerFilter::interfaceID() { + return IID_ITileRawChannelBuilderWienerFilter; +} + + +#define TILE_WienerFilterBUILDERVERBOSE false + +// initialize Wiener coefficients ( ZeroBias <mu>=40 ) +const float TileRawChannelBuilderWienerFilter::m_gfcWiener[8] = { + -0.194081, + 0.467562, + -1.25846, + 2.18528, + -1.2648, + 0.63341, + -0.276375, + -14.2093 +}; +const float TileRawChannelBuilderWienerFilter::m_ofcWiener[7][8] = { + {0.04382, -0.999, 0.5212, 1.179, -0.8052, 0.4123, -0.1756, -5.393}, + {0.1806, 0.02539, -1.253, 2.187, -1.254, 0.6214, -0.2732, -14.53}, + {-0.3267, 0.6014, -1.346, 2.224, -1.258, 0.6043, -0.276, -14.41}, + {-0.1914, 0.4693, -1.246, 2.169, -1.243, 0.6191, -0.2633, -15.15}, + {-0.1831, 0.4734, -1.253, 2.241, -1.424, 0.9173, -0.5364, -15.43}, + {-0.2245, 0.5627, -1.398, 2.433, -1.722, 1.204, -0.5617, -15.26}, + {-0.2367, 0.617, -1.502, 2.539, -1.729, 0.7092, -0.1435, -14.15} +}; + +TileRawChannelBuilderWienerFilter::TileRawChannelBuilderWienerFilter(const std::string& type, + const std::string& name, const IInterface *parent) + : TileRawChannelBuilder(type, name, parent) + , m_bunchCrossingTool("BunchCrossingTool") + , m_nSignal(0) + , m_nNegative(0) + , m_nCenter(0) + , m_nConst(0) + , m_nSamples(0) + , m_t0SamplePosition(0) + , m_maxTime(0.0) + , m_minTime(0.0) +{ + //declare interfaces + declareInterface< TileRawChannelBuilder >( this ); + declareInterface< TileRawChannelBuilderWienerFilter >(this); + + m_rawChannelContainerKey = "TileRawChannelWiener"; + + //declare properties + declareProperty("BunchCrossingTool", m_bunchCrossingTool); + declareProperty("MaxIterations",m_maxIterations = 5); + declareProperty("PedestalMode",m_pedestalMode = 17); + declareProperty("TimeForConvergence",m_timeForConvergence = 0.5); + declareProperty("ConfTB",m_confTB = false); + declareProperty("Minus1Iteration",m_minus1Iter = false); + declareProperty("AmplitudeCorrection",m_correctAmplitude = false); + declareProperty("TimeCorrection", m_correctTimeNI = false); + declareProperty("BestPhase",m_bestPhase = false); + declareProperty("EmulateDSP",m_emulateDsp = false); + declareProperty("MC",m_isMC = false); +} + + +TileRawChannelBuilderWienerFilter::~TileRawChannelBuilderWienerFilter() { +} + + +StatusCode TileRawChannelBuilderWienerFilter::initialize() { + + ATH_MSG_INFO( "initialize()" ); + + m_rChType = TileFragHash::WienerFilterOffline; // type for offline Wiener Filter + + // init in superclass + ATH_CHECK( TileRawChannelBuilder::initialize() ); + + if (m_maxIterations != 1) m_correctTimeNI = false; + + // bits 12-15 - various options + // if (m_correctTimeNI) m_bsflags |= 0x1000; + if (m_correctAmplitude) m_bsflags |= 0x2000; + if (m_maxIterations > 1) m_bsflags |= 0x4000; + if (m_bestPhase) m_bsflags |= 0x8000; + + ATH_MSG_DEBUG( " MaxIterations=" << m_maxIterations + << " PedestalMode=" << m_pedestalMode + << " TimeForConvergence=" << m_timeForConvergence + << " ConfTB=" << m_confTB + << " Minus1Iteration=" << m_minus1Iter + << " AmplitudeCorrection=" << m_correctAmplitude + << " TimeCorrection=" << m_correctTimeNI + << " Best Phase " << m_bestPhase ); + + m_nSamples = m_tileInfo->NdigitSamples(); + m_t0SamplePosition = m_tileInfo->ItrigSample(); + m_maxTime = 25 * (m_nSamples - m_t0SamplePosition - 1); + m_minTime = -25 * m_t0SamplePosition; + ATH_MSG_DEBUG(" NSamples=" << m_nSamples + << " T0Sample=" << m_t0SamplePosition + << " minTime=" << m_minTime + << " maxTime=" << m_maxTime ); + + if (m_pedestalMode % 10 > 2 && m_nSamples != m_pedestalMode % 10) { + if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Changing PedestalMode from " << m_pedestalMode; + m_pedestalMode = (m_pedestalMode / 10) * 10 + m_nSamples; + if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " to " << m_pedestalMode << endmsg; + } + + if (m_nSamples != 7 && (m_pedestalMode == 71 || m_pedestalMode == 7621)) { + ATH_MSG_ERROR( "Incompatable pedestal mode [" << m_pedestalMode + << "] and number of samples [" << m_nSamples << "]" ); + return StatusCode::FAILURE; + } + + m_nSignal = 0; + m_nNegative = 0; + m_nCenter = 0; + m_nConst = 0; + + //=== get TileCondToolNoiseSample + ATH_CHECK( m_tileToolNoiseSample.retrieve() ); + + if (m_bestPhase) { + //=== get TileToolTiming + // TileToolTiming can be disabled in the TileRawChannelBuilder + if (!m_tileToolTiming.isEnabled()) { + m_tileToolTiming.enable(); + } + ATH_CHECK( m_tileToolTiming.retrieve() ); + } + + //=== get TrigBunchCrossingTool + if (m_isMC) { + m_bunchCrossingTool.setTypeAndName("Trig::MCBunchCrossingTool/BunchCrossingTool"); + } else { // is data + m_bunchCrossingTool.setTypeAndName("Trig::LHCBunchCrossingTool/BunchCrossingTool"); + } + ATH_CHECK( m_bunchCrossingTool.retrieve() ); + + ATH_MSG_INFO( "initialization completed" ); + + return StatusCode::SUCCESS; +} + +StatusCode TileRawChannelBuilderWienerFilter::finalize() { + + if (msgLvl(MSG::VERBOSE)) { + msg(MSG::VERBOSE) << "Counters: Signal=" << m_nSignal + << " Constant=" << m_nConst + << " Total=" << m_nSignal + m_nConst << endmsg; + } + + ATH_MSG_DEBUG( "Finalizing" ); + + return StatusCode::SUCCESS; +} + + +TileRawChannel * TileRawChannelBuilderWienerFilter::rawChannel(const TileDigits* digits) { + + ++m_chCounter; + + double pedestal = 0.; + double energy = 0.; + double time = 0.; + double chi2 = 0.; + m_digits = digits->samples(); + m_digits.erase(m_digits.begin(),m_digits.begin()+m_firstSample); + m_digits.resize(m_nSamples); + const HWIdentifier adcId = digits->adc_HWID(); + int gain = m_tileHWID->adc(adcId); + + ATH_MSG_VERBOSE( "Building Raw Channel, with WienerFilter, HWID:" << m_tileHWID->to_string(adcId) + << " gain=" << gain ); + + int ros = m_tileHWID->ros(adcId); + int drawer = m_tileHWID->drawer(adcId); + int channel = m_tileHWID->channel(adcId); + + chi2 = filter(ros, drawer, channel, gain, pedestal, energy, time); + + unsigned int drawerIdx = TileCalibUtils::getDrawerIdx(ros, drawer); + + if (m_calibrateEnergy) { + energy = m_tileToolEmscale->doCalibCis(drawerIdx, channel, gain, energy); + } + + if (msgLvl(MSG::VERBOSE)) { + msg(MSG::VERBOSE) << "Creating WienerFilter RawChannel" + << " a=" << energy + << " t=" << time + << " ped=" << pedestal + << " q=" << chi2 << endmsg; + + msg(MSG::VERBOSE) << "digits:"; + + for (unsigned int i = 0; i < m_digits.size(); ++i) + msg(MSG::VERBOSE) << " " << m_digits[i]; + + msg(MSG::VERBOSE) << " " << endmsg; + } + + DataPool<TileRawChannel> tileRchPool(m_dataPoollSize); + TileRawChannel *rawCh = tileRchPool.nextElementPtr(); + rawCh->assign (adcId, + energy, + time, + chi2, + pedestal); + + if (m_correctTime + && (time != 0 + && time < m_maxTime + && time > m_minTime)) { + + time -= m_tileToolTiming->getSignalPhase(drawerIdx, channel, gain); + rawCh->insertTime(time); + + ATH_MSG_VERBOSE( "Correcting time, new time=" << rawCh->time() ); + + } + + if (TileID::HIGHGAIN == gain) { + ++m_nChH; + m_RChSumH += energy; + } else { + ++m_nChL; + m_RChSumL += energy; + } + + return rawCh; +} + + +int TileRawChannelBuilderWienerFilter::findMaxDigitPosition() { + + ATH_MSG_VERBOSE( " findMaxDigitPosition()" ); + + int iMaxDigit = 0; + float maxDigit = 0.; + bool saturated = false; + + for (unsigned int i = 0; i < m_digits.size(); i++) { + if (m_digits[i] > 1022.99) saturated = true; + if (maxDigit < m_digits[i]) { + maxDigit = m_digits[i]; + iMaxDigit = i; + } + } + + if (msgLvl(MSG::VERBOSE)) { + for (unsigned int i = 0; i < m_digits.size(); i++) { + msg(MSG::VERBOSE) << " " << m_digits[i]; + } + + msg(MSG::VERBOSE) << "; Max: digit[" << iMaxDigit << "]=" << maxDigit << endmsg; + + if (saturated) msg(MSG::VERBOSE) << " Samples saturated" << endmsg; + } + + return iMaxDigit; +} + + +float TileRawChannelBuilderWienerFilter::getPedestal(int ros, int drawer, int channel, int gain) { + float pedestal = 0.; + + switch (m_pedestalMode) { + case -1: + // use pedestal from conditions DB + pedestal = m_tileToolNoiseSample->getPed(TileCalibUtils::getDrawerIdx(ros, drawer), channel, gain); + break; + case 7: + pedestal = m_digits[6]; + break; + case 9: + pedestal = m_digits[8]; + break; + case 12: + pedestal = .5 * (m_digits[0] + m_digits[1]); + break; + case 17: + pedestal = .5 * (m_digits[0] + m_digits[6]); + break; + case 19: + pedestal = .5 * (m_digits[0] + m_digits[8]); + break; + case 71: + pedestal = std::min(m_digits[0], m_digits[6]); + break; + case 7621: + pedestal = 0.5 * std::min(m_digits[0] + m_digits[1], m_digits[5] + m_digits[6]); + break; + default: + pedestal = m_digits[0]; + break; + } + + ATH_MSG_VERBOSE("getPedestal(): pedestal=" << pedestal); + + return pedestal; +} + +int TileRawChannelBuilderWienerFilter::getBCIDIndex() { + int bcidIndex = -1; + + const xAOD::EventInfo* eventInfo(0); + + if (evtStore()->retrieve(eventInfo).isSuccess()) { + int bcid = eventInfo->bcid(); + int distFront = m_bunchCrossingTool->distanceFromFront(bcid, Trig::IBunchCrossingTool::BunchCrossings); + int distTail = m_bunchCrossingTool->distanceFromTail(bcid, Trig::IBunchCrossingTool::BunchCrossings); + + ATH_MSG_VERBOSE( "EventInfo loaded: " + << " BCID=" << bcid + << " DistFront=" << distFront + << " DistTail=" << distTail); + + if (distFront == -1 || distTail == -1) { + bcidIndex = -1; + } else if (distFront < 3 && distTail > distFront) { + bcidIndex = distFront; + } else if (distTail < 3 && distTail < distFront) { + bcidIndex = 6 - distTail; + } else { + bcidIndex = 3; + } + } else { + ATH_MSG_VERBOSE("EventInfo not available"); + } + + ATH_MSG_VERBOSE("getBCIDIndex(): bcidIndex=" << bcidIndex); + + return bcidIndex; +} + + +double TileRawChannelBuilderWienerFilter::filter(int ros, int drawer, int channel + , int &gain, double &pedestal, double &litude, double &time) { + + ATH_MSG_VERBOSE( "filter()" ); + + amplitude = 0.; + time = 0.; + double chi2 = 0.; + + auto minMaxDigits = std::minmax_element(m_digits.begin(), m_digits.end()); + float minDigit = *minMaxDigits.first; + float maxDigit = *minMaxDigits.second; + + if (maxDigit - minDigit < 0.01) { // constant value in all samples + + pedestal = minDigit; + chi2 = 0.; + ATH_MSG_VERBOSE( "CASE NO SIGNAL: maxdig-mindig = " << maxDigit << "-" << minDigit + << " = " << maxDigit - minDigit ); + + m_nConst++; + + } else { + float weights[8]; + memset(weights, 0, sizeof(weights)); + + int bcidIndex = getBCIDIndex(); + + if (ros > 1 && channel == 1 && bcidIndex >= 0) { + ATH_MSG_VERBOSE( "Switch to optimal mode"); + for (int wi = 0; wi < 8; wi++) { + weights[wi] = m_ofcWiener[bcidIndex][wi]; + } + } else { + ATH_MSG_VERBOSE( "Switch to general mode"); + for (int wi = 0; wi < 8; wi++) { + weights[wi] = m_gfcWiener[wi]; + } + } + + for (unsigned int i = 0; i < m_digits.size(); i++) { + amplitude += weights[i] * m_digits[i]; // Wiener coefs + } + amplitude += weights[7]; // Wiener bias + + double phase = 0.; + pedestal = getPedestal(ros, drawer, channel, gain); + + chi2 = compute(ros, drawer, channel, gain, pedestal, amplitude, time, phase); + + m_nSignal++; + } + + return chi2; +} + +double TileRawChannelBuilderWienerFilter::compute(int ros, int drawer, int channel, int gain, + double &pedestal, double &litude, double &time, double& phase) { + + ATH_MSG_VERBOSE( "compute();" + << " ros=" << ros + << " drawer=" << drawer + << " channel=" << channel + << " gain=" << gain + << " pedestal=" << pedestal + << " amplitude=" << amplitude + << " time=" << time + << " phase=" << phase ); + + double chi2 = 0.; + return chi2; +} diff --git a/TileCalorimeter/TileRecUtils/src/components/TileRecUtils_entries.cxx b/TileCalorimeter/TileRecUtils/src/components/TileRecUtils_entries.cxx index fa0c5cccd0f6c29f1fb2a216af9ff32446b53a7c..db545ac7fa8b97226b4898ea4876d03186f1452d 100644 --- a/TileCalorimeter/TileRecUtils/src/components/TileRecUtils_entries.cxx +++ b/TileCalorimeter/TileRecUtils/src/components/TileRecUtils_entries.cxx @@ -6,6 +6,7 @@ #include "TileRecUtils/TileRawChannelBuilderQIEFilter.h" #include "TileRecUtils/TileRawChannelBuilderManyAmps.h" #include "TileRecUtils/TileRawChannelBuilderMF.h" +#include "TileRecUtils/TileRawChannelBuilderWienerFilter.h" #include "TileRecUtils/TileCellBuilder.h" #include "TileRecUtils/TileCellBuilderFromHit.h" #include "TileRecUtils/TileCellFakeProb.h" @@ -36,6 +37,7 @@ DECLARE_COMPONENT( TileRawChannelBuilderOpt2Filter ) DECLARE_COMPONENT( TileRawChannelBuilderQIEFilter ) DECLARE_COMPONENT( TileRawChannelBuilderManyAmps ) DECLARE_COMPONENT( TileRawChannelBuilderMF ) +DECLARE_COMPONENT( TileRawChannelBuilderWienerFilter ) DECLARE_COMPONENT( TileCellBuilder ) DECLARE_COMPONENT( TileCellBuilderFromHit ) DECLARE_COMPONENT( TileCellFakeProb )