ATLAS Offline Software
LArDigitOscillationCorrTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 
7 #include "GaudiKernel/MsgStream.h"
8 #include "GaudiKernel/IIncidentSvc.h"
9 #include "CLHEP/Units/SystemOfUnits.h"
10 
15 
16 #include <sstream>
17 #include <iostream>
18 #include <fstream>
19 #include <iomanip>
20 #include <cmath>
21 
22 
23 using CLHEP::hertz;
24 using CLHEP::ns;
25 
26 
28  const std::string& name,
29  const IInterface* parent)
30  : base_class(type, name, parent),
31  m_priority(1400),
32  m_nSigma(3.0),
33  m_eventPhase(0),
34  m_omega(1.024e6*hertz)
35 {
36  declareProperty("BeginRunPriority",m_priority);
37  declareProperty("SignalCutInSigma",m_nSigma);
38  declareProperty("Omega",m_omega);
39 }
40 
42 {
43  ATH_MSG_DEBUG ( "LArDigitOscillationCorrTool initialize() begin" );
44 
45  ATH_CHECK( detStore()->retrieve(m_lar_on_id,"LArOnlineID") );
46 
47  const CaloCell_ID* idHelper = nullptr;
48  ATH_CHECK( detStore()->retrieve (idHelper, "CaloCell_ID") );
49  m_emId = idHelper->em_idHelper();
50  m_fcalId = idHelper->fcal_idHelper();
51  m_hecId = idHelper->hec_idHelper();
52 
54 
55  SmartIF<IIncidentSvc> incSvc{Gaudi::svcLocator()->service("IncidentSvc")};
56  ATH_CHECK(incSvc.isValid());
57 
58  //start listening to "BeginRun"
59  incSvc->addListener(this, "BeginRun", m_priority);
60 
61  ATH_MSG_DEBUG ( "LArDigitOscillationCorrTool initialize() end" );
62  return StatusCode::SUCCESS;
63 }
64 
66 {
67  ATH_MSG_DEBUG ( "In execute: calculate event phase" );
68 
69  //Pointer to conditions data objects to get pedetestals
70  ATH_MSG_DEBUG ( "Retrieving pedestal " );
71  const ILArPedestal* larPedestal=nullptr;
72  ATH_CHECK( detStore()->retrieve(larPedestal) );
73 
74  //Pointer to conditions data objects to get channel phases
75  ATH_MSG_DEBUG ( "Retrieving channel phases and amplituides " );
76  const ILArH6Oscillation* larH6Oscillations = nullptr;
77  ATH_CHECK( detStore()->retrieve(larH6Oscillations) );
78  unsigned int iHECChan(0);
79  std::vector<std::vector<short> > theSamples;
80  std::vector<float> thePedestals;
81  std::vector<double> theChannelAmplitudes;
82  std::vector<double> theChannelPhases;
83  std::vector<float> theRMSValues;
84 
85  //Now all data is available, start loop over Digit Container
86  ATH_MSG_DEBUG ( "Loop over Digit Container with size <" << theDC.size() << ">" );
87 
89  const LArOnOffIdMapping* cabling{*cablingHdl};
90  if(!cabling) {
91  ATH_MSG_ERROR("Do not have mapping object " << m_cablingKey.key() );
92  return StatusCode::FAILURE;
93  }
94 
95 
96  for (unsigned int i=0;i<theDC.size();i++) {
97  //Get data from LArDigit
98  const LArDigit *theDigit = theDC[i];
99  const std::vector<short>& samples = theDigit->samples();
100  const HWIdentifier chid = theDigit->channelID();
101  const CaloGain::CaloGain gain = theDigit->gain();
102 
103  // use only HEC channels
104  const Identifier id = cabling->cnvToIdentifier(chid);
105  bool isHEC = false;
106 
107  if (m_hecId->is_lar_hec(id)) {
108  isHEC = true;
109  }
110 
111  if ( isHEC) {
112  float DBpedestalRMS=larPedestal->pedestalRMS(chid,gain);
113  if (DBpedestalRMS <= (1.0+LArElecCalib::ERRORCODE)) {
114  ATH_MSG_DEBUG ( "No pedestal RMS found for this cell. Exiting ...." );
115  return StatusCode::FAILURE;
116  }
117  // log << MSG::DEBUG << "Retriving channelPhase " << endmsg;
118  const double& DBchannelPhase=larH6Oscillations->channelPhase(chid);
119  // log << MSG::DEBUG << "Retriving channelAmplitude " << endmsg;
120  const double& DBchannelAmplitude=larH6Oscillations->channelAmplitude(chid);
121 
122  if( DBpedestalRMS > 0 && DBchannelAmplitude>0 ) {
123  theRMSValues.push_back(DBpedestalRMS);
124  theSamples.push_back(samples);
125  theChannelPhases.push_back(DBchannelPhase);
126  theChannelAmplitudes.push_back(DBchannelAmplitude);
127  float DBpedestal=larPedestal->pedestal(chid,gain);
128  if (DBpedestal >= (1.0+LArElecCalib::ERRORCODE))
129  thePedestals.push_back(DBpedestal);
130  else {
131  ATH_MSG_DEBUG ( "No valid pedestal found for this cell. Exiting ...." );
132  return StatusCode::FAILURE;
133  }
134  iHECChan++;
135  }
136  }
137  }
138 
139  ATH_MSG_DEBUG ( "Start calculating the event phase" );
140  double lchimin = -1;
141  int nTotSamples = 0;
142  for(double myEventPhase = -M_PI; myEventPhase<M_PI; myEventPhase+=0.1) {
143  double lchitest = 0;
144  for (unsigned int i=0;i<iHECChan;i++) {
145  unsigned int nSamples = theSamples[i].size();
146  for(unsigned int j=0;j<nSamples;j++) {
147  // exclude all samples (and the previous and next sample) which
148  //have more than nSigma worth of absolute signal
149  if ( j == 0 || (std::abs(theSamples[i][j] - theSamples[i][0]) <= m_nSigma*theRMSValues[i]
150  && ( j <= 0 || std::abs(theSamples[i][j-1] - theSamples[i][0]) <= m_nSigma*theRMSValues[i])
151  && ( j >= nSamples-1 || std::abs(theSamples[i][j+1] - theSamples[i][0]) <= m_nSigma*theRMSValues[i])) ) {
152 
153  if ( lchimin < 0 )
154  nTotSamples ++;
155  lchitest += pow((theSamples[i][j] - thePedestals[i]
156  - theChannelAmplitudes[i]*sin(j*25*ns*m_omega
157  +theChannelPhases[i]+myEventPhase))
158  /theRMSValues[i],2);
159  }
160  }
161  }
162  if ( lchitest < lchimin || lchimin < 0 ) {
163  lchimin = lchitest;
164  m_eventPhase = myEventPhase;
165  }
166  }
167 
168  ATH_MSG_DEBUG ( "Ending eventphase calculation. Phase = <" << m_eventPhase << ">, Number of Channels used = <" << iHECChan << ">, Average Number of Samples Used = <"
169  << (iHECChan > 0 ? (double)nTotSamples/iHECChan : 0)
170  << ">" );
171  return StatusCode::SUCCESS;
172 }
173 
174 
176 {
177  ATH_MSG_DEBUG ( "In execute: correct LAr Digit" );
178 
179  //Pointer to conditions data objects to get channel phases
180  const ILArH6Oscillation* larH6Oscillations = nullptr;
181  ATH_CHECK( detStore()->retrieve(larH6Oscillations) );
182 
183  std::vector<double> theChannelAmplitudes;
184  std::vector<double> theChannelPhases;
185 
186  //Now all data is available, start loop over Digit Container
187  ATH_MSG_DEBUG ( "Loop over Digit Container " );
188 
190  const LArOnOffIdMapping* cabling{*cablingHdl};
191  if(!cabling) {
192  ATH_MSG_ERROR("Do not have mapping object " << m_cablingKey.key() );
193  return StatusCode::FAILURE;
194  }
195 
196  for (unsigned int i=0;i<theDC.size();i++) {
197  //Get data from LArDigit
198  LArDigit * theDigit = theDC[i];
199  const std::vector<short>& samples = theDigit->samples();
200  unsigned int nSamples = samples.size();
201  const HWIdentifier chid = theDigit->channelID();
202  const CaloGain::CaloGain gain = theDigit->gain();
203 
204  // correct only HEC channels
205  const Identifier id = cabling->cnvToIdentifier(chid);
206  bool isHEC = false;
207 
208  if (m_hecId->is_lar_hec(id)) {
209  isHEC = true;
210  }
211 
212  if ( isHEC) {
213  const double& DBchannelPhase=larH6Oscillations->channelPhase(chid);
214  const double& DBchannelAmplitude=larH6Oscillations->channelAmplitude(chid);
215  ATH_MSG_DEBUG ( "The HWId is " << chid
216  << ", the offline Id is " << m_hecId->show_to_string(id,nullptr,'/')
217  << ", the ChannelAmplitude is " << DBchannelAmplitude
218  << ", the ChannelPhase is " << DBchannelPhase );
219  ATH_MSG_DEBUG ( "m_omega value is " << m_omega );
220  if( DBchannelAmplitude>0 ) {
221  std::vector<short> new_samples(nSamples);
222  for(unsigned int j=0;j<nSamples;j++)
223  new_samples[j] = (short)(samples[j] + 0.5
224  - DBchannelAmplitude*sin(j*25*ns*m_omega+DBchannelPhase+m_eventPhase));
225  HWIdentifier new_id(chid);
226  CaloGain::CaloGain new_gain(gain);
227  LArDigit * theNewDigit = new LArDigit(new_id,new_gain,new_samples);
228  theDC[i] = theNewDigit;
229  }
230  }
231  }
232  return StatusCode::SUCCESS;
233 }
234 
235 
236 void LArDigitOscillationCorrTool::handle(const Incident& /* inc*/ )
237 {
238  ATH_MSG_DEBUG ( "LArDigitOscillationCorrTool handle()" );
239 
240  this->retrieveDB().ignore();
241 }
242 
243 // *** retrieve subfactors from the DB ***
245 {
246  ATH_MSG_DEBUG ( "in retrieveDB() " );
247  return StatusCode::SUCCESS;
248 }
249 
250 
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
ILArPedestal::pedestal
virtual float pedestal(const HWIdentifier &id, int gain) const =0
LArDigitOscillationCorrTool::correctLArDigits
StatusCode correctLArDigits(LArDigitContainer &theDC) override
Definition: LArDigitOscillationCorrTool.cxx:175
LArDigitOscillationCorrTool::m_priority
int m_priority
Definition: LArDigitOscillationCorrTool.h:52
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
LArDigitOscillationCorrTool::retrieveDB
StatusCode retrieveDB()
Definition: LArDigitOscillationCorrTool.cxx:244
CaloCell_ID::em_idHelper
const LArEM_ID * em_idHelper() const
access to EM idHelper
Definition: CaloCell_ID.h:63
LArDigitOscillationCorrTool::m_lar_on_id
const LArOnlineID * m_lar_on_id
Definition: LArDigitOscillationCorrTool.h:60
ILArPedestal
Definition: ILArPedestal.h:12
CaloCondBlobAlgs_fillNoiseFromASCII.gain
gain
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:110
ReadCellNoiseFromCool.cabling
cabling
Definition: ReadCellNoiseFromCool.py:154
LArDigit::samples
const std::vector< short > & samples() const
Definition: LArDigit.h:78
M_PI
#define M_PI
Definition: ActiveFraction.h:11
LArDigitOscillationCorrTool::LArDigitOscillationCorrTool
LArDigitOscillationCorrTool(const std::string &type, const std::string &name, const IInterface *parent)
Definition: LArDigitOscillationCorrTool.cxx:27
LArDigitOscillationCorrTool::m_hecId
const LArHEC_ID * m_hecId
Definition: LArDigitOscillationCorrTool.h:59
SG::VarHandleKey::key
const std::string & key() const
Return the StoreGate ID for the referenced object.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:141
HWIdentifier
Definition: HWIdentifier.h:13
LArDigitOscillationCorrTool::initialize
virtual StatusCode initialize() override
Definition: LArDigitOscillationCorrTool.cxx:41
LArDigitOscillationCorrTool::m_eventPhase
double m_eventPhase
Definition: LArDigitOscillationCorrTool.h:54
LArDigitOscillationCorrTool::m_nSigma
double m_nSigma
Definition: LArDigitOscillationCorrTool.h:53
CaloCell_ID.h
CaloCell_ID::hec_idHelper
const LArHEC_ID * hec_idHelper() const
access to HEC idHelper
Definition: CaloCell_ID.h:69
python.SystemOfUnits.hertz
int hertz
Definition: SystemOfUnits.py:125
LArDigitOscillationCorrTool::m_fcalId
const LArFCAL_ID * m_fcalId
Definition: LArDigitOscillationCorrTool.h:58
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
LArDigit
Liquid Argon digit base class.
Definition: LArDigit.h:25
lumiFormat.i
int i
Definition: lumiFormat.py:85
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
LArDigitOscillationCorrTool::m_omega
double m_omega
Definition: LArDigitOscillationCorrTool.h:55
test_pyathena.parent
parent
Definition: test_pyathena.py:15
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
AtlasDetectorID::is_lar_hec
bool is_lar_hec(Identifier id) const
Definition: AtlasDetectorID.h:829
CaloCell_ID
Helper class for offline cell identifiers.
Definition: CaloCell_ID.h:34
ILArH6Oscillation.h
python.PyKernel.detStore
detStore
Definition: PyKernel.py:41
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
ILArH6Oscillation::channelAmplitude
virtual const double & channelAmplitude(const HWIdentifier &id) const =0
access to channel amplitude index by Identifier
SG::CondHandleKey::initialize
StatusCode initialize(bool used=true)
LArDigitContainer.h
CaloGain::CaloGain
CaloGain
Definition: CaloGain.h:11
LArDigit::gain
CaloGain::CaloGain gain() const
Definition: LArDigit.h:72
LArDigitOscillationCorrTool::handle
virtual void handle(const Incident &) override
Definition: LArDigitOscillationCorrTool.cxx:236
AtlasDetectorID::show_to_string
std::string show_to_string(Identifier id, const IdContext *context=0, char sep='.') const
or provide the printout in string form
Definition: AtlasDetectorID.cxx:574
LArDigitOscillationCorrTool::calculateEventPhase
StatusCode calculateEventPhase(const LArDigitContainer &theDC) override
Definition: LArDigitOscillationCorrTool.cxx:65
ILArPedestal.h
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
ILArH6Oscillation::channelPhase
virtual const double & channelPhase(const HWIdentifier &id) const =0
LArDigitContainer
Container class for LArDigit.
Definition: LArDigitContainer.h:24
python.SystemOfUnits.ns
int ns
Definition: SystemOfUnits.py:130
LArDigit::channelID
const HWIdentifier & channelID() const
Definition: LArDigit.h:69
LArElecCalib::ERRORCODE
@ ERRORCODE
Definition: LArCalibErrorCode.h:17
LArDigits2NtupleDumper.nSamples
nSamples
Definition: LArDigits2NtupleDumper.py:70
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
LArDigitOscillationCorrTool::m_cablingKey
SG::ReadCondHandleKey< LArOnOffIdMapping > m_cablingKey
Definition: LArDigitOscillationCorrTool.h:62
ILArH6Oscillation
Definition: ILArH6Oscillation.h:14
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
ILArPedestal::pedestalRMS
virtual float pedestalRMS(const HWIdentifier &id, int gain) const =0
access to RMS of Pedestal index by Identifier, and gain setting
LArDigitOscillationCorrTool.h
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
LArDigitOscillationCorrTool::m_emId
const LArEM_ID * m_emId
Definition: LArDigitOscillationCorrTool.h:57
CaloCell_ID::fcal_idHelper
const LArFCAL_ID * fcal_idHelper() const
access to FCAL idHelper
Definition: CaloCell_ID.h:75
LArOnOffIdMapping
Definition: LArOnOffIdMapping.h:20
Identifier
Definition: IdentifierFieldParser.cxx:14