1- void runCalibration ( TString lInputFileName = "results/AR_544122_glauberNBD_ancestorMode2_hFT0C_BCs.root" , double anchorPointPercentage = 90.0 , double matchRange = 200.0 , bool doNpartNcoll = false){
2- TFile * file = new TFile (lInputFileName .Data (), "READ" );
1+ void runCalibration (TString lInputFileName = "results/AR_544122_glauberNBD_ancestorMode2_hFT0C_BCs.root" , double anchorPointPercentage = 90.0 , double matchRange = 200.0 , bool doNpartNcoll = false)
2+ {
3+ TFile * file = new TFile (lInputFileName .Data (), "READ" );
34 file -> ls ();
4-
5- TH1F * hData = (TH1F * ) file -> Get ("hV0MUltraFine" );
6- TH1F * hGlauberParameters = (TH1F * ) file -> Get ("hGlauberParameters" );
7- TH1F * hGlauberFitRange = (TH1F * ) file -> Get ("hGlauberFitRange" );
5+
6+ TH1F * hData = (TH1F * )file -> Get ("hV0MUltraFine" );
7+ TH1F * hGlauberParameters = (TH1F * )file -> Get ("hGlauberParameters" );
8+ TH1F * hGlauberFitRange = (TH1F * )file -> Get ("hGlauberFitRange" );
89 hData -> SetName ("hData" );
9- TH1F * hStitched = (TH1F * ) hData -> Clone ("hStitched" );
10- TH1F * hFit = (TH1F * ) file -> Get ("hGlauber" );
11-
12- TCanvas * c1 = new TCanvas ("c1" , "" , 800 , 600 );
10+ TH1F * hStitched = (TH1F * )hData -> Clone ("hStitched" );
11+ TH1F * hFit = (TH1F * )file -> Get ("hGlauber" );
12+
13+ TCanvas * c1 = new TCanvas ("c1" , "" , 800 , 600 );
1314 c1 -> SetLeftMargin (0.17 );
1415 c1 -> SetBottomMargin (0.17 );
1516 c1 -> SetRightMargin (0.15 );
1617 c1 -> SetTopMargin (0.05 );
17- c1 -> SetTicks (1 ,1 );
18+ c1 -> SetTicks (1 , 1 );
1819 c1 -> SetLogz ();
1920 c1 -> SetFrameFillStyle (0 );
2021 c1 -> SetFillStyle (0 );
21-
22-
23- cout <<"Data bin width: " <<hData -> GetBinWidth (1 )<<endl ;
24- cout <<"Fit bin width: " <<hFit -> GetBinWidth (1 )<<endl ;
25- cout <<"Match range to use: " <<matchRange <<endl ;
26-
22+
23+ cout << "Data bin width: " << hData -> GetBinWidth (1 ) << endl ;
24+ cout << "Fit bin width: " << hFit -> GetBinWidth (1 ) << endl ;
25+ cout << "Match range to use: " << matchRange << endl ;
26+
2727 //____________________________________________
28- double anchorPointFraction = anchorPointPercentage / 100.f ;
28+ double anchorPointFraction = anchorPointPercentage / 100.f ;
2929 double anchorPoint = -1 ; // the anchor point value in raw
30-
30+
3131 //____________________________________________
3232 // doing partial integration up to certain point for finding anchor point bin
33- for (int ii = 1 ; ii < hData -> GetNbinsX ()+ 1 ; ii ++ ){
33+ for (int ii = 1 ; ii < hData -> GetNbinsX () + 1 ; ii ++ ) {
3434 // renormalize data curve
35- int bin1 = ii + 1 ;
36- int bin2 = hData -> FindBin ( hData -> GetBinLowEdge (ii + 1 ) + matchRange + 1e-3 );
37- double matchRangeData = hData -> Integral ( bin1 , bin2 );
38- double matchRangeFit = hFit -> Integral ( bin1 , bin2 );
39-
35+ int bin1 = ii + 1 ;
36+ int bin2 = hData -> FindBin (hData -> GetBinLowEdge (ii + 1 ) + matchRange + 1e-3 );
37+ double matchRangeData = hData -> Integral (bin1 , bin2 );
38+ double matchRangeFit = hFit -> Integral (bin1 , bin2 );
39+
4040 // rescale fit to match in the vicinity of the region we're at
41- hFit -> Scale (matchRangeData /matchRangeFit );
42-
43- double integralFit = hFit -> Integral (1 ,ii );
44- double integralData = hData -> Integral (ii + 1 ,hData -> GetNbinsX ()+ 1 );
45- double integralAll = integralFit + integralData ;
46-
47- cout <<"at bin #" <<ii <<", integrated up to " <<hData -> GetBinLowEdge (ii + 1 )<<" fraction above this value is: " <<integralData /integralAll <<endl ;
48- anchorPoint = hData -> GetBinLowEdge (ii + 1 );
49-
50- if (integralData /integralAll < anchorPointFraction ) break ;
41+ hFit -> Scale (matchRangeData / matchRangeFit );
42+
43+ double integralFit = hFit -> Integral (1 , ii );
44+ double integralData = hData -> Integral (ii + 1 , hData -> GetNbinsX () + 1 );
45+ double integralAll = integralFit + integralData ;
46+
47+ cout << "at bin #" << ii << ", integrated up to " << hData -> GetBinLowEdge (ii + 1 ) << " fraction above this value is: " << integralData / integralAll << endl ;
48+ anchorPoint = hData -> GetBinLowEdge (ii + 1 );
49+
50+ if (integralData / integralAll < anchorPointFraction )
51+ break ;
5152 }
52-
53+
5354 //____________________________________________
54- for (int ii = 1 ; ii < hData -> GetNbinsX ()+ 1 ; ii ++ ){
55+ for (int ii = 1 ; ii < hData -> GetNbinsX () + 1 ; ii ++ ) {
5556 // renormalize data curve
56- if ( hData -> GetBinCenter (ii ) < anchorPoint ) hStitched -> SetBinContent (ii , hFit -> GetBinContent (ii ));
57+ if (hData -> GetBinCenter (ii ) < anchorPoint )
58+ hStitched -> SetBinContent (ii , hFit -> GetBinContent (ii ));
5759 }
58-
59-
60- cout <<"Anchor point determined to be: " <<anchorPoint <<endl ;
61- cout <<"Preparing stitched histogram ... " <<endl ;
62-
60+
61+ cout << "Anchor point determined to be: " << anchorPoint << endl ;
62+ cout << "Preparing stitched histogram ... " << endl ;
63+
6364 hFit -> SetLineColor (kRed );
6465 hStitched -> SetLineColor (kBlue );
65-
66+
6667 hData -> GetYaxis ()-> SetTitleSize (0.055 );
6768 hData -> GetXaxis ()-> SetTitleSize (0.055 );
6869 hData -> GetYaxis ()-> SetLabelSize (0.04 );
@@ -71,75 +72,74 @@ void runCalibration( TString lInputFileName = "results/AR_544122_glauberNBD_ance
7172 hData -> Draw ("hist" );
7273 hFit -> Draw ("hist same" );
7374 hStitched -> Draw ("hist same" );
74-
75- //All fine, let's try the calibrator
76- multCalibrator * lCalib = new multCalibrator ("lCalib" );
75+
76+ // All fine, let's try the calibrator
77+ multCalibrator * lCalib = new multCalibrator ("lCalib" );
7778 lCalib -> SetAnchorPointPercentage (100.0f );
7879 lCalib -> SetAnchorPointRaw (-1e-6 );
79-
80- //Set standard Pb-Pb boundaries
80+
81+ // Set standard Pb-Pb boundaries
8182 lCalib -> SetStandardOnePercentBoundaries ();
82-
83+
8384 TString calibFileName = lInputFileName .Data ();
8485 calibFileName .ReplaceAll ("glauberNBD" , "calibration" );
85- TFile * fileCalib = new TFile (calibFileName .Data (), "RECREATE" );
86-
87- TH1F * hCalib = lCalib -> GetCalibrationHistogram (hStitched , "hCalib" );
88-
89- TCanvas * c2 = new TCanvas ("c2" , "" , 800 , 600 );
86+ TFile * fileCalib = new TFile (calibFileName .Data (), "RECREATE" );
87+
88+ TH1F * hCalib = lCalib -> GetCalibrationHistogram (hStitched , "hCalib" );
89+
90+ TCanvas * c2 = new TCanvas ("c2" , "" , 800 , 600 );
9091 c2 -> SetLeftMargin (0.17 );
9192 c2 -> SetBottomMargin (0.17 );
9293 c2 -> SetRightMargin (0.15 );
9394 c2 -> SetTopMargin (0.05 );
94- c2 -> SetTicks (1 ,1 );
95+ c2 -> SetTicks (1 , 1 );
9596 // c2->SetLogz();
9697 c2 -> SetFrameFillStyle (0 );
9798 c2 -> SetFillStyle (0 );
98-
99+
99100 hCalib -> GetYaxis ()-> SetTitleSize (0.055 );
100101 hCalib -> GetXaxis ()-> SetTitleSize (0.055 );
101102 hCalib -> GetYaxis ()-> SetLabelSize (0.04 );
102103 hCalib -> GetXaxis ()-> SetLabelSize (0.04 );
103104 hCalib -> SetTitle ("" );
104105 hCalib -> Draw ();
105-
106-
106+
107107 fileCalib -> cd ();
108108
109109 hData -> Write ();
110110 hCalib -> Write ();
111111 hStitched -> Write ();
112112 hFit -> Write ();
113-
114- if (doNpartNcoll ){
115- cout << "Will now attempt to calculate % -> Np, Nc map..." << endl ;
116-
117- TProfile * hProfileNpart = new TProfile ("hProfileNpart" , "" , 100 , 0 , 100 );
118- TProfile * hProfileNcoll = new TProfile ("hProfileNcoll" , "" , 100 , 0 , 100 );
119- TH2F * h2dNpart = new TH2F ("h2dNpart" , "" , 100 , 0 , 100 , 500 , -0.5f , 499.5f );
120- TH2F * h2dNcoll = new TH2F ("h2dNcoll" , "" , 100 , 0 , 100 , 3000 , -0.5f , 2999.5 );
121-
113+
114+ if (doNpartNcoll ) {
115+ cout << "Will now attempt to calculate % -> Np, Nc map..." << endl ;
116+
117+ TProfile * hProfileNpart = new TProfile ("hProfileNpart" , "" , 100 , 0 , 100 );
118+ TProfile * hProfileNcoll = new TProfile ("hProfileNcoll" , "" , 100 , 0 , 100 );
119+ TH2F * h2dNpart = new TH2F ("h2dNpart" , "" , 100 , 0 , 100 , 500 , -0.5f , 499.5f );
120+ TH2F * h2dNcoll = new TH2F ("h2dNcoll" , "" , 100 , 0 , 100 , 3000 , -0.5f , 2999.5 );
121+
122122 // Replay
123- multGlauberNBDFitter * g = new multGlauberNBDFitter ("lglau" );
124- TF1 * fitfunc = g -> GetGlauberNBD ();
125-
126- //Step 1: open the (Npart, Ncoll) pair information, provide
127- TFile * fbasefile = new TFile ("basehistos.root" ,"READ" );
128- TH2D * hNpNc = (TH2D * ) fbasefile -> Get ("hNpNc" );
123+ multGlauberNBDFitter * g = new multGlauberNBDFitter ("lglau" );
124+ TF1 * fitfunc = g -> GetGlauberNBD ();
125+
126+ // Step 1: open the (Npart, Ncoll) pair information, provide
127+ TFile * fbasefile = new TFile ("basehistos.root" , "READ" );
128+ TH2D * hNpNc = (TH2D * )fbasefile -> Get ("hNpNc" );
129129 g -> SetNpartNcollCorrelation (hNpNc );
130130 g -> InitializeNpNc ();
131-
131+
132132 fitfunc -> SetParameter (0 , hGlauberParameters -> GetBinContent (1 ));
133133 fitfunc -> SetParameter (1 , hGlauberParameters -> GetBinContent (2 ));
134134 fitfunc -> SetParameter (2 , hGlauberParameters -> GetBinContent (3 ));
135135 fitfunc -> SetParameter (3 , hGlauberParameters -> GetBinContent (4 ));
136136 fitfunc -> SetParameter (4 , hGlauberParameters -> GetBinContent (5 ));
137-
138- Double_t lMax = hData -> GetBinLowEdge ( hData -> GetNbinsX () + 1 );
139-
137+
138+ Double_t lMax = hData -> GetBinLowEdge (hData -> GetNbinsX () + 1 );
139+
140140 // uncomment if Np Nc needed -> Warning, slow!
141- g -> CalculateAvNpNc ( hProfileNpart , hProfileNcoll , h2dNpart , h2dNcoll , hCalib , 0 , lMax );
142-
141+ g -> CalculateAvNpNc (hProfileNpart , hProfileNcoll , h2dNpart , h2dNcoll , hCalib , 0 , lMax );
142+
143143 hProfileNpart -> Write ();
144144 hProfileNcoll -> Write ();
145145 h2dNpart -> Write ();
0 commit comments