227 B2INFO(
"Starting eclNOptimalAlgorithm");
233 const int lastBin = inputParameters->GetNbinsX();
234 const double scale = inputParameters->GetBinContent(lastBin);
235 for (
int ib = 1; ib < lastBin; ib++) {
236 double param = inputParameters->GetBinContent(ib);
237 inputParameters->SetBinContent(ib, param / scale);
238 inputParameters->SetBinError(ib, 0.);
244 const double groupNum = groupNumberOfEachCellID->GetBinContent(ic);
245 groupNumberOfEachCellID->SetBinContent(ic, groupNum / scale);
246 groupNumberOfEachCellID->SetBinError(ic, 0.);
254 TFile* histFile =
new TFile(
"eclNOptimalAlgorithm.root",
"recreate");
255 inputParameters->Write();
256 groupNumberOfEachCellID->Write();
257 entriesPerThetaIdEnergy->Write();
258 mcEnergyDiff->Write();
262 const int nCrystalGroups = (int)(inputParameters->GetBinContent(1) + 0.5);
263 const int nEnergies = (int)(inputParameters->GetBinContent(2) + 0.5);
264 const int nLeakReg = 3;
265 auto generatedE = create2Dvector<float>(nLeakReg, nEnergies);
267 for (
int ireg = 0; ireg < nLeakReg; ireg++) {
268 for (
int ie = 0; ie < nEnergies; ie++) {
270 generatedE[ireg][ie] = inputParameters->GetBinContent(bin);
273 B2INFO(
"nCrystalGroups = " << nCrystalGroups);
278 const int iExp = exprun[0].first;
279 const int iRun = exprun[0].second;
280 B2INFO(
"Experiment, run = " << iExp <<
", " << iRun);
298 B2INFO(
"eclNOptimalAlgorithm: reading previous payload");
300 std::vector<float> eUpperBoundariesFwdPrev = existingECLnOptimal->getUpperBoundariesFwd();
301 std::vector<float> eUpperBoundariesBrlPrev = existingECLnOptimal->getUpperBoundariesBrl();
302 std::vector<float> eUpperBoundariesBwdPrev = existingECLnOptimal->getUpperBoundariesBwd();
303 TH2F nOptimal2DPrev = existingECLnOptimal->getNOptimal();
304 nOptimal2DPrev.SetName(
"nOptimal2DPrev");
306 const int nPrevE = eUpperBoundariesFwdPrev.size();
307 B2INFO(
"Energy upper boundaries from previous payload for each region");
308 for (
int ie = 0; ie < nPrevE; ie++) {
309 B2INFO(
" " << ie <<
" " << eUpperBoundariesFwdPrev[ie] <<
" " << eUpperBoundariesBrlPrev[ie] <<
" " << eUpperBoundariesBwdPrev[ie]);
314 nOptimal2DPrev.Write();
318 std::vector< std::vector<int> > initialNOptimal;
319 initialNOptimal.resize(nCrystalGroups, std::vector<int>(nEnergies, 0));
322 TH2F* nInitialPerGroup =
new TH2F(
"nInitialPerGroup",
"initial nCrys, energy point vs group number;group;energy point",
323 nCrystalGroups, 0., nCrystalGroups, nEnergies, 0., nEnergies);
324 TH2F* generatedEPerGroup =
new TH2F(
"generatedEPerGroup",
"Generated energy, energy point vs group number;group;energy point",
325 nCrystalGroups, 0., nCrystalGroups, nEnergies, 0., nEnergies);
326 TH1F* firstCellIDPerGroup =
new TH1F(
"firstCellIDPerGroup",
"First cellID of group;group", nCrystalGroups, 0., nCrystalGroups);
327 TH1F* nCrystalsPerGroup =
new TH1F(
"nCrystalsPerGroup",
"Number of crystals per group;group", nCrystalGroups, 0., nCrystalGroups);
328 std::vector<int> crystalsPerGroup;
329 crystalsPerGroup.resize(nCrystalGroups, 0);
332 for (
int cellID = 1; cellID <= nCrystalsTotal; cellID++) {
335 int iGroup = (int)(0.5 + groupNumberOfEachCellID->GetBinContent(cellID));
338 std::vector<float> eUpperBoundariesPrev = eUpperBoundariesBrlPrev;
341 eUpperBoundariesPrev = eUpperBoundariesFwdPrev;
344 eUpperBoundariesPrev = eUpperBoundariesBwdPrev;
351 for (
int ie = 0; ie < nEnergies; ie++) {
352 float energy = generatedE[iRegion][ie];
354 while (energy > eUpperBoundariesPrev[iEnergy]) {iEnergy++;}
355 initialNOptimal[iGroup][ie] = (int)(0.5 + nOptimal2DPrev.GetBinContent(iGroup + 1, iEnergy + 1));
357 B2INFO(
"ie iEnergy mismatch: cellID " << cellID <<
" ie " << ie <<
" iEnergy " << iEnergy);
361 generatedEPerGroup->SetBinContent(iGroup + 1, ie + 1, energy);
362 nInitialPerGroup->SetBinContent(iGroup + 1, ie + 1, initialNOptimal[iGroup][ie]);
366 if (crystalsPerGroup[iGroup] == 0) {
367 firstCellIDPerGroup->SetBinContent(iGroup + 1, cellID);
368 firstCellIDPerGroup->SetBinError(iGroup + 1, 0.);
370 crystalsPerGroup[iGroup]++;
374 for (
int ig = 0; ig < nCrystalGroups; ig++) {
375 nCrystalsPerGroup->SetBinContent(ig + 1, crystalsPerGroup[ig]);
376 nCrystalsPerGroup->SetBinError(ig + 1, 0.);
379 nInitialPerGroup->Write();
380 generatedEPerGroup->Write();
381 firstCellIDPerGroup->Write();
382 nCrystalsPerGroup->Write();
383 B2INFO(
"Wrote initial diagnostic histograms");
387 for (
int ig = 0; ig < nCrystalGroups; ig++) {
388 for (
int ie = 0; ie < nEnergies; ie++) {
389 std::string name =
"eSum_" + std::to_string(ig) +
"_" + std::to_string(ie);
391 name =
"biasSum_" + std::to_string(ig) +
"_" + std::to_string(ie);
398 B2INFO(
"Wrote eSum and biasSum histograms");
405 TH2F* nOptimalPerGroup =
new TH2F(
"nOptimalPerGroup",
"nOptimal;group;energy point", nCrystalGroups, 0., nCrystalGroups, nEnergies,
407 TH2F* changeInNOptimal =
new TH2F(
"changeInNOptimal",
"nOptimal minus nInitial;group;energy point", nCrystalGroups, 0.,
408 nCrystalGroups, nEnergies, 0., nEnergies);
409 TH2F* binsInFit =
new TH2F(
"binsInFit",
"Number of bins used in Novo fit, energy point vs group number;group;energy point",
410 nCrystalGroups, 0., nCrystalGroups, nEnergies, 0., nEnergies);
411 TH2F* maxEntriesPerHist =
new TH2F(
"maxEntriesPerHist",
412 "Max entries in histogram bin, energy point vs group number;group;energy point", nCrystalGroups, 0., nCrystalGroups, nEnergies, 0.,
414 TH2F* peakPerGroup =
new TH2F(
"peakPerGroup",
"peak energy, energy point vs group number;group;energy point", nCrystalGroups, 0.,
415 nCrystalGroups, nEnergies, 0., nEnergies);
416 TH2F* effSigmaPerGroup =
new TH2F(
"effSigmaPerGroup",
"fit effective sigma, energy point vs group number;group;energy point",
417 nCrystalGroups, 0., nCrystalGroups, nEnergies, 0., nEnergies);
418 TH2F* etaPerGroup =
new TH2F(
"etaPerGroup",
"fit eta, energy point vs group number;group;energy point", nCrystalGroups, 0.,
419 nCrystalGroups, nEnergies, 0., nEnergies);
420 TH2F* fitProbPerGroup =
new TH2F(
"fitProbPerGroup",
"fit probability, energy point vs group number;group;energy point",
421 nCrystalGroups, 0., nCrystalGroups, nEnergies, 0., nEnergies);
422 TH2F* resolutionPerGroup =
new TH2F(
"resolutionPerGroup",
"resolution/(peak-bias), energy point vs group number;group;energy point",
423 nCrystalGroups, 0., nCrystalGroups, nEnergies, 0., nEnergies);
424 TH2F* fracContainedEPerGroup =
new TH2F(
"fracContainedEPerGroup",
425 "peak fraction of energy contained in nOpt crystals;group;energy point", nCrystalGroups, 0., nCrystalGroups, nEnergies, 0.,
428 TH2F* biasPerGroup =
new TH2F(
"biasPerGroup",
"bias (GeV), energy point vs group number;group;energy point", nCrystalGroups, 0.,
429 nCrystalGroups, nEnergies, 0., nEnergies);
430 TH2F* sigmaBiasPerGroup =
new TH2F(
"sigmaBiasPerGroup",
"sigma bias (GeV), energy point vs group number;group;energy point",
431 nCrystalGroups, 0., nCrystalGroups, nEnergies, 0., nEnergies);
433 TH2F* existingResolutionPerGroup =
new TH2F(
"existingResolutionPerGroup",
434 "existing resolution/peak, energy point vs group number;group;energy point", nCrystalGroups, 0., nCrystalGroups, nEnergies, 0.,
441 TF1* func =
new TF1(
"eclNOptimalNovo", eclNOptimalNovo, 0.4, 1.4, 4);
442 func->SetLineColor(kRed);
444 for (
int ig = 0; ig < nCrystalGroups; ig++) {
445 for (
int ie = 0; ie < nEnergies; ie++) {
448 std::string name =
"eSum_" + std::to_string(ig) +
"_" + std::to_string(ie);
452 std::string biasName =
"biasSum_" + std::to_string(ig) +
"_" + std::to_string(ie);
456 const double genEnergy = generatedEPerGroup->GetBinContent(ig + 1, ie + 1);
461 const int nCrysBins = eSum->GetNbinsX();
462 const int nCrysMax = nCrysBins - 1;
465 TH1D* eSumOpt{
nullptr};
467 double bestFractionalResolution = 999.;
468 double existingFractionalResolution = 999.;
469 double biasForNOpt = 0.;
470 double biasSigmaForNOpt = 0.;
476 const int initialnCrysSumToFit = initialNOptimal[ig][ie];
477 int nCrysSumToFit = initialnCrysSumToFit;
478 while (nCrysSumToFit > 0) {
480 TH1D* hEnergy = (TH1D*)eSum->ProjectionY(
"hEnergy", nCrysSumToFit, nCrysSumToFit);
481 TString newName = name +
"_" + std::to_string(nCrysSumToFit);
482 hEnergy->SetName(newName);
486 const double minESumEntries = 800.;
487 if (hEnergy->GetEntries() < minESumEntries) {
488 B2INFO(
"Insufficient entries in eSum: " << hEnergy->GetEntries() <<
" for group " << ig <<
" energy point " << ie);
490 B2INFO(
"closed histFile; quitting");
495 const double minPeakValue = 300.;
496 while (hEnergy->GetMaximum() < minPeakValue) {hEnergy->Rebin(2);}
500 double fitFraction = 0.5;
502 int iHi = hEnergy->GetNbinsX();
503 const int maxBinsToFit = 59;
504 while ((iHi - iLo) > maxBinsToFit) {
505 std::vector<int> iBins = eclNOptimalFitRange(hEnergy, fitFraction);
508 if ((iHi - iLo) > maxBinsToFit) {hEnergy->Rebin(2);}
512 const int nBinsX = hEnergy->GetNbinsX();
513 const double fitlow = hEnergy->GetBinLowEdge(iLo);
514 const double fithigh = hEnergy->GetBinLowEdge(iHi + 1);
515 double sum = hEnergy->Integral(iLo, iHi);
516 B2INFO(
"group " << ig <<
" ie " << ie <<
" name " << name <<
" newName " << newName <<
" nBinsX " << nBinsX <<
" 50% target: " <<
517 fitFraction * hEnergy->GetEntries() <<
" iLo: " << iLo <<
" iHi: " << iHi <<
" sum: " << sum);
521 double normalization = hEnergy->GetMaximum();
522 const int iMax = hEnergy->GetMaximumBin();
523 double peak = hEnergy->GetBinLowEdge(iMax);
524 double effSigma = 0.5 * (fithigh - fitlow);
527 func->SetParameters(normalization, peak, effSigma, eta);
528 func->SetParNames(
"normalization",
"peak",
"effSigma",
"eta");
531 func->SetParLimits(1, fitlow, fithigh);
532 func->SetParLimits(2, 0., 2.*effSigma);
533 func->SetParLimits(3, -1, 3);
536 if (nCrysSumToFit == 1) {
537 const double etaForOneCrystal = 0.95;
538 func->SetParameter(3, etaForOneCrystal);
539 func->SetParLimits(3, etaForOneCrystal, etaForOneCrystal);
543 hEnergy->Fit(func,
"LIEQ",
"", fitlow, fithigh);
549 double biasSigma = 0.;
550 if (nCrysSumToFit < nCrysBins) {
551 TH1D* hBias = (TH1D*)biasSum->ProjectionY(
"hBias", nCrysSumToFit, nCrysSumToFit);
553 std::vector<int> jBins = eclNOptimalFitRange(hBias, fitFraction);
554 const double lowEdge = hBias->GetBinLowEdge(jBins[0]);
555 const double highEdge = hBias->GetBinLowEdge(jBins[1] + 1);
556 bias = 0.5 * (lowEdge + highEdge);
557 biasSigma = 0.5 * (highEdge - lowEdge);
564 std::vector<int> iBins = eclNOptimalFitRange(hEnergy, fitFraction);
565 const double resolution = eclNOptimalResolution(hEnergy, iBins[0], iBins[1]);
570 peak = func->GetParameter(1);
571 const double fractionalResolution = resolution / (peak - bias / genEnergy);
575 if (fractionalResolution < bestFractionalResolution and nCrysSumToFit != nCrysBins) {
576 bestFractionalResolution = fractionalResolution;
577 nOpt = nCrysSumToFit;
579 nFitBins = 1 + iHi - iLo;
581 biasSigmaForNOpt = biasSigma;
582 }
else if (nCrysSumToFit == nCrysBins) {
583 existingFractionalResolution = fractionalResolution;
590 const double resTolerance = 1.05;
591 if (nCrysSumToFit == initialnCrysSumToFit) {
594 if (nCrysSumToFit > 1) {
599 }
else if (nCrysSumToFit < initialnCrysSumToFit) {
603 if (fractionalResolution < resTolerance* bestFractionalResolution and nCrysSumToFit > 1) {
605 }
else if (initialnCrysSumToFit != nCrysMax) {
606 nCrysSumToFit = initialnCrysSumToFit + 1;
608 nCrysSumToFit = nCrysBins;
610 }
else if (nCrysSumToFit < nCrysBins) {
614 if (fractionalResolution < resTolerance * bestFractionalResolution and nCrysSumToFit < nCrysMax) {
617 nCrysSumToFit = nCrysBins;
619 }
else if (nCrysSumToFit == nCrysBins) {
626 B2INFO(
" ig: " << ig <<
" ie: " << ie <<
" initial nOpt: " << initialnCrysSumToFit <<
" final nOpt: " << nOpt);
632 TF1* funcOpt = (TF1*)eSumOpt->GetFunction(
"eclNOptimalNovo");
634 nOptimalPerGroup->SetBinContent(ig + 1, ie + 1, nOpt);
636 changeInNOptimal->SetBinContent(ig + 1, ie + 1, nOpt - initialnCrysSumToFit);
637 changeInNOptimal->SetBinError(ig + 1, ie + 1, 0.);
639 biasPerGroup->SetBinContent(ig + 1, ie + 1, biasForNOpt);
640 biasPerGroup->SetBinError(ig + 1, ie + 1, 0.);
642 sigmaBiasPerGroup->SetBinContent(ig + 1, ie + 1, biasSigmaForNOpt);
643 sigmaBiasPerGroup->SetBinError(ig + 1, ie + 1, 0.);
645 binsInFit->SetBinContent(ig + 1, ie + 1, nFitBins);
646 binsInFit->SetBinError(ig + 1, ie + 1, 0.);
648 maxEntriesPerHist->SetBinContent(ig + 1, ie + 1, eSumOpt->GetMaximum());
649 maxEntriesPerHist->SetBinError(ig + 1, ie + 1, 0.);
651 const double peakOpt = funcOpt->GetParameter(1);
652 const double peakOptUnc = funcOpt->GetParError(1);
653 peakPerGroup->SetBinContent(ig + 1, ie + 1, peakOpt);
654 peakPerGroup->SetBinError(ig + 1, ie + 1, peakOptUnc);
656 const double genE = generatedEPerGroup->GetBinContent(ig + 1, ie + 1);
657 const double peakCor = peakOpt - biasForNOpt / genE;
658 fracContainedEPerGroup->SetBinContent(ig + 1, ie + 1, peakCor);
659 fracContainedEPerGroup->SetBinError(ig + 1, ie + 1, peakOptUnc);
661 effSigmaPerGroup->SetBinContent(ig + 1, ie + 1, funcOpt->GetParameter(2));
662 effSigmaPerGroup->SetBinError(ig + 1, ie + 1, funcOpt->GetParError(2));
664 etaPerGroup->SetBinContent(ig + 1, ie + 1, funcOpt->GetParameter(3));
665 etaPerGroup->SetBinError(ig + 1, ie + 1, funcOpt->GetParError(3));
667 fitProbPerGroup->SetBinContent(ig + 1, ie + 1, funcOpt->GetProb());
668 fitProbPerGroup->SetBinError(ig + 1, ie + 1, 0.);
670 resolutionPerGroup->SetBinContent(ig + 1, ie + 1, bestFractionalResolution);
671 resolutionPerGroup->SetBinError(ig + 1, ie + 1, 0.);
673 existingResolutionPerGroup->SetBinContent(ig + 1, ie + 1, existingFractionalResolution);
674 existingResolutionPerGroup->SetBinError(ig + 1, ie + 1, 0.);
686 changeInNOptimal->Write();
687 biasPerGroup->Write();
688 sigmaBiasPerGroup->Write();
690 maxEntriesPerHist->Write();
691 peakPerGroup->Write();
692 effSigmaPerGroup->Write();
693 etaPerGroup->Write();
694 fitProbPerGroup->Write();
695 resolutionPerGroup->Write();
696 existingResolutionPerGroup->Write();
697 fracContainedEPerGroup->Write();
699 B2INFO(
"Wrote fit summary histograms");
707 TH2F* fracContainedAdjacent[2];
708 TH2F* biasAdjacent[2];
709 const TString updown[2] = {
"lower",
"higher"};
710 for (
int ia = 0; ia < 2; ia++) {
711 TString name =
"fracContainedAdjacent_" + std::to_string(ia);
712 TString title =
"peak fraction of energy contained in nOpt crystals, " + updown[ia] +
" E;group;energy point";
713 fracContainedAdjacent[ia] =
new TH2F(name, title, nCrystalGroups, 0., nCrystalGroups, nEnergies, 0., nEnergies);
715 name =
"biasAdjacent_" + std::to_string(ia);
716 title =
"bias (GeV) in nOpt crystals, " + updown[ia] +
" E;group;energy point";
717 biasAdjacent[ia] =
new TH2F(name, title, nCrystalGroups, 0., nCrystalGroups, nEnergies, 0., nEnergies);
722 for (
int ig = 0; ig < nCrystalGroups; ig++) {
723 for (
int ie = 0; ie < nEnergies; ie++) {
726 const int nOpt = nOptimalPerGroup->GetBinContent(ig + 1, ie + 1);
730 for (
int ia = 0; ia < 2; ia++) {
731 const int ieAdj = ie + 2 * ia - 1;
732 double eFracPeakAdj = fracContainedEPerGroup->GetBinContent(ig + 1, ie + 1);
733 double biasAdj = biasPerGroup->GetBinContent(ig + 1, ie + 1);
736 if (ieAdj >= 0 and ieAdj < nEnergies) {
740 std::string biasName =
"biasSum_" + std::to_string(ig) +
"_" + std::to_string(ieAdj);
742 TH1D* hBias = (TH1D*)biasSum->ProjectionY(
"hBias", nOpt, nOpt);
745 double fitFraction = 0.683;
746 std::vector<int> jBins = eclNOptimalFitRange(hBias, fitFraction);
747 const double lowEdge = hBias->GetBinLowEdge(jBins[0]);
748 const double highEdge = hBias->GetBinLowEdge(jBins[1] + 1);
749 biasAdj = 0.5 * (lowEdge + highEdge);
753 std::string name =
"eSum_" + std::to_string(ig) +
"_" + std::to_string(ieAdj);
755 TH1D* hEnergy = (TH1D*)eSum->ProjectionY(
"hEnergy", nOpt, nOpt);
756 TString newName = name +
"_" + std::to_string(nOpt);
757 hEnergy->SetName(newName);
759 B2INFO(
"fit adj ig = " << ig <<
" ie = " << ie <<
" ia = " << ia <<
" " << newName <<
" " <<
" entries = " << hEnergy->GetEntries()
760 <<
" GetMaxiumum = " << hEnergy->GetMaximum());
763 const double minPeakValue = 300.;
764 while (hEnergy->GetMaximum() < minPeakValue) {hEnergy->Rebin(2);}
769 int iHi = hEnergy->GetNbinsX();
770 const int maxBinsToFit = 59;
771 while ((iHi - iLo) > maxBinsToFit) {
772 std::vector<int> iBins = eclNOptimalFitRange(hEnergy, fitFraction);
775 if ((iHi - iLo) > maxBinsToFit) {hEnergy->Rebin(2);}
779 const double fitlow = hEnergy->GetBinLowEdge(iLo);
780 const double fithigh = hEnergy->GetBinLowEdge(iHi + 1);
781 double normalization = hEnergy->GetMaximum();
782 const int iMax = hEnergy->GetMaximumBin();
783 double peak = hEnergy->GetBinLowEdge(iMax);
784 double effSigma = 0.5 * (fithigh - fitlow);
787 func->SetParameters(normalization, peak, effSigma, eta);
788 func->SetParNames(
"normalization",
"peak",
"effSigma",
"eta");
791 func->SetParLimits(1, fitlow, fithigh);
792 func->SetParLimits(2, 0., 2.*effSigma);
793 func->SetParLimits(3, -1, 3);
797 const double etaForOneCrystal = 0.95;
798 func->SetParameter(3, etaForOneCrystal);
799 func->SetParLimits(3, etaForOneCrystal, etaForOneCrystal);
803 hEnergy->Fit(func,
"LIEQ",
"", fitlow, fithigh);
806 const double peakAdj = func->GetParameter(1);
807 const double genE = generatedEPerGroup->GetBinContent(ig + 1, ieAdj + 1);
808 eFracPeakAdj = peakAdj - biasAdj / genE;
812 fracContainedAdjacent[ia]->SetBinContent(ig + 1, ie + 1, eFracPeakAdj);
813 biasAdjacent[ia]->SetBinContent(ig + 1, ie + 1, biasAdj);
816 B2INFO(
" ig " << ig <<
" ie " << ie <<
" eFrac (3): "
817 << fracContainedAdjacent[0]->GetBinContent(ig + 1, ie + 1) <<
" "
818 << fracContainedEPerGroup->GetBinContent(ig + 1, ie + 1) <<
" "
819 << fracContainedAdjacent[1]->GetBinContent(ig + 1, ie + 1)
821 << biasAdjacent[0]->GetBinContent(ig + 1, ie + 1) <<
" "
822 << biasPerGroup->GetBinContent(ig + 1, ie + 1) <<
" "
823 << biasAdjacent[1]->GetBinContent(ig + 1, ie + 1)
831 fracContainedAdjacent[0]->Write();
832 fracContainedAdjacent[1]->Write();
833 biasAdjacent[0]->Write();
834 biasAdjacent[1]->Write();
836 B2INFO(
"Wrote adjacent energy point summary histograms");
843 auto boundaryE = create2Dvector<float>(nLeakReg, nEnergies);
844 for (
int ireg = 0; ireg < nLeakReg; ireg++) {
845 B2INFO(
"Generated energies and boundaries for region = " << ireg);
846 for (
int ie = 0; ie < nEnergies - 1; ie++) {
847 double logE0 = log(generatedE[ireg][ie]);
848 double logE1 = log(generatedE[ireg][ie + 1]);
849 double logBoundary = 0.5 * (logE0 + logE1);
850 boundaryE[ireg][ie] = exp(logBoundary);
851 B2INFO(
" " << ie <<
" " << generatedE[ireg][ie] <<
" " << boundaryE[ireg][ie] <<
" GeV");
855 boundaryE[ireg][nEnergies - 1] = 9999.;
856 B2INFO(
" " << nEnergies - 1 <<
" " << generatedE[ireg][nEnergies - 1] <<
" " << boundaryE[ireg][nEnergies - 1] <<
" GeV");
860 std::vector<int> groupNumber;
862 const int iGroup = (int)(0.5 + groupNumberOfEachCellID->GetBinContent(cellID));
863 groupNumber.push_back(iGroup);
867 TH2F nOptimal2D(
"nOptimal2D",
"Optimal nCrys, energy bin vs group number", nCrystalGroups, 0, nCrystalGroups, nEnergies, 0.,
869 for (
int ig = 0; ig < nCrystalGroups; ig++) {
870 for (
int ie = 0; ie < nEnergies; ie++) {
871 double nOpt = nOptimalPerGroup->GetBinContent(ig + 1, ie + 1);
872 nOptimal2D.SetBinContent(ig + 1, ie + 1, nOpt);
878 TH2F peakFracEnergy(
"peakFracEnergy",
"peak contained energy over generated E, energy bin vs group number", 3 * nCrystalGroups, 0,
879 nCrystalGroups, nEnergies, 0., nEnergies);
880 TH2F bias(
"bias",
"median bias (GeV), energy bin vs group number", 3 * nCrystalGroups, 0, nCrystalGroups, nEnergies, 0., nEnergies);
881 TH2F logPeakEnergy(
"logPeakEnergy",
"log peak Energy (GeV), energy bin vs group number", 3 * nCrystalGroups, 0, nCrystalGroups,
882 nEnergies, 0., nEnergies);
883 for (
int ig = 0; ig < nCrystalGroups; ig++) {
884 for (
int ie = 0; ie < nEnergies; ie++) {
885 const int iy = ie + 1;
886 const int ix = 1 + 3 * ig;
891 double peakAdj = fracContainedEPerGroup->GetBinContent(ig + 1, ie + 1);
892 const double genE = generatedEPerGroup->GetBinContent(ig + 1, ie + 1);
893 double logE = log(genE * peakAdj);
895 peakFracEnergy.SetBinContent(ix, iy, peakAdj);
896 logPeakEnergy.SetBinContent(ix, iy, logE);
899 double peakAdj0 = peakAdj;
902 peakAdj0 = fracContainedAdjacent[0]->GetBinContent(ig + 1, ie + 1);
903 const double genE0 = generatedEPerGroup->GetBinContent(ig + 1, ie);
904 logE0 = log(genE0 * peakAdj0);
907 peakFracEnergy.SetBinContent(ix + 1, iy, peakAdj0);
908 logPeakEnergy.SetBinContent(ix + 1, iy, logE0);
911 double peakAdj1 = peakAdj;
913 if (ie < nEnergies - 1) {
914 peakAdj1 = fracContainedAdjacent[1]->GetBinContent(ig + 1, ie + 1);
915 const double genE1 = generatedEPerGroup->GetBinContent(ig + 1, ie + 2);
916 logE1 = log(genE1 * peakAdj1);
919 peakFracEnergy.SetBinContent(ix + 2, iy, peakAdj1);
920 logPeakEnergy.SetBinContent(ix + 2, iy, logE1);
923 double medianBias = biasPerGroup->GetBinContent(ig + 1, ie + 1);
924 bias.SetBinContent(ix, iy, medianBias);
925 medianBias = biasAdjacent[0]->GetBinContent(ig + 1, ie + 1);
926 bias.SetBinContent(ix + 1, iy, medianBias);
927 medianBias = biasAdjacent[1]->GetBinContent(ig + 1, ie + 1);
928 bias.SetBinContent(ix + 2, iy, medianBias);
935 peakFracEnergy.Write();
936 logPeakEnergy.Write();
945 for (
int ireg = 0; ireg < nLeakReg; ireg++) {
946 std::vector<float> eUpperBoundaries;
947 for (
int ie = 0; ie < nEnergies; ie++) {eUpperBoundaries.push_back(boundaryE[ireg][ie]);}
950 }
else if (ireg == 1) {
970 B2RESULT(
"eclNOptimalAlgorithm: successfully stored payload ECLnOptimal");