CMS 3D CMS Logo
TPedValues.cc
Go to the documentation of this file.
1 #include "CalibCalorimetry/EcalPedestalOffsets/interface/TPedValues.h"
2 #include "FWCore/MessageLogger/interface/MessageLogger.h"
3 
4 #include "TAxis.h"
5 #include "TF1.h"
6 #include "TGraphErrors.h"
7 #include <cassert>
8 #include <cmath>
9 #include <iostream>
10 
11 void reset(double vett[256]) {
12  for (int i = 0; i < 256; ++i)
13  vett[i] = 0.;
14 }
15 
16 TPedValues::TPedValues(double RMSmax, int bestPedestal) : m_bestPedestal(bestPedestal), m_RMSmax(RMSmax) {
17  LogDebug("EcalPedOffset") << "entering TPedValues ctor ...";
18  for (int i = 0; i < 1700; ++i)
19  endcapCrystalNumbers[i] = 0;
20 }
21 
22 TPedValues::TPedValues(const TPedValues &orig) {
23  LogDebug("EcalPedOffset") << "entering TPedValues copyctor ...";
24  m_bestPedestal = orig.m_bestPedestal;
25  m_RMSmax = orig.m_RMSmax;
26 
27  for (int gain = 0; gain < 3; ++gain)
28  for (int crystal = 0; crystal < 1700; ++crystal)
29  for (int DAC = 0; DAC < 256; ++DAC)
30  m_entries[gain][crystal][DAC] = orig.m_entries[gain][crystal][DAC];
31 
32  for (int i = 0; i < 1700; ++i)
33  endcapCrystalNumbers[i] = orig.endcapCrystalNumbers[i];
34 }
35 
36 TPedValues::~TPedValues() {}
37 
38 void TPedValues::insert(const int gainId, const int crystal, const int DAC, const int pedestal, const int endcapIndex) {
39  // assert (gainId > 0) ;
40  // assert (gainId < 4) ;
41  if (gainId <= 0 || gainId >= 4) {
42  edm::LogWarning("EcalPedOffset") << "WARNING : TPedValues : gainId " << gainId << " does not exist, entry skipped";
43  return;
44  }
45  // assert (crystal > 0) ;
46  // assert (crystal <= 1700) ;
47  if (crystal <= 0 || crystal > 1700) {
48  edm::LogWarning("EcalPedOffset") << "WARNING : TPedValues : crystal " << crystal
49  << " does not exist, entry skipped";
50  return;
51  }
52  // assert (DAC >= 0) ;
53  // assert (DAC < 256) ;
54  if (DAC < 0 || DAC >= 256) {
55  edm::LogWarning("EcalPedOffset") << "WARNING : TPedValues : DAC value " << DAC << " is out range, entry skipped";
56  return;
57  }
58  m_entries[gainId - 1][crystal - 1][DAC].insert(pedestal);
59  endcapCrystalNumbers[crystal - 1] = endcapIndex;
60  return;
61 }
62 
63 TPedResult TPedValues::terminate(const int &DACstart, const int &DACend) const {
64  assert(DACstart >= 0);
65  assert(DACend <= 256);
66  // checkEntries (DACstart, DACend) ;
67 
68  TPedResult bestDAC;
70  for (int gainId = 1; gainId < 4; ++gainId) {
72  for (int crystal = 0; crystal < 1700; ++crystal) {
74  double delta = 1000;
75  int dummyBestDAC = -1;
76  bool hasDigis = false;
78  for (int DAC = DACstart; DAC < DACend; ++DAC) {
79  double average = m_entries[gainId - 1][crystal][DAC].average();
80  if (average == -1)
81  continue;
82  hasDigis = true;
83  if (m_entries[gainId - 1][crystal][DAC].RMSSq() > m_RMSmax * m_RMSmax)
84  continue;
85  if (fabs(average - m_bestPedestal) < delta && average > 1) {
86  delta = fabs(average - m_bestPedestal);
87  dummyBestDAC = DAC;
88  }
89  }
90 
91  bestDAC.m_DACvalue[gainId - 1][crystal] = dummyBestDAC;
92 
93  if ((dummyBestDAC == (DACend - 1) || dummyBestDAC == -1) && hasDigis) {
94  int gainHuman;
95  if (gainId == 1)
96  gainHuman = 12;
97  else if (gainId == 2)
98  gainHuman = 6;
99  else if (gainId == 3)
100  gainHuman = 1;
101  else
102  gainHuman = -1;
103  edm::LogError("EcalPedOffset") << " TPedValues : cannot find best DAC value for channel: "
104  << endcapCrystalNumbers[crystal] << " gain: " << gainHuman;
105  }
106 
107  } // loop over crystals
108  } // loop over gains
109  return bestDAC;
110 }
111 
112 int TPedValues::checkEntries(const int &DACstart, const int &DACend) const {
113  assert(DACstart >= 0);
114  assert(DACend <= 256);
115  int returnCode = 0;
117  for (int gainId = 1; gainId < 4; ++gainId) {
119  for (int crystal = 0; crystal < 1700; ++crystal) {
121  for (int DAC = DACstart; DAC < DACend; ++DAC) {
122  double average = m_entries[gainId - 1][crystal][DAC].average();
123  if (average == -1) {
124  ++returnCode;
126  /*
127  std::cerr << "[TPedValues][checkEntries] WARNING!"
128  << "\tgainId " << gainId
129  << "\tcrystal " << crystal+1
130  << "\tDAC " << DAC
131  << " : pedestal measurement missing"
132  << std::endl ;
133  */
134  }
135  /*
136  std::cout << "[pietro][RMS]: " <<
137  m_entries[gainId-1][crystal][DAC].RMS () //FIXME
138  << "\t" <<
139  m_entries[gainId-1][crystal][DAC].RMSSq () //FIXME
140  << "\t" << DAC //FIXME
141  << "\t" << gainId //FIXME
142  << "\t" << crystal << std::endl ; //FIXME
143  */
144  }
145  } // loop over crystals
146  } // loop over gains
147  return returnCode;
148 }
149 
151 int TPedValues::makePlots(TFile *rootFile,
152  const std::string &dirName,
153  const double maxSlope,
154  const double minSlope,
155  const double maxChi2OverNDF) const {
156  using namespace std;
157  // prepare the ROOT file
158  if (!rootFile->cd(dirName.c_str())) {
159  rootFile->mkdir(dirName.c_str());
160  rootFile->cd(dirName.c_str());
161  }
162 
163  // loop over the crystals
164  for (int xtl = 0; xtl < 1700; ++xtl) {
165  // loop over the gains
166  for (int gain = 0; gain < 3; ++gain) {
167  vector<double> asseX;
168  vector<double> sigmaX;
169  vector<double> asseY;
170  vector<double> sigmaY;
171  asseX.reserve(256);
172  sigmaX.reserve(256);
173  asseY.reserve(256);
174  sigmaY.reserve(256);
175  // loop over DAC values
176  for (int dac = 0; dac < 256; ++dac) {
177  double average = m_entries[gain][xtl][dac].average();
178  if (average > -1) {
179  double rms = m_entries[gain][xtl][dac].RMS();
180  asseX.push_back(dac);
181  sigmaX.push_back(0);
182  asseY.push_back(average);
183  sigmaY.push_back(rms);
184  }
185  } // loop over DAC values
186  if (!asseX.empty()) {
187  int lastBin = 0;
188  while (lastBin < (int)asseX.size() - 1 && asseY[lastBin + 1] > 0 &&
189  (asseY[lastBin + 1] - asseY[lastBin + 2]) != 0)
190  lastBin++;
191 
192  int fitRangeEnd = (int)asseX[lastBin];
193  int kinkPt = 64;
194  if (fitRangeEnd < 66)
195  kinkPt = fitRangeEnd - 4;
196  TGraphErrors graph(asseX.size(), &(*asseX.begin()), &(*asseY.begin()), &(*sigmaX.begin()), &(*sigmaY.begin()));
197  char funct[120];
198  sprintf(funct, "(x<%d)*([0]*x+[1])+(x>=%d)*([2]*x+[3])", kinkPt, kinkPt);
199  TF1 fitFunction("fitFunction", funct, asseX[0], fitRangeEnd);
200  fitFunction.SetLineColor(2);
201 
202  char name[120];
203  int gainHuman;
204  if (gain == 0)
205  gainHuman = 12;
206  else if (gain == 1)
207  gainHuman = 6;
208  else if (gain == 2)
209  gainHuman = 1;
210  else
211  gainHuman = -1;
212  sprintf(name, "XTL%04d_GAIN%02d", endcapCrystalNumbers[xtl], gainHuman);
213  graph.GetXaxis()->SetTitle("DAC value");
214  graph.GetYaxis()->SetTitle("Average pedestal ADC");
215  graph.Fit(&fitFunction, "RWQ");
216  graph.Write(name);
217 
218  double slope1 = fitFunction.GetParameter(0);
219  double slope2 = fitFunction.GetParameter(2);
220 
221  if (fitFunction.GetChisquare() / fitFunction.GetNDF() > maxChi2OverNDF ||
222  fitFunction.GetChisquare() / fitFunction.GetNDF() < 0 || slope1 > 0 || slope2 > 0 ||
223  ((slope1 < -29 || slope1 > -18) && slope1 < 0) || ((slope2 < -29 || slope2 > -18) && slope2 < 0)) {
224  edm::LogError("EcalPedOffset") << "TPedValues : TGraph for channel:" << endcapCrystalNumbers[xtl]
225  << " gain:" << gainHuman << " is not linear;"
226  << " slope of line1:" << fitFunction.GetParameter(0)
227  << " slope of line2:" << fitFunction.GetParameter(2) << " reduced chi-squared:"
228  << fitFunction.GetChisquare() / fitFunction.GetNDF();
229  }
230  // LogDebug("EcalPedOffset") << "TPedValues : TGraph for channel:" <<
231  // xtl+1 << " gain:"
232  // << gainHuman << " has " << asseX.size() << " points...back is:" <<
233  // asseX.back()
234  // << " and front+1 is:" << asseX.front()+1;
235  if ((asseX.back() - asseX.front() + 1) != asseX.size())
236  edm::LogError("EcalPedOffset") << "TPedValues : Pedestal average not found "
237  << "for all DAC values scanned in channel:" << endcapCrystalNumbers[xtl]
238  << " gain:" << gainHuman;
239  }
240  } // loop over the gains
241  } // (loop over the crystals)
242 
243  return 0;
244 }
245 
246 // Look up the crystal number in the EE schema and return it
247 int TPedValues::getCrystalNumber(int xtal) const { return endcapCrystalNumbers[xtal]; }
ecalLiteDTU::gainId
constexpr int gainId(sample_type sample)
get the gainId (2 bits)
Definition: EcalLiteDTUSample.h:14
mps_fire.i
i
Definition: mps_fire.py:428
BeamSpotPI::sigmaX
Definition: BeamSpotPayloadInspectorHelper.h:34
MessageLogger.h
TPedValues::terminate
TPedResult terminate(const int &DACstart=0, const int &DACend=256) const
calculate the offset values for all the containers
Definition: TPedValues.cc:63
TPedResult
Definition: TPedResult.h:15
TSinglePedEntry::RMS
double RMS() const
get the RMS of the inserted values
Definition: TSinglePedEntry.cc:30
cms::cuda::assert
assert(be >=bs)
TPedValues::insert
void insert(const int gainId, const int crystal, const int DAC, const int pedestal, const int endcapIndex)
add a single value
Definition: TPedValues.cc:38
TPedValues::m_bestPedestal
int m_bestPedestal
Definition: TPedValues.h:52
SiStripPI::rms
Definition: SiStripPayloadInspectorHelper.h:169
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
TPedValues::TPedValues
TPedValues(double RMSmax=2, int bestPedestal=200)
ctor
Definition: TPedValues.cc:16
indexGen.rootFile
rootFile
Definition: indexGen.py:92
TPedValues
Definition: TPedValues.h:19
BeamSpotPI::sigmaY
Definition: BeamSpotPayloadInspectorHelper.h:35
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
TPedValues::endcapCrystalNumbers
int endcapCrystalNumbers[1700]
Definition: TPedValues.h:56
dumpMFGeometry_cfg.delta
delta
Definition: dumpMFGeometry_cfg.py:25
TPedValues::makePlots
int makePlots(TFile *rootFile, const std::string &dirName, const double maxSlopeAllowed, const double minSlopeAllowed, const double maxChi2OverNDF) const
create a plot of the DAC pedestal trend
Definition: TPedValues.cc:151
createfilelist.int
int
Definition: createfilelist.py:10
runEdmFileComparison.returnCode
returnCode
Definition: runEdmFileComparison.py:263
average
Definition: average.py:1
EcalCondDBWriter_cfi.pedestal
pedestal
Definition: EcalCondDBWriter_cfi.py:49
edm::LogError
Log< level::Error, false > LogError
Definition: MessageLogger.h:123
TPedValues::m_entries
TSinglePedEntry m_entries[3][1700][256]
Definition: TPedValues.h:50
cms::cuda::for
for(int i=first, nt=offsets[nh];i< nt;i+=gridDim.x *blockDim.x)
Definition: HistoContainer.h:27
TPedValues::getCrystalNumber
int getCrystalNumber(int xtal) const
return the index from the table
Definition: TPedValues.cc:247
PedestalClient_cfi.gain
gain
Definition: PedestalClient_cfi.py:37
TPedValues::checkEntries
int checkEntries(const int &DACstart=0, const int &DACend=256) const
check whether the scan was complete in the range
Definition: TPedValues.cc:112
TPedValues::~TPedValues
~TPedValues()
dtor
Definition: TPedValues.cc:36
std
Definition: JetResolutionObject.h:76
TPedValues.h
Transient container Store all the pedestal values depending on the gain and pedestal offset $Date: $R...
TSinglePedEntry::insert
void insert(const int &pedestal)
add a single value
Definition: TSinglePedEntry.cc:18
TrackerOfflineValidation_Dqm_cff.dirName
dirName
Definition: TrackerOfflineValidation_Dqm_cff.py:55
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
reset
void reset(double vett[256])
Definition: TPedValues.cc:11
TPedResult::m_DACvalue
int m_DACvalue[3][1700]
Definition: TPedResult.h:25
edm::Log
Definition: MessageLogger.h:70
TPedValues::m_RMSmax
double m_RMSmax
Definition: TPedValues.h:53
funct
Definition: Abs.h:5
TSinglePedEntry::average
double average() const
get the average of the inserted values
Definition: TSinglePedEntry.cc:24
Generated for CMSSW Reference Manual by doxygen 1.8.16