ATLAS Offline Software
Loading...
Searching...
No Matches
ZDCNLCalibration Class Reference

#include <ZDCNLCalibration.h>

Collaboration diagram for ZDCNLCalibration:

Public Member Functions

 ZDCNLCalibration (const std::string &file, int maxNLPower=3, bool useGRL=true, int debugLevel=0)
virtual ~ZDCNLCalibration ()
std::array< float, 4 > FindSNPeaks (size_t LBLow, size_t LBHigh, size_t side)
std::pair< float, float > FindSNRange (size_t LBLow, size_t LBHigh, size_t side)
std::pair< std::pair< float, float >, std::pair< float, float > > FindSNTwoNRanges (size_t LBLow, size_t LBHigh, size_t side)
void SetDefaultCalibration (size_t side, const CalibData &calib)
void Calibrate (size_t side, const std::string &calibInput, const std::string &calibOutput, size_t LBLow, size_t LBHigh, std::array< int, 4 > maxPowerModule, const std::vector< std::pair< double, double > > &nNeutERange, bool excludeHE, float heSumThresh, float HEDeweight)
void TestCalibration (int side, const std::string &calibName)
TH1 * GetTestSNHist ()
TH1 * GetTestFracHist (size_t module)
TTree * GetTestTree ()

Public Attributes

std::array< std::map< std::string, CalibData >, 2 > m_calibrations
bool m_haveTest {}
TH1D * m_testCalibSNHist {}
std::array< TH1D *, 4 > m_testCalibHEFracHist
std::array< TH1D *, 4 > m_testCalibEnergyHist
TTree * m_testTree {}

Private Types

typedef std::multimap< unsigned int, std::pair< unsigned int, unsigned int > > LBEvtMap

Private Member Functions

void FillLumiBlockEvtMap ()
void FillMinimizationData (TMatrixD &minimMatrix, TVectorD &minimVector, std::array< int, 4 > maxPowerModule, float HEDeweight, const std::vector< std::vector< double > > &sums1DVec, const std::vector< double > &sumsHE, const std::vector< std::vector< double > > &sums2DVec, const std::vector< double > &sumsHE2D)
void AddCalibration (size_t side, const std::string &tag, const CalibData &calib)
CalibData GetCalibration (size_t side, const std::string &tag)
void AddToSums (std::vector< double > &sums1D, std::vector< double > &sums2D, float *amps)

Static Private Member Functions

static double CalculateEnergy (const float *moduleAmps, const CalibData &calib)

Private Attributes

TFile * m_TFile
TTree * m_tree
size_t m_maxNLPower
bool m_useGRL
int m_debugLevel
const float m_SNEnergy
const std::vector< float > m_HEFraction
UInt_t runNumber
UInt_t eventNumber
UInt_t lumiBlock
UInt_t bcid
UInt_t passBits
Float_t zdc_ZdcAmp [2]
UInt_t zdc_ZdcModuleMask
Float_t zdc_ZdcModuleAmp [2][4]
Bool_t L1_ZDC_A
Bool_t L1_ZDC_C
Bool_t L1_ZDC_AND
Bool_t L1_ZDC_A_C
TBranch * b_runNumber
TBranch * b_eventNumber
TBranch * b_lumiBlock
TBranch * b_bcid
TBranch * b_passBits
TBranch * b_zdc_ZdcAmp
TBranch * b_zdc_ZdcModuleMask
TBranch * b_zdc_ZdcModuleAmp
TBranch * b_L1_ZDC_A
TBranch * b_L1_ZDC_C
TBranch * b_L1_ZDC_AND
TBranch * b_L1_ZDC_A_C
LBEvtMap m_LumiBlockEvtMap

Detailed Description

Definition at line 55 of file ZDCNLCalibration.h.

Member Typedef Documentation

◆ LBEvtMap

typedef std::multimap<unsigned int, std::pair<unsigned int, unsigned int> > ZDCNLCalibration::LBEvtMap
private

Definition at line 107 of file ZDCNLCalibration.h.

Constructor & Destructor Documentation

◆ ZDCNLCalibration()

ZDCNLCalibration::ZDCNLCalibration ( const std::string & file,
int maxNLPower = 3,
bool useGRL = true,
int debugLevel = 0 )

◆ ~ZDCNLCalibration()

virtual ZDCNLCalibration::~ZDCNLCalibration ( )
inlinevirtual

Definition at line 122 of file ZDCNLCalibration.h.

122{}

Member Function Documentation

◆ AddCalibration()

void ZDCNLCalibration::AddCalibration ( size_t side,
const std::string & tag,
const CalibData & calib )
private

Definition at line 51 of file ZDCNLCalibration.cxx.

52{
53 //cppcheck-suppress mismatchingContainers
54 std::map<std::string, CalibData>::iterator iter = m_calibrations.at(side).find(tag);
55 if (iter != m_calibrations[side].end()) (*iter).second = calib;//*iter = std::pair<std::string, CalibData>(tag, calib);
56 else {
57 m_calibrations[side].insert(std::pair<std::string, CalibData>(tag, calib));
58 }
59}
std::array< std::map< std::string, CalibData >, 2 > m_calibrations

◆ AddToSums()

