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
211
212 size_t numWeights = std::accumulate(maxPowerModule.begin(), maxPowerModule.end(), 0);
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
227
228
229
232
233 std::vector<int> numNeutEvent(numNeutMax, 0);
234 int numHEEvent = 0;
235
236
237
238 std::string calibName = (calibInput != "" ? calibInput : "default");
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
253 if ((testMask & 0xf) != 0xf) continue;
254
256
257 if (triggerSame && sumAmp > heSumThresh && !excludeHE) {
259 numHEEvent++;
260 }
261
262
263
264 if (triggerOpp && sumAmp > nNeutERange[0].first && sumAmp < nNeutERange[0].second) {
266 numNeutEvent[0]++;
267 }
268
269 for (size_t iNNeut = 1; iNNeut < numNeutMax; iNNeut++) {
270 if (sumAmp > nNeutERange[iNNeut].first && sumAmp < nNeutERange[iNNeut].second) {
272 numNeutEvent[iNNeut]++;
273 }
274 }
275 }
276 }
277
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
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
314
315 size_t numFitParams = numWeights;
316 TMatrixD minimMatrix(numFitParams, numFitParams);
317 TVectorD minimVector(numFitParams);
318
320 sumsNeutVec, sumsHE, sums2DNeutVec, sums2DHE);
321
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;
350 calibData.
LBEnd = LBHigh;
351
352 int weightIndex = 0;
353 for (size_t module : {0, 1, 2, 3}) {
355
356 for (
size_t power = 0; power < maxPowerModule[
module]; power++) {
357 calibData.
weights.at(module).at(power) = minimVector(weightIndex++);
358
360 std::cout <<
"Weight for module,power " <<
module <<
", " << power <<
" = " << calibData.
weights.at(module).at(power)
361 << std::endl;
362 }
363 }
364 }
365
366
367
369 }
370}
void AddCalibration(size_t side, const std::string &tag, const CalibData &calib)
CalibData GetCalibration(size_t side, const std::string &tag)
LBEvtMap m_LumiBlockEvtMap
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)
std::array< std::vector< float >, 4 > weights