Skip to content

Commit 75216ea

Browse files
bquilainbquilain
committed
Bug found by S. Izumiyama in the time residual PDF production which makes a slight bump in the -100ns region. I re-wrote ProducePDF.cpp to correctly rebin the histogram and evaluate the spline. ProducePDF and PDF in inputs has been modified
1 parent 5fa1df9 commit 75216ea

File tree

2 files changed

+59
-80
lines changed

2 files changed

+59
-80
lines changed

‎inputs/timePDF_DRnew_Large.root

-2.43 KB
Binary file not shown.

‎macros/ProducePDF.cpp

+59-80
Original file line numberDiff line numberDiff line change
@@ -1,55 +1,9 @@
1-
#include <iostream>
2-
#include <sstream>
3-
#include <fstream>
4-
#include <iomanip>
5-
#include <vector>
6-
#include <map>
7-
#include <TROOT.h>
8-
#include <TApplication.h>
9-
#include <TStyle.h>
10-
#include <TFile.h>
11-
#include <TTree.h>
12-
#include <TCanvas.h>
13-
#include <TChain.h>
14-
#include <TBranch.h>
15-
#include <TF1.h>
16-
#include <TH1.h>
17-
#include <TH2.h>
18-
#include <TH3.h>
19-
#include <TMath.h>
20-
#include <TGraph.h>
21-
#include <TGraphErrors.h>
22-
#include <TMinuit.h>
23-
#include <TFitter.h>
24-
#include <TLegend.h>
25-
#include <TGraph2D.h>
26-
#include <TSpline.h>
27-
#include <getopt.h>
28-
29-
using namespace std;
30-
31-
int main(int argc, char **argv){
1+
void ProducePDF(){
322
TString ofname = "timePDF.root";
333
TString ifname = "Analyze_hkhybridmpmt10pc14374100Hz_4.2kHzbl_10MeV.root"; // Originally used in B.Quilain codes
344

355
int opt;
366

37-
while ((opt = getopt(argc, argv, "f:o:h")) != -1){
38-
switch (opt) {
39-
case 'f':
40-
ifname = TString(optarg);
41-
break;
42-
case 'o':
43-
ofname = TString(optarg);
44-
break;
45-
default:
46-
cerr << "Usage: " << argv[0] << " [-f input] [-o output]" << endl;
47-
break;
48-
}
49-
}
50-
51-
52-
537
const int nFiles=1;
548
TFile * _f[nFiles];
559
const int nPMTtypes = 2;
@@ -715,51 +669,73 @@ int main(int argc, char **argv){
715669
fOut->cd();
716670
splineExpoConv[f][i]->Write(splineExpoConv[f][i]->GetTitle());
717671

718-
const int nBins=270;//150;//Number of big bins
719-
double xPos[nBins];
720-
double yPos[nBins];
721-
double drPos[nBins];
672+
cout<<"Start to merge bins in higher size bins to optimize PDF"<<endl;
673+
vector <double> xPos;xPos.clear();
674+
vector <double> yPos;yPos.clear();
675+
vector <double> drPos;drPos.clear();
722676
int iBinActive=0;
723677
int nRebins=10;//Number of small bin gathered in a big one
724678
double xAverage=0;
725679
double yAverage=0;
726680
double drAverage=0;
727681
int binLimit=HitTimeTOFProfile[f][i]->FindBin(-limitFitGausExpo[i]);//HitTimeTOFProfile[f][i]->FindBin(limitFitGausExpo[i]-5);//Lower limit of the active big bin
728-
729-
682+
const int nLimits = 8;
683+
double minValue = -limitFitGausExpo[i];
684+
double maxValue = HitTimeTOFProfile[f][i]->GetXaxis()->GetXmax();
685+
double lowLimits[nLimits] = {minValue, -100., -30., -10., 8., 30., 100.,maxValue};
686+
double binLowLimits[nLimits];
687+
for(int il = 0;il<nLimits;il++){
688+
binLowLimits[il] = HitTimeTOFProfile[f][i]->FindBin(lowLimits[il]);
689+
}
690+
int rebinLowLimits[nLimits] = {100, 10, 5, 2, 10, 30, 100, 100};
691+
int currentLimit = 0;
692+
int nBinsAverage = 0;
693+
int nBins=0;//150;//Number of big bins
694+
//int iBinActive = 0;
695+
//In between the different limits, the rebinning factor will be different.
696+
//We should first ensure in which region we are. Based on this, we use a different rebinning value.
697+
//Then, we will loop over the bins which are in the correct region. We will have an internal counter. As soon as this internal counter reach the rebinning value, we stop and average.
698+
//Another way of reaching the factor could be to reach another regiob. We should take this into account.
699+
730700
//for(int ibinx=HitTimeTOFProfile[f][i]->FindBin(limitFitGausExpo[i]-5);ibinx<=HitTimeTOFProfile[f][i]->GetNbinsX();ibinx++){
701+
//Loop over all bins of the histogram
731702
for(int ibinx=HitTimeTOFProfile[f][i]->FindBin(-limitFitGausExpo[i]);ibinx<=HitTimeTOFProfile[f][i]->GetNbinsX();ibinx++){
732-
xAverage+=HitTimeTOFProfile[f][i]->GetBinCenter(ibinx);
733-
yAverage+=HitTimeTOFProfile[f][i]->GetBinContent(ibinx);
734-
drAverage+=HitTimeTOFDR[f][i]->GetBinContent(ibinx);
735-
736-
//if(HitTimeTOFProfile[f][i]->GetBinCenter(ibinx)<-3) nRebins=1;
737-
if(HitTimeTOFProfile[f][i]->GetBinCenter(ibinx)<-100) nRebins=100;
738-
else if(HitTimeTOFProfile[f][i]->GetBinCenter(ibinx)<-30) nRebins=10;
739-
else if(HitTimeTOFProfile[f][i]->GetBinCenter(ibinx)<-10) nRebins=5;
740-
else if(HitTimeTOFProfile[f][i]->GetBinCenter(ibinx)<8) nRebins=2;
741-
else if(HitTimeTOFProfile[f][i]->GetBinCenter(ibinx)<30) nRebins=10;
742-
else if(HitTimeTOFProfile[f][i]->GetBinCenter(ibinx)<100) nRebins=30;
743-
else nRebins =100;
744-
//cout<<"Rebins="<<nRebins<<", bin limit="<<ibinx-binLimit<<", average x="<<xAverage<<endl;
745-
if((ibinx-binLimit)%nRebins==(nRebins-1)){
746-
xPos[iBinActive]=xAverage/nRebins;
747-
drPos[iBinActive]=drAverage/nRebins;
748-
yPos[iBinActive]=(yAverage/nRebins);
749-
cout<<"Position = "<<xPos[iBinActive]<<", value="<<yPos[iBinActive]<<", DR = "<<drPos[iBinActive]<<endl;
750-
iBinActive++;
703+
cout<<"Bin = "<<ibinx<<" i.e. value of low edge = "<< HitTimeTOFProfile[f][i]->GetBinLowEdge(ibinx)<<", current low edge limit bin = "<<binLowLimits[currentLimit]<<", nBinsAverage = "<<nBinsAverage<<", rebinning factor = "<<rebinLowLimits[currentLimit]<<endl;
704+
if((ibinx >= binLowLimits[currentLimit] && ibinx < binLowLimits[currentLimit+1]) && nBinsAverage < rebinLowLimits[currentLimit]){
705+
//Sum over bin content to average them until we reach the last bin
706+
xAverage+=HitTimeTOFProfile[f][i]->GetBinCenter(ibinx);
707+
yAverage+=HitTimeTOFProfile[f][i]->GetBinContent(ibinx);
708+
drAverage+=HitTimeTOFDR[f][i]->GetBinContent(ibinx);
709+
nBinsAverage++;
710+
}
711+
if( (ibinx+1 >= binLowLimits[currentLimit+1]) || (nBinsAverage >= rebinLowLimits[currentLimit]) ){
712+
//If we reached the number of binning to rebin: store information.
713+
//Same if we reached the limit of the region to rebin of a given factor.
714+
//What if both happens at the same time, or worse: we reach limit of bins on bin n, and at n+1, we overcome the limit. In that case, we do not have anyting to fill our average? So, we understand that if we are in the last bin below the limit, we should store and then pass to the next step of limit. So, the check of the limit should always be on the next bin.
715+
xPos.push_back(xAverage/nBinsAverage);
716+
drPos.push_back(drAverage/nBinsAverage);
717+
yPos.push_back((yAverage/nBinsAverage));
718+
cout<<"Position = "<<xPos[nBins]<<", value="<<yPos[nBins]<<", DR = "<<drPos[nBins]<<endl;
719+
nBins++;
720+
nBinsAverage = 0;
751721
xAverage=0;
752722
yAverage=0;
753723
drAverage=0;
754-
binLimit=ibinx+1;
755-
}
756-
if(iBinActive>=nBins){
757-
cout<<"Add more bins"<<endl;
758-
break;
724+
if(ibinx+1 >= binLowLimits[currentLimit+1]){
725+
currentLimit++;
726+
}
759727
}
760728
}
761-
762-
graphExpoQueue[f][i] = new TGraph(nBins,xPos,yPos);
729+
double nBins_graph = nBins;
730+
double * xPos_graph = new double[nBins];
731+
double * yPos_graph = new double[nBins];
732+
double * drPos_graph = new double[nBins];
733+
for(int ibin = 0; ibin<nBins;ibin++){
734+
xPos_graph[ibin] = xPos.at(ibin);
735+
yPos_graph[ibin] = yPos.at(ibin);
736+
drPos_graph[ibin] = drPos.at(ibin);
737+
}
738+
graphExpoQueue[f][i] = new TGraph(nBins_graph,xPos_graph,yPos_graph);
763739

764740
//HitTimeTOFProfile[f][i]->Rebin(4);
765741
//ExpoQueue[f][i] = new TF1(Form("ExpoQueue%d_%d",f,i),"[0]*TMath::Exp(-[1]*x)",limitFitGausExpo[i],500);
@@ -773,12 +749,15 @@ int main(int argc, char **argv){
773749
//splineExpoQueue[f][i] = new TSpline3(Form("splineExpoQueue%d_%d",f,i),limitFitGausExpo[i],500.,ExpoQueue[f][i],300);
774750
splineExpoQueue[f][i] = new TSpline3(Form("splineExpoQueue%d_%d",f,i),graphExpoQueue[f][i]);
775751
splineExpoQueue[f][i]->SetLineColor(kCyan);
752+
graphExpoQueue[f][i]->SetMarkerColor(kCyan);
753+
graphExpoQueue[f][i]->SetMarkerStyle(20);
754+
graphExpoQueue[f][i]->Draw("Psame");
776755
splineExpoQueue[f][i]->Draw("lcsame");
777756
fOut->cd();
778757
splineExpoQueue[f][i]->Write(splineExpoQueue[f][i]->GetTitle());
779758

780759

781-
graphDR[f][i] = new TGraph(nBins,xPos,drPos);
760+
graphDR[f][i] = new TGraph(nBins_graph,xPos_graph,drPos_graph);
782761
//splineDR[f][i] = new TSpline3(Form("splineDR%d_%d",f,i),limitFitGausExpo[i],500.,DR[f][i],300);
783762
splineDR[f][i] = new TSpline3(Form("splineDR%d_%d",f,i),graphDR[f][i]);
784763
splineDR[f][i]->SetLineColor(kCyan);

0 commit comments

Comments
 (0)