void ZDCNLCalibration::AddToSums ( std::vector< double > & sums1D,
std::vector< double > & sums2D,
float * amps )
inlineprivate

Definition at line 161 of file ZDCNLCalibration.h.

162 {
163 size_t index1D = 0;
164 size_t index2D = 0;
165
166 for (size_t module1 : {0, 1, 2, 3}) {
167 double amp1 = amps[module1];
168 double amp1Pow = amp1;
169
170 for (size_t power1 = 0; power1 < m_maxNLPower; power1++) {
171 for (size_t module2 : {0, 1, 2, 3}) {
172 double amp2 = amps[module2];
173 double amp2Pow = amp2;
174
175 for (size_t power2 = 0; power2 < m_maxNLPower; power2++) {
176 sums2D.at(index2D++) += amp1Pow*amp2Pow;
177 amp2Pow *= amp2;
178 }
179 }
180
181 sums1D.at(index1D++) += amp1Pow;
182 amp1Pow *= amp1;
183 }
184 }
185 }

◆ CalculateEnergy()

double ZDCNLCalibration::CalculateEnergy ( const float * moduleAmps,
const CalibData & calib )
inlinestaticprivate

Definition at line 187 of file ZDCNLCalibration.h.

188 {
189 double energy = 0;
190
191 for (size_t module : {0, 1, 2, 3}) {
192 float amp = moduleAmps[module];
193 float ampPow = amp;
194
195 for (size_t power = 0; power < calib.weights[module].size(); power++) {
196 energy += calib.weights[module][power]*ampPow;
197 ampPow *= amp;
198 }
199 }
200
201 return energy;
202 }

◆ Calibrate()

void ZDCNLCalibration::Calibrate ( size_t side,
const std::string & calibInput,
const std::string & calibOutput,
size_t LBLow,
size_t LBHigh,
std::array< int, 4 > maxPowerModule,
const std::vector< std::pair< double, double > > & nNeutERange,
bool excludeHE,
float heSumThresh,
float HEDeweight )

Definition at line 197 of file ZDCNLCalibration.cxx.

201{
202 std::for_each(maxPowerModule.begin(), maxPowerModule.end(),
203 [](int power) {
204 if (power == 0) {
205 std::cout << "Found max power = 0, quitting" << std::endl;
206 return;
207 }
208 } );
209
210 // Set up container for sums that we will use to calcualting the weights
211 //
212 size_t numWeights = std::accumulate(maxPowerModule.begin(), maxPowerModule.end(), 0);
213 size_t arrayDim = 4*m_maxNLPower;
214
215 size_t arrayDim2D = arrayDim*arrayDim;
216
217 size_t numNeutMax = nNeutERange.size();
218
219 std::vector<double> sumsHE(arrayDim, 0);
220 std::vector<double> sums2DHE(arrayDim2D, 0);
221
222 std::vector<std::vector<double> > sumsNeutVec(numNeutMax, std::vector<double>(arrayDim, 0));
223 std::vector<std::vector<double> > sums2DNeutVec(numNeutMax, std::vector<double>(arrayDim2D, 0));
224
225 LBEvtMap::const_iterator LBStart = m_LumiBlockEvtMap.lower_bound(LBLow);
226 LBEvtMap::const_iterator LBEnd = m_LumiBlockEvtMap.upper_bound(LBHigh);
227
228 // Select which trigger we're going to look at based on the side
229 //
230 const bool& triggerOpp = (side == 1 ? L1_ZDC_C : L1_ZDC_A);
231 const bool& triggerSame = (side == 1 ? L1_ZDC_A : L1_ZDC_C);
232
233 std::vector<int> numNeutEvent(numNeutMax, 0);
234 int numHEEvent = 0;
235
236 //
237 //
238 std::string calibName = (calibInput != "" ? calibInput : "default");
239 CalibData calib = GetCalibration(side, calibName);
240
241 for ( LBEvtMap::const_iterator iter = LBStart; iter != LBEnd; ++iter) {
242 unsigned int entryStart = iter->second.first;
243 unsigned int entryEnd = iter->second.second;
244
245 for (unsigned int entry = entryStart; entry < entryEnd; ++entry) {
246 Long64_t nb = m_tree->GetEntry(entry);
247 if (nb < 0) {
248 std::cout << "Error reading entry " << entry << ", quitting" << std::endl;
249 throw std::runtime_error("Error reading entry for ZDCNLCalibration");
250 }
251
252 int testMask = zdc_ZdcModuleMask >> side*4;
253 if ((testMask & 0xf) != 0xf) continue;
254
255 double sumAmp = CalculateEnergy(zdc_ZdcModuleAmp[side], calib);
256
257 if (triggerSame && sumAmp > heSumThresh && !excludeHE) {
258 AddToSums(sumsHE, sums2DHE, zdc_ZdcModuleAmp[side]);
259 numHEEvent++;
260 }
261
262 // For 1N we require the opposite side trigger
263 //
264 if (triggerOpp && sumAmp > nNeutERange[0].first && sumAmp < nNeutERange[0].second) {
265 AddToSums(sumsNeutVec[0], sums2DNeutVec[0], zdc_ZdcModuleAmp[side]);
266 numNeutEvent[0]++;
267 }
268
269 for (size_t iNNeut = 1; iNNeut < numNeutMax; iNNeut++) {
270 if (sumAmp > nNeutERange[iNNeut].first && sumAmp < nNeutERange[iNNeut].second) {
271 AddToSums(sumsNeutVec[iNNeut], sums2DNeutVec[iNNeut], zdc_ZdcModuleAmp[side]);
272 numNeutEvent[iNNeut]++;
273 }
274 }
275 }
276 }
277
278 if (m_debugLevel >= 0) {
279 std::cout << "Statistics for side " << side << ", lb range = " << LBLow << ", " << LBHigh
280 << " : # high energy events = " << numHEEvent;
281
282 for (size_t iNNeut = 0; iNNeut < numNeutMax; iNNeut++) {
283 std::cout << ", # of " + std::to_string(iNNeut) + " N events = " << numNeutEvent[iNNeut];
284 }
285
286 std::cout << std::endl;
287 }
288
289
290 // Normalize the averages
291 //
292 for (size_t iNNeut = 0; iNNeut < numNeutMax; iNNeut++) {
293 if (numNeutEvent[iNNeut] == 0) {
294 std::cout << "Found no " + std::to_string(iNNeut) + "N events, quitting" << std::endl;
295 return;
296 }
297
298 std::for_each(sumsNeutVec[iNNeut].begin(), sumsNeutVec[iNNeut].end(), [&](double& sum) {sum /= numNeutEvent[iNNeut];});
299 std::for_each(sums2DNeutVec[iNNeut].begin(), sums2DNeutVec[iNNeut].end(), [&](double& sum) {sum /= numNeutEvent[iNNeut];});
300 }
301
302 if (!excludeHE) {
303 if (numHEEvent == 0) {
304 std::cout << "Found no high-energy events, quitting" << std::endl;
305 return;
306 }
307
308 std::for_each(sumsHE.begin(), sumsHE.end(), [&](double& sum) {sum /= numHEEvent;});
309 std::for_each(sums2DHE.begin(), sums2DHE.end(), [&](double& sum) {sum /= numHEEvent;});
310 }
311
312
313 // Now we have the averages, set up to do the minimization
314 //
315 size_t numFitParams = numWeights;
316 TMatrixD minimMatrix(numFitParams, numFitParams);
317 TVectorD minimVector(numFitParams);
318
319 FillMinimizationData(minimMatrix, minimVector, maxPowerModule, HEDeweight,
320 sumsNeutVec, sumsHE, sums2DNeutVec, sums2DHE);
321
322 if (m_debugLevel >= 2) {
323 std::cout << "Dumping matrix equation: " << std::endl;
324
325 for (size_t index1 = 0; index1 < numFitParams; index1++) {
326 std::cout << "| ";
327
328 for (size_t index2 = 0; index2 < numFitParams; index2++) {
329 std::cout << std::scientific << std::setw(9) << std::setprecision(2) << minimMatrix(index1, index2) << " ";
330 }
331
332 std::cout << "| | w_" << index1 << " = | " << minimVector(index1) << std::endl;
333 }
334 std::cout << std::endl;
335 }
336
337
338 TDecompLU lu(minimMatrix);
339 lu.Decompose();
340
341 bool status = lu.Solve(minimVector);
342
343 if (!status) {
344 std::cout << "Minimization failed for side " << side << ", lb range "
345 << LBLow << ", " << LBHigh << endl;
346 }
347 else {
348 CalibData calibData;
349 calibData.LBStart = LBLow;
350 calibData.LBEnd = LBHigh;
351
352 int weightIndex = 0;
353 for (size_t module : {0, 1, 2, 3}) {
354 calibData.weights[module].assign(maxPowerModule[module], 0);
355
356 for (size_t power = 0; power < maxPowerModule[module]; power++) {
357 calibData.weights.at(module).at(power) = minimVector(weightIndex++);
358
359 if (m_debugLevel >= 0) {
360 std::cout << "Weight for module,power " << module << ", " << power << " = " << calibData.weights.at(module).at(power)
361 << std::endl;
362 }
363 }
364 }
365
366 // Now store this calibration
367 //
368 AddCalibration(side, calibOutput, calibData);
369 }
370}
void AddCalibration(size_t side, const std::string &tag, const CalibData &calib)
CalibData GetCalibration(size_t side, const std::string &tag)
void FillMinimizationData(TMatrixD &minimMatrix, TVectorD &minimVector, std::array< int, 4 > maxPowerModule, float HEDeweight, const std::vector< std::vector< double > > &sums1DVec, const std::vector< double > &sumsHE, const std::vector< std::vector< double > > &sums2DVec, const std::vector< double > &sumsHE2D)
static double CalculateEnergy(const float *moduleAmps, const CalibData &calib)
Float_t zdc_ZdcModuleAmp[2][4]
void AddToSums(std::vector< double > &sums1D, std::vector< double > &sums2D, float *amps)
status
Definition merge.py:16
unsigned int LBEnd
unsigned int LBStart
std::array< std::vector< float >, 4 > weights

◆ FillLumiBlockEvtMap()

void ZDCNLCalibration::FillLumiBlockEvtMap ( )
private

Definition at line 17 of file ZDCNLCalibration.cxx.

18{
19 Long64_t numEntries = m_tree->GetEntriesFast();
20
21 unsigned int countPassed = 0;
22
23 int lastLumiBlock = -1;
24 int LBStartEntry = -1;
25
26 for (Long64_t entry = 0; entry < numEntries; entry++) {
27 Long64_t numBytes = m_tree->GetEntry(entry);
28
29 if (m_useGRL && passBits !=0) {
30 //std::cout << "lumiBlock = " << lumiBlock << " passBits = " << passBits << std::endl;
31 continue; // reject based on GRL or errors
32 }
33
34 if (lumiBlock != lastLumiBlock) {
35 if (lastLumiBlock > -1) {
36 m_LumiBlockEvtMap.insert(LBEvtMap::value_type(lumiBlock, std::pair<unsigned int, unsigned int>(LBStartEntry, entry - 1)));
37 }
38
39 LBStartEntry = entry;
40 lastLumiBlock = lumiBlock;
41 }
42
43 countPassed++;
44 }
45
46 if (m_debugLevel >= 0) {
47 std::cout << "From a total of " << numEntries << " entries, " << countPassed << " were selected for analysis" << std::endl;
48 }
49}

◆ FillMinimizationData()

void ZDCNLCalibration::FillMinimizationData ( TMatrixD & minimMatrix,
TVectorD & minimVector,
std::array< int, 4 > maxPowerModule,
float HEDeweight,
const std::vector< std::vector< double > > & sums1DVec,
const std::vector< double > & sumsHE,
const std::vector< std::vector< double > > & sums2DVec,
const std::vector< double > & sumsHE2D )
private

Definition at line 484 of file ZDCNLCalibration.cxx.

488{
489 size_t numWeights = minimVector.GetNrows();
490 size_t numColumns = minimMatrix.GetNcols();
491 size_t numRows = minimMatrix.GetNrows();
492 size_t numNeutMax = sums1DVec.size();
493
494 if (numRows != numWeights || numColumns != numWeights){
495 throw std::runtime_error("ZDCNLCalibration::FillMinimizationData; incompatible row/col numbers");
496 }
497
498 // if (m_debugLevel > 0) std::cout << "Filling single neutorn entries in matrix" << std::endl;
499
500 double sumFracSqr = 0;
501 for (size_t module : {0, 1, 2, 3}) {
502 sumFracSqr += m_HEFraction.at(module)*m_HEFraction.at(module);
503 }
504
505 // Start with the single neutron contributions
506 //
507 {
508 size_t index1 = 0;
509
510 for (size_t module1 : {0, 1, 2, 3}) {
511 for (size_t power1 = 0; power1 < maxPowerModule[module1]; power1++) {
512
513 size_t index2 = 0;
514 size_t sumIndexOffset = (module1*m_maxNLPower + power1)*4*m_maxNLPower;
515
516 for (size_t module2 : {0, 1, 2, 3}) {
517 for (size_t power2 = 0; power2 < maxPowerModule[module2]; power2++) {
518
519 // The N neutron contributions
520 //
521 unsigned int sum2DIndex = sumIndexOffset + power2;
522 double element = 0;
523 std::for_each(sums2DVec.begin(), sums2DVec.end(), [&] (const std::vector<double>& sums) {element += 2*sums.at(sum2DIndex);});
524
525 // The TN contribution
526 //
527 // element += 2*sumsTN2D.at(sumIndexOffset + power2);
528
529 // The HE fraction contribution
530 //
531 double sum = sumsHE2D.at(sumIndexOffset + power2);
532 sum *= HEDeweight;
533
534 if (module2 == module1) element += 2*sum;
535 element -= 2*m_HEFraction.at(module2)*sum;
536 element -= 2*m_HEFraction.at(module1)*sum;
537 element += 2*sumFracSqr*sum;
538
539 minimMatrix(index1, index2) += element;
540
541 if (m_debugLevel > 2) {
542 std::cout << "Filling module 1, power 1 = " << module1 << ", " << power1 << ", module 2, power 2 = "
543 << module2 << ", " << power2 << " with sum at index " << sumIndexOffset + power2
544 << ", value = " << minimMatrix(index1, index2) << std::endl;
545 }
546
547 index2++;
548 }
549
550 sumIndexOffset += m_maxNLPower;
551 }
552
553 size_t sum1DIndex = module1*m_maxNLPower + power1;
554
555 minimVector[index1] = 0;
556 for (unsigned int iNNeut = 0; iNNeut < numNeutMax; iNNeut++) {
557 minimVector[index1] += 2*((double) iNNeut + 1)*m_SNEnergy*sums1DVec[iNNeut].at(sum1DIndex);
558 }
559
560 if (m_debugLevel > 2) {
561 std::cout << "Filling vector module , power = " << module1 << ", " << power1 << " with sum at index " << sum1DIndex
562 << ", value = " << minimVector[index1] << std::endl;
563 }
564
565 // minimVector[index1] = 2*m_SNEnergy*sumsSN.at(sum1DIndex) + 2*2*m_SNEnergy*sumsTN.at(sum1DIndex);
566 index1++;
567 }
568 }
569 }
570
571 // if (m_debugLevel > 0) std::cout << "Filling high energy entries in matrix" << std::endl;
572
573 // {
574
575 // size_t sumIndexOffset = 0;
576 // size_t index1 = 0;
577
578 // for (size_t module1 : {0, 1, 2, 3}) {
579 // for (size_t power1 = 0; power1 < maxPowerModule[module1]; power1++) {
580 // size_t index2 = 0;
581
582 // for (size_t module2 : {0, 1, 2, 3}) {
583 // for (size_t power2 = 0; power2 < maxPowerModule[module2]; power2++) {
584
585 // double element = 0;
586 // double sum = sumsHE2D.at(sumIndexOffset + power2);
587
588 // // Now assemble to terms contributing to this matrix element
589 // //
590 // if (module2 == module1) element += 2*sum;
591 // element -= 2*m_HEFraction.at(module2)*sum;
592 // element -= 2*m_HEFraction.at(module1)*sum;
593 // element += 2*sumFracSqr*sum;
594
595 // minimMatrix(index1, index2) += m_HEDeweight*element;
596 // index2++;
597 // }
598
599 // sumIndexOffset += m_maxNLPower;
600 // }
601
602 // index1++;
603 // }
604 // }
605 // }
606
607}
const std::vector< float > m_HEFraction
const float m_SNEnergy

◆ FindSNPeaks()

std::array< float, 4 > ZDCNLCalibration::FindSNPeaks ( size_t LBLow,
size_t LBHigh,
size_t side )

◆ FindSNRange()

std::pair< float, float > ZDCNLCalibration::FindSNRange ( size_t LBLow,
size_t LBHigh,
size_t side )

Definition at line 72 of file ZDCNLCalibration.cxx.

73{
74 static TH1D* snHist = new TH1D("ZDCNLCalibSNHist", "", 200, 0, 2000);
75 snHist->Reset();
76
77 LBEvtMap::const_iterator LBStart = m_LumiBlockEvtMap.lower_bound(LBLow);
78 LBEvtMap::const_iterator LBEnd = m_LumiBlockEvtMap.upper_bound(LBHigh);
79
80 // Select which trigger we're going to look at based on the side
81 //
82 const bool& trigger = (side == 1 ? L1_ZDC_C : L1_ZDC_A);
83
84 for ( LBEvtMap::const_iterator iter = LBStart; iter != LBEnd; ++iter) {
85 unsigned int entryStart = iter->second.first;
86 unsigned int entryEnd = iter->second.second;
87
88 for (unsigned int entry = entryStart; entry < entryEnd; entry++) {
89 Long64_t nb = m_tree->GetEntry(entry);
90 if (nb < 0) {
91 std::cout << "Error reading entry " << entry << ", quitting" << std::endl;
92 throw std::runtime_error("Error reading entry for ZDCNLCalibration");
93 }
94
95 if (!trigger) continue;
96
97 int testMask = zdc_ZdcModuleMask >> side*4;
98 if ((testMask & 0xf) != 0xf) continue;
99
100 float sumAmp = zdc_ZdcAmp[side];
101 snHist->Fill(sumAmp);
102 }
103 }
104
105 int ibin = snHist->GetMaximumBin();
106 double center = snHist->GetBinCenter(ibin);
107 double maxCounts = snHist->GetBinContent(ibin);
108
109 static TF1* gaussFunc = new TF1("gaussFunc", "gaus", 100, 2000);
110
111 gaussFunc->SetParameters(maxCounts, center, 100);
112 snHist->Fit("gaussFunc", "", "", center*0.7, center*1.3);
113
114 double gaussMean = gaussFunc->GetParameter(1);
115 double gaussWidth = gaussFunc->GetParameter(2);
116
117 double cutLow = gaussMean -2.*gaussWidth;
118 double cutHigh = gaussMean + 2*gaussWidth;
119
120 if (m_debugLevel >= 0) {
121 std::cout << "SN cut range = " << cutLow << " to " << cutHigh << std::endl;
122 }
123
124 return std::pair<float, float>(cutLow, cutHigh);
125}

◆ FindSNTwoNRanges()

std::pair< std::pair< float, float >, std::pair< float, float > > ZDCNLCalibration::FindSNTwoNRanges ( size_t LBLow,
size_t LBHigh,
size_t side )

Definition at line 128 of file ZDCNLCalibration.cxx.

129{
130 static TH1D* sntnHist = new TH1D("ZDCNLCalibSNTNHist", "", 200, 0, 2000);
131 sntnHist->Reset();
132
133 LBEvtMap::const_iterator LBStart = m_LumiBlockEvtMap.lower_bound(LBLow);
134 LBEvtMap::const_iterator LBEnd = m_LumiBlockEvtMap.upper_bound(LBHigh);
135
136 // Select which trigger we're going to look at based on the side
137 //
138 const bool& trigger = (side == 1 ? L1_ZDC_C : L1_ZDC_A);
139
140 for ( LBEvtMap::const_iterator iter = LBStart; iter != LBEnd; ++iter) {
141 unsigned int entryStart = iter->second.first;
142 unsigned int entryEnd = iter->second.second;
143
144 for (unsigned int entry = entryStart; entry < entryEnd; ++entry) {
145 Long64_t nb = m_tree->GetEntry(entry);
146 if (nb < 0) {
147 std::cout << "Error reading entry " << entry << ", quitting" << std::endl;
148 throw std::runtime_error("Error reading entry for ZDCNLCalibration");
149 }
150
151 if (!trigger) continue;
152
153 int testMask = zdc_ZdcModuleMask >> side*4;
154 if ((testMask & 0xf) != 0xf) continue;
155
156 float sumAmp = zdc_ZdcAmp[side];
157 sntnHist->Fill(sumAmp);
158 }
159 }
160
161 int ibin = sntnHist->GetMaximumBin();
162 double center = sntnHist->GetBinCenter(ibin);
163 double maxCounts = sntnHist->GetBinContent(ibin);
164
165 static TF1* twogaussFunc = new TF1("twogaussFunc", "[0]*exp(-pow((x-[1])/[2],2)/2.) + [3]*exp(-pow((x-[4])/[5],2)/2.)", 100, 2000);
166
167 twogaussFunc->SetParameters(maxCounts, center, 100, maxCounts/2, 2*center, 200);
168 sntnHist->Fit("twogaussFunc", "", "", center*0.7, 2*center*1.3);
169
170 double gaussMean1 = twogaussFunc->GetParameter(1);
171 double gaussMean2 = twogaussFunc->GetParameter(4);
172
173 double gaussWidth1 = std::abs(twogaussFunc->GetParameter(2));
174 double gaussWidth2 = std::abs(twogaussFunc->GetParameter(5));
175
176 double cutLow1 = gaussMean1 -2.*gaussWidth1;
177 double cutHigh1 = gaussMean1 + 2*gaussWidth1
178;
179 double cutLow2 = gaussMean2 -2.*gaussWidth2;
180 double cutHigh2 = gaussMean2 + 2*gaussWidth2;
181
182 if (cutHigh1 > cutLow2) {
183 double middle = 0.5*(cutHigh1 + cutLow2);
184 cutHigh1 = middle;
185 cutLow2 = middle;
186 }
187
188 if (m_debugLevel >= 0) {
189 std::cout << "SN cut range = " << cutLow1 << " to " << cutHigh1 << std::endl;
190 std::cout << "TN cut range = " << cutLow2 << " to " << cutHigh2 << std::endl;
191 }
192
193 return std::pair<std::pair<float, float>, std::pair<float, float> >(std::pair<float, float>(cutLow1, cutHigh1),
194 std::pair<float, float>(cutLow2, cutHigh2));
195}

◆ GetCalibration()

CalibData ZDCNLCalibration::GetCalibration ( size_t side,
const std::string & tag )
private

Definition at line 61 of file ZDCNLCalibration.cxx.

62{
63 CalibData null;
64
65 auto iter = m_calibrations.at(side).find(tag);
66 if (iter != m_calibrations.at(side).end()) {
67 return iter->second;
68 }
69 else return null;
70}

◆ GetTestFracHist()

TH1 * ZDCNLCalibration::GetTestFracHist ( size_t module)
inline

Definition at line 142 of file ZDCNLCalibration.h.

142{return (m_haveTest ? m_testCalibHEFracHist.at(module) : 0);}
std::array< TH1D *, 4 > m_testCalibHEFracHist

◆ GetTestSNHist()

TH1 * ZDCNLCalibration::GetTestSNHist ( )
inline

Definition at line 141 of file ZDCNLCalibration.h.

141{return (m_haveTest ? m_testCalibSNHist : 0);}

◆ GetTestTree()

TTree * ZDCNLCalibration::GetTestTree ( )
inline

Definition at line 144 of file ZDCNLCalibration.h.

144{return m_testTree;}

◆ SetDefaultCalibration()

void ZDCNLCalibration::SetDefaultCalibration ( size_t side,
const CalibData & calib )
inline

Definition at line 130 of file ZDCNLCalibration.h.

130 {
131 AddCalibration(side, "default", calib);
132 }

◆ TestCalibration()

void ZDCNLCalibration::TestCalibration ( int side,
const std::string & calibName )

Definition at line 372 of file ZDCNLCalibration.cxx.

373{
374 std::string calibName = (name != "" ? name : "default");
375 CalibData calib = GetCalibration(side, calibName);
376
377 if (m_debugLevel >= 0) {
378 std::cout << "Testing calibration for side " << side << "name: " << calibName << std::endl;
379 }
380
381 m_testCalibSNHist = new TH1D("testCalibSNHist" , "", 4000, 0, 200000);
382 for (int module : {0, 1, 2, 3}) {
383 std::string fracHistName = "testCalibHEFracHist_" + std::to_string(module);
384 std::string energyHistName = "testCalibEnergyHist_" + std::to_string(module);
385
386 m_testCalibHEFracHist[module] = new TH1D(fracHistName.c_str(), "", 200, 0., 1.);
387 m_testCalibEnergyHist[module] = new TH1D(energyHistName.c_str(), "", 1000, 0., 100000);
388 }
389
390 static std::array<float, 4> treeRawAmp;
391 static std::array<float, 4> treeCalibAmp;
392 static float treeRawSum, treeCalibSum, treeRawSumOpp;
393 static unsigned int treeEventNumber;
394 static bool passedSN, passedTN, passedHE;
395
396 if (m_testTree != 0) delete m_testTree;
397
398 m_testTree = new TTree("ZDCNLCalibTest", "ZDCNLCalibTest");
399 m_testTree->SetDirectory(0);
400
401 m_testTree->Branch("EventNumber", &treeEventNumber, "eventNumber/I");
402 m_testTree->Branch("lumiBlock", &lumiBlock, "lumiBlock/I");
403 m_testTree->Branch("L1_ZDC_A", &L1_ZDC_A, "L1_ZDC_A/b");
404 m_testTree->Branch("L1_ZDC_C", &L1_ZDC_C, "L1_ZDC_C/b");
405 m_testTree->Branch("L1_ZDC_A_C", &L1_ZDC_A_C, "L1_ZDC_A_C/b");
406
407 m_testTree->Branch("RawSum", &treeRawSum, "RawSum/f");
408 m_testTree->Branch("CalibSum", &treeCalibSum, "CalibSum/f");
409 m_testTree->Branch("RawAmp", &treeRawAmp[0], "RawAmp[4]/f");
410 m_testTree->Branch("CalibAmp", &treeCalibAmp[0], "CalibAmp[4]/f");
411 m_testTree->Branch("RawSumOpp", &treeRawSumOpp, "RawSumOpp/f");
412
413 LBEvtMap::const_iterator LBStart = m_LumiBlockEvtMap.lower_bound(calib.LBStart);
414 LBEvtMap::const_iterator LBEnd = m_LumiBlockEvtMap.upper_bound(calib.LBEnd);
415
416 int count = 0;
417 for ( LBEvtMap::const_iterator iter = LBStart; iter != LBEnd; ++iter) {
418 unsigned int entryStart = iter->second.first;
419 unsigned int entryEnd = iter->second.second;
420
421 for (unsigned int entry = entryStart; entry < entryEnd; ++entry) {
422 Long64_t nb = m_tree->GetEntry(entry);
423 if (nb < 0) {
424 std::cout << "Error reading entry " << entry << ", quitting" << std::endl;
425 throw std::runtime_error("Error reading entry for ZDCNLCalibration");
426 }
427
428 int testMask = zdc_ZdcModuleMask >> side*4;
429 if ((testMask & 0xf) != 0xf) continue;
430
431 double sumAmp = zdc_ZdcAmp[side];
432
433 // Calculate the calibrated energy
434 //
435 float sumCalibAmp = 0;
436 float weightIndex = 0;
437
438 std::array<double, 4> calibAmpModule = {0, 0, 0, 0};
439
440 if (m_debugLevel > 2 && count < 20) {
441 std::cout << "TestCalibration:: dumping calibration results for event " << count << std::endl;
442 }
443
444 for (size_t module : {0, 1, 2, 3}) {
445 float amp = zdc_ZdcModuleAmp[side][module];
446 float ampPow = amp;
447
448 for (size_t power = 0; power < calib.weights[module].size(); power++) {
449 calibAmpModule[module] += calib.weights[module][power]*ampPow;
450 ampPow *= amp;
451 }
452
453 double energy = calibAmpModule[module];
454 sumCalibAmp += energy;
455 m_testCalibEnergyHist[module]->Fill(energy);
456
457 treeRawAmp[module] = amp;
458 treeCalibAmp[module] = energy;
459
460 if (m_debugLevel > 2 && count < 20) {
461 std::cout << "Module " << module << " : raw, calibrated energy = " << std::fixed << amp << ", " << energy << std::endl;
462 }
463 }
464
465 if (m_debugLevel > 2 && count < 20) {
466 std::cout << "Total: raw, calibrated energy = " << std::fixed << sumAmp << ", " << sumCalibAmp << std::endl;
467 }
468
469 treeRawSum = sumAmp;
470 treeRawSumOpp = (side == 0 ? zdc_ZdcAmp[1] : zdc_ZdcAmp[0]);
471 treeCalibSum = sumCalibAmp;
472
473 m_testTree->Fill();
474
475 m_testCalibSNHist->Fill(sumCalibAmp);
476
477 count++;
478 }
479 }
480
481 m_haveTest = true;
482}
std::array< TH1D *, 4 > m_testCalibEnergyHist
int count(std::string s, const std::string &regx)
count how many occurances of a regx are in a string
Definition hcg.cxx:146

Member Data Documentation

◆ b_bcid

TBranch* ZDCNLCalibration::b_bcid
private

Definition at line 90 of file ZDCNLCalibration.h.

◆ b_eventNumber

TBranch* ZDCNLCalibration::b_eventNumber
private

Definition at line 88 of file ZDCNLCalibration.h.

◆ b_L1_ZDC_A

TBranch* ZDCNLCalibration::b_L1_ZDC_A
private

Definition at line 99 of file ZDCNLCalibration.h.

◆ b_L1_ZDC_A_C

TBranch* ZDCNLCalibration::b_L1_ZDC_A_C
private

Definition at line 102 of file ZDCNLCalibration.h.

◆ b_L1_ZDC_AND

TBranch* ZDCNLCalibration::b_L1_ZDC_AND
private

Definition at line 101 of file ZDCNLCalibration.h.

◆ b_L1_ZDC_C

TBranch* ZDCNLCalibration::b_L1_ZDC_C
private

Definition at line 100 of file ZDCNLCalibration.h.

◆ b_lumiBlock

TBranch* ZDCNLCalibration::b_lumiBlock
private

Definition at line 89 of file ZDCNLCalibration.h.

◆ b_passBits

TBranch* ZDCNLCalibration::b_passBits
private

Definition at line 92 of file ZDCNLCalibration.h.

◆ b_runNumber

TBranch* ZDCNLCalibration::b_runNumber
private

Definition at line 87 of file ZDCNLCalibration.h.

◆ b_zdc_ZdcAmp

TBranch* ZDCNLCalibration::b_zdc_ZdcAmp
private

Definition at line 94 of file ZDCNLCalibration.h.

◆ b_zdc_ZdcModuleAmp

TBranch* ZDCNLCalibration::b_zdc_ZdcModuleAmp
private

Definition at line 97 of file ZDCNLCalibration.h.

◆ b_zdc_ZdcModuleMask

TBranch* ZDCNLCalibration::b_zdc_ZdcModuleMask
private

Definition at line 96 of file ZDCNLCalibration.h.

◆ bcid

UInt_t ZDCNLCalibration::bcid
private

Definition at line 73 of file ZDCNLCalibration.h.

◆ eventNumber

UInt_t ZDCNLCalibration::eventNumber
private

Definition at line 71 of file ZDCNLCalibration.h.

◆ L1_ZDC_A

Bool_t ZDCNLCalibration::L1_ZDC_A
private

Definition at line 81 of file ZDCNLCalibration.h.

◆ L1_ZDC_A_C

Bool_t ZDCNLCalibration::L1_ZDC_A_C
private

Definition at line 84 of file ZDCNLCalibration.h.

◆ L1_ZDC_AND

Bool_t ZDCNLCalibration::L1_ZDC_AND
private

Definition at line 83 of file ZDCNLCalibration.h.

◆ L1_ZDC_C

Bool_t ZDCNLCalibration::L1_ZDC_C
private

Definition at line 82 of file ZDCNLCalibration.h.

◆ lumiBlock

UInt_t ZDCNLCalibration::lumiBlock
private

Definition at line 72 of file ZDCNLCalibration.h.

◆ m_calibrations

std::array<std::map<std::string, CalibData>, 2> ZDCNLCalibration::m_calibrations

Definition at line 112 of file ZDCNLCalibration.h.

◆ m_debugLevel

int ZDCNLCalibration::m_debugLevel
private

Definition at line 62 of file ZDCNLCalibration.h.

◆ m_haveTest

bool ZDCNLCalibration::m_haveTest {}

Definition at line 114 of file ZDCNLCalibration.h.

114{};

◆ m_HEFraction

const std::vector<float> ZDCNLCalibration::m_HEFraction
private

Definition at line 67 of file ZDCNLCalibration.h.

◆ m_LumiBlockEvtMap

LBEvtMap ZDCNLCalibration::m_LumiBlockEvtMap
private

Definition at line 108 of file ZDCNLCalibration.h.

◆ m_maxNLPower

size_t ZDCNLCalibration::m_maxNLPower
private

Definition at line 60 of file ZDCNLCalibration.h.

◆ m_SNEnergy

const float ZDCNLCalibration::m_SNEnergy
private

Definition at line 66 of file ZDCNLCalibration.h.

◆ m_testCalibEnergyHist

std::array<TH1D*, 4> ZDCNLCalibration::m_testCalibEnergyHist

Definition at line 117 of file ZDCNLCalibration.h.

◆ m_testCalibHEFracHist

std::array<TH1D*, 4> ZDCNLCalibration::m_testCalibHEFracHist

Definition at line 116 of file ZDCNLCalibration.h.

◆ m_testCalibSNHist

TH1D* ZDCNLCalibration::m_testCalibSNHist {}

Definition at line 115 of file ZDCNLCalibration.h.

115{};

◆ m_testTree

TTree* ZDCNLCalibration::m_testTree {}

Definition at line 118 of file ZDCNLCalibration.h.

118{};

◆ m_TFile

TFile* ZDCNLCalibration::m_TFile
private

Definition at line 57 of file ZDCNLCalibration.h.

◆ m_tree

TTree* ZDCNLCalibration::m_tree
private

Definition at line 58 of file ZDCNLCalibration.h.

◆ m_useGRL

bool ZDCNLCalibration::m_useGRL
private

Definition at line 61 of file ZDCNLCalibration.h.

◆ passBits

UInt_t ZDCNLCalibration::passBits
private

Definition at line 74 of file ZDCNLCalibration.h.

◆ runNumber

UInt_t ZDCNLCalibration::runNumber
private

Definition at line 70 of file ZDCNLCalibration.h.

◆ zdc_ZdcAmp

Float_t ZDCNLCalibration::zdc_ZdcAmp[2]
private

Definition at line 76 of file ZDCNLCalibration.h.

◆ zdc_ZdcModuleAmp

Float_t ZDCNLCalibration::zdc_ZdcModuleAmp[2][4]
private

Definition at line 79 of file ZDCNLCalibration.h.

◆ zdc_ZdcModuleMask

UInt_t ZDCNLCalibration::zdc_ZdcModuleMask
private

Definition at line 78 of file ZDCNLCalibration.h.


The documentation for this class was generated from the following files: