]>
Commit | Line | Data |
---|---|---|
47747ea2 | 1 | /* $Id$ */ |
2 | ||
3 | #include "TFile.h" | |
6aa36480 | 4 | #include "TCanvas.h" |
47747ea2 | 5 | #include "TDirectory.h" |
6 | #include "TNamed.h" | |
7 | #include "TF1.h" | |
8 | #include "TAttLine.h" | |
9 | #include "TAttMarker.h" | |
20fb7b55 | 10 | #include "TH1.h" |
47747ea2 | 11 | #include "TH1F.h" |
12 | #include "TH2F.h" | |
13 | #include "TSystem.h" | |
14 | #include "TStyle.h" | |
47747ea2 | 15 | #include "TMath.h" |
16 | #include "TList.h" | |
17 | #include <TString.h> | |
18 | #include <TObject.h> | |
19 | #include <TROOT.h> | |
20 | #include <THashList.h> | |
21 | #include <Rtypes.h> | |
22 | #include <TPRegexp.h> | |
23 | #include "TFitResult.h" | |
d1f44d99 | 24 | #include <TMap.h> |
25 | #include <TObjString.h> | |
6aa36480 | 26 | #include "TPad.h" |
20fb7b55 | 27 | #include "TAxis.h" |
47747ea2 | 28 | |
29 | namespace RawProduction { | |
30 | ||
a9e9e511 | 31 | bool ignoreErrors = false; |
04f85d20 | 32 | const int centBinVersion = 2; |
47747ea2 | 33 | |
34 | // Fit range | |
35 | Double_t rangeMin=0.05 ; | |
36 | Double_t rangeMax=0.3 ; | |
a9e9e511 | 37 | |
47747ea2 | 38 | // Fit limits |
39 | double lowerMass = 0.110; double upperMass = 0.145; | |
40 | double lowerWidth = 0.001; double upperWidth = 0.012; | |
41 | ||
42 | Double_t PeakPosition(Double_t pt){ | |
43 | //Fit to the measured peak position | |
44 | return 4.99292e-003*exp(-pt*9.32300e-001)+1.34944e-001 ; | |
45 | } | |
46 | //----------------------------------------------------------------------------- | |
47 | Double_t PeakWidth(Double_t pt){ | |
48 | //fit to the measured peak width | |
49 | Double_t a=0.0068, b=0.0025, c=0.000319 ; | |
50 | return TMath::Sqrt(a*a+b*b/pt/pt+c*c*pt*pt) ; | |
51 | } | |
52 | ||
53 | // Pt bin parameters | |
54 | Int_t nPtBins=0; | |
55 | Double_t ptBinEdges[1000] = {0}; | |
56 | int fool; | |
57 | double GetPtBin(int bin); | |
58 | void MakePtBins(); | |
a9e9e511 | 59 | |
60 | // various parameters | |
47747ea2 | 61 | const Double_t kMean=0.136 ; //Approximate peak position to facilitate error estimate |
62 | const char format[] = ".pdf"; // say if you want .pdf | |
a9e9e511 | 63 | |
47747ea2 | 64 | class Input; |
65 | ||
66 | TH1* GetHistogram_cent(Input& input, const TString& name, int centrality); | |
67 | TH1* MergeHistogram_cent(Input& input, const TString& name, int newCentIndex, int fromCentIndex, int toCentIndex); | |
68 | int kNCentralityBins = 3; | |
69 | ||
70 | Double_t CGausPol1(Double_t * x, Double_t * par); | |
71 | Double_t CGausPol2(Double_t * x, Double_t * par); | |
72 | Double_t CGausPol0(Double_t * x, Double_t * par); | |
73 | Double_t CPol1(Double_t * x, Double_t * par); | |
74 | Double_t CPol2(Double_t * x, Double_t * par); | |
a9e9e511 | 75 | |
76 | ||
d1f44d99 | 77 | // Output Bin |
a53ea5d3 | 78 | class TriggerBin { |
d1f44d99 | 79 | public: |
a53ea5d3 | 80 | TriggerBin(const TString& trigger = "kCentral" ); |
81 | TriggerBin(const char* trigger); | |
3374aa5b | 82 | const TString& Dir() const {return fDir;} |
d1f44d99 | 83 | const TString& Key() const {return fKey;} |
84 | const TString& Trigger() const {return fTrigger;} | |
a9e9e511 | 85 | protected: |
d1f44d99 | 86 | TString fTrigger; |
3374aa5b | 87 | TString fDir; |
d1f44d99 | 88 | TString fKey; |
89 | }; | |
47747ea2 | 90 | // Object describing the input for the macros |
91 | class Input { | |
92 | public: | |
a53ea5d3 | 93 | Input(const TString& fileName, const TriggerBin& inputBin, TString listPath = ""); |
47747ea2 | 94 | TH1 * GetHistogram(const char* name=""); |
a53ea5d3 | 95 | const TriggerBin& Bin() const {return fBin;} |
47747ea2 | 96 | private: |
97 | static TFile* fFile; | |
98 | TList* fList; | |
a53ea5d3 | 99 | TriggerBin fBin; |
47747ea2 | 100 | }; |
a9e9e511 | 101 | |
d1f44d99 | 102 | // Output Bin |
a53ea5d3 | 103 | class TriCenPidBin : public TriggerBin { |
d1f44d99 | 104 | public: |
a53ea5d3 | 105 | TriCenPidBin(Int_t centrality, const TString& pid, const TString& trigger); |
d1f44d99 | 106 | Int_t Centrality() const {return fCentrality;} |
107 | const TString& PID() const {return fPID;} | |
108 | private: | |
109 | Int_t fCentrality; | |
110 | TString fPID; | |
111 | }; | |
112 | // Object describing the output of the macros | |
113 | class Output { | |
114 | public: | |
6b5d1b53 | 115 | Output(const TString& fileName = "RawProduction.root", const char* options = "UPDATE"); |
a53ea5d3 | 116 | TH1* GetHistogram(const TString& name, const TriggerBin& inBin); |
baacba22 | 117 | TH1* GetHistogram(const TString& name); |
a53ea5d3 | 118 | void SetDir(const TriggerBin& inBin); |
d1f44d99 | 119 | void Write(); |
120 | private: | |
d1f44d99 | 121 | TFile* fFile; |
d1f44d99 | 122 | }; |
47747ea2 | 123 | |
5e8d81ec | 124 | void MakePi0Fit(Input& input, const TriggerBin& outBin, Output& output); |
125 | void MakePi0FitTCP(Input& input, const TriCenPidBin& outBin, Output& output); | |
126 | ||
127 | void MakePi0Fit(Input& input, const TriggerBin& outBin, Output& output) { | |
128 | MakePtBins(); | |
129 | Printf("\nMakePi0Fit(%s)", outBin.Key().Data()); | |
130 | output.SetDir(outBin); | |
131 | ||
132 | for(int m1i = 1; m1i <=3; ++m1i) { | |
133 | for(int m2i = m1i; m2i <=3; ++m2i) { | |
134 | TStringToken names("A;M;W;p0;p1;p2", ";"); | |
135 | TStringToken titles("Amplitude;Mass;Width;p0;p1;p2", ";"); | |
136 | while( names.NextToken() && titles.NextToken() ) { | |
137 | TString name(Form("%s_M%i%i", names.Data(), m1i, m2i)); | |
138 | TString title(Form("%s_M%i%i", titles.Data(), m1i, m2i)); | |
139 | new TH1D(name.Data(), title.Data(), nPtBins,ptBinEdges); | |
140 | new TH1D(Form("%s_error", name.Data()), title.Data(), nPtBins,ptBinEdges); | |
141 | } | |
142 | } | |
143 | } | |
144 | ||
145 | TF1 * pol2Gaus = new TF1("pol2Gaus", CGausPol2, 0., 1., 6); | |
146 | pol2Gaus->SetParNames("A", "m_{0}", "#sigma", "a_{0}", "a_{1}", "a_{2}"); | |
147 | for(int m1i = 1; m1i <=3; ++m1i) { | |
148 | for(int m2i = m1i; m2i <=3; ++m2i) { | |
149 | TH2F* hPi0MNN = (TH2F*) input.GetHistogram(Form("hPi0M%i%i", m1i, m2i)); | |
150 | for(Int_t ptBin=1; ptBin<=nPtBins; ptBin++){ | |
151 | Int_t imin = hPi0MNN->GetYaxis()->FindBin(ptBinEdges[ptBin-1]+0.0001); | |
152 | Int_t imax = hPi0MNN->GetYaxis()->FindBin(ptBinEdges[ptBin]+0.0001) -1; | |
153 | TH1D* hPi0MNNProj = hPi0MNN->ProjectionX(Form("pt%03i_hPi0M%i%i",ptBin, m1i, m2i), imin, imax); | |
154 | hPi0MNNProj->SetTitle(Form("M_{#gamma#gamma}, M%i%i, %.1f<p_{T}<%.1f, %s", m1i, m2i, ptBinEdges[ptBin-1], ptBinEdges[ptBin], outBin.Trigger().Data())); | |
155 | hPi0MNNProj->Sumw2(); | |
156 | int entries = hPi0MNNProj->Integral(hPi0MNNProj->FindBin(lowerMass), hPi0MNNProj->FindBin(upperMass)); | |
157 | Printf("pt bin %i, %3.1f to %3.1f, entries: %i", ptBin, ptBinEdges[ptBin-1], ptBinEdges[ptBin], entries); | |
158 | if( entries < 10 ) continue; | |
159 | ||
160 | ||
161 | // Error Fix | |
162 | for(Int_t ib=1; ib<=hPi0MNNProj->GetNbinsX();ib++){if(hPi0MNNProj ->GetBinContent(ib)==0)hPi0MNNProj ->SetBinError(ib,1.);} | |
163 | ||
164 | pol2Gaus->SetParameters(entries, 0.134, 0.006, entries, 0, 0); | |
165 | pol2Gaus->SetParLimits(1, lowerMass, upperMass); | |
166 | pol2Gaus->SetParLimits(2, lowerWidth, upperWidth); | |
167 | TFitResultPtr resPtr = hPi0MNNProj->Fit(pol2Gaus, "MSQ", "", rangeMin, rangeMax); | |
168 | Int_t error = int(resPtr) != 4000; | |
169 | error = error || TMath::Abs(pol2Gaus->GetParameter(1) - lowerMass) < 0.0001 || TMath::Abs(pol2Gaus->GetParameter(1) - upperMass) <0.0001; | |
170 | error = error || TMath::Abs(pol2Gaus->GetParameter(2) - lowerWidth) < 0.0001 || TMath::Abs(pol2Gaus->GetParameter(2) - upperWidth) <0.0001; | |
171 | ||
172 | if(error) { | |
173 | ((TH1D*) output.GetHistogram(Form("A_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(0)); | |
174 | ((TH1D*) output.GetHistogram(Form("M_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(1)); | |
175 | ((TH1D*) output.GetHistogram(Form("W_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(2)); | |
176 | ((TH1D*) output.GetHistogram(Form("p0_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(3)); | |
177 | ((TH1D*) output.GetHistogram(Form("p1_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(4)); | |
178 | ((TH1D*) output.GetHistogram(Form("p2_M%i%i_error", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(5)); | |
179 | ((TH1D*) output.GetHistogram(Form("A_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(0)); | |
180 | ((TH1D*) output.GetHistogram(Form("M_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(1)); | |
181 | ((TH1D*) output.GetHistogram(Form("W_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(2)); | |
182 | ((TH1D*) output.GetHistogram(Form("p0_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(3)); | |
183 | ((TH1D*) output.GetHistogram(Form("p1_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(4)); | |
184 | ((TH1D*) output.GetHistogram(Form("p2_M%i%i_error", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(5)); | |
185 | } else { | |
186 | ((TH1D*) output.GetHistogram(Form("A_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(0)); | |
187 | ((TH1D*) output.GetHistogram(Form("M_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(1)); | |
188 | ((TH1D*) output.GetHistogram(Form("W_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(2)); | |
189 | ((TH1D*) output.GetHistogram(Form("p0_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(3)); | |
190 | ((TH1D*) output.GetHistogram(Form("p1_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(4)); | |
191 | ((TH1D*) output.GetHistogram(Form("p2_M%i%i", m1i, m2i), outBin))->SetBinContent(ptBin, pol2Gaus->GetParameter(5)); | |
192 | ((TH1D*) output.GetHistogram(Form("A_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(0)); | |
193 | ((TH1D*) output.GetHistogram(Form("M_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(1)); | |
194 | ((TH1D*) output.GetHistogram(Form("W_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(2)); | |
195 | ((TH1D*) output.GetHistogram(Form("p0_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(3)); | |
196 | ((TH1D*) output.GetHistogram(Form("p1_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(4)); | |
197 | ((TH1D*) output.GetHistogram(Form("p2_M%i%i", m1i, m2i), outBin))->SetBinError(ptBin, pol2Gaus->GetParError(5)); | |
198 | } | |
199 | } | |
200 | } | |
201 | } | |
202 | } | |
203 | void MakePi0FitTCP(Input& input, const TriCenPidBin& outBin, Output& output) { | |
47747ea2 | 204 | MakePtBins(); |
20fb7b55 | 205 | Printf("\nMakePi0Fit(%s)", outBin.Key().Data()); |
6aa36480 | 206 | output.SetDir(outBin); |
a9e9e511 | 207 | |
47747ea2 | 208 | TH1F * hTotSelEvents = (TH1F*) input.GetHistogram("hTotSelEvents"); |
209 | TH2F * hCentrality = (TH2F*) input.GetHistogram("hCenPHOSCells"); | |
210 | TH1D * hCentralityX = hCentrality->ProjectionX(); | |
20fb7b55 | 211 | TH2F *hPi0 = (TH2F*)GetHistogram_cent(input, Form("hPi0%s", outBin.PID().Data()), outBin.Centrality()); |
212 | TH2F *hPi0Mix = (TH2F*)GetHistogram_cent(input, Form("hMiPi0%s", outBin.PID().Data()), outBin.Centrality()); | |
47747ea2 | 213 | |
214 | printf("TotSelEvents (4): %.0f \n", hTotSelEvents->GetBinContent(4)) ; | |
215 | printf("Centrality: %.0f \n", hCentralityX->Integral()) ; | |
a9e9e511 | 216 | |
47747ea2 | 217 | if( !hPi0 || !hPi0Mix ) { |
218 | Printf(Form("no histogram(0x%p, 0x%p), returning", hPi0, hPi0Mix)); | |
219 | return; | |
220 | } | |
221 | ||
222 | // for temp convas for drawing/monitoring | |
6aa36480 | 223 | TCanvas* canvas = new TCanvas("cMakePi0Fit", Form("MakePi0Fit Canvas, %s", outBin.Key().Data()),10,10,1200,800); |
a9e9e511 | 224 | |
225 | ||
47747ea2 | 226 | // Peak Parameterization |
227 | // 1. Linear Bg | |
228 | TF1 * funcRatioFit1 = new TF1("funcRatioFit1",CGausPol1,0.,1.,5) ; | |
229 | funcRatioFit1->SetParNames("A", "m_{0}", "#sigma", "a_{0}", "a_{1}"); | |
230 | funcRatioFit1->SetLineWidth(2) ; | |
231 | funcRatioFit1->SetLineColor(2) ; | |
232 | // 2. Quadratic Bg | |
233 | TF1 * funcRatioFit2 = new TF1("funcRatioFit2",CGausPol2,0.,1.,6) ; | |
234 | funcRatioFit2->SetParNames("A", "m_{0}", "#sigma", "a_{0}", "a_{1}", "a_{2}"); | |
235 | funcRatioFit2->SetLineWidth(2) ; | |
236 | funcRatioFit2->SetLineColor(4) ; | |
237 | funcRatioFit2->SetLineStyle(2) ; | |
238 | // Other | |
239 | TF1 * fgs = new TF1("gs",CGausPol0,0.,1.,4) ; | |
240 | fgs->SetParNames("A", "m_{0}", "#sigma", "B") ; | |
241 | fgs->SetLineColor(2) ; | |
242 | fgs->SetLineWidth(1) ; | |
243 | TF1 * fbg1 = new TF1("bg1",CPol1,0.,1.,2) ; | |
244 | TF1 * fbg2 = new TF1("bg2",CPol2,0.,1.,3) ; | |
245 | ||
246 | // Adding Histograms: | |
247 | // 1. Linear Bg | |
248 | TStringToken names("mr1;mr1r;sr1;sr1r;ar1;br1;yr1;yr1int", ";"); | |
249 | TStringToken titles("Mass;Mass, Ratio Fit;Width;Width, Ratio Fit;a;b;Raw Yield; Raw Yield, integrated", ";"); | |
6aa36480 | 250 | while( names.NextToken() && titles.NextToken() ) { |
251 | new TH1D(names.Data(), titles.Data(), nPtBins,ptBinEdges); | |
252 | new TH1D(Form("%s_error", names.Data()), titles.Data(), nPtBins,ptBinEdges); | |
253 | } | |
47747ea2 | 254 | // 2. Quadratic Bg |
255 | TStringToken names2("mr2;mr2r;sr2;sr2r;ar2;br2;cr2;yr2;yr2int", ";"); | |
256 | TStringToken titles2("Mass;Mass, Ratio Fit;Width;Width, Ratio Fit;a;b;c;Raw Yield; Raw Yield, integrated", ";"); | |
6aa36480 | 257 | while( names2.NextToken() && titles2.NextToken() ) { |
258 | new TH1D(names2.Data(), titles2.Data(), nPtBins,ptBinEdges); | |
259 | new TH1D(Form("%s_error", names2.Data()), titles2.Data(), nPtBins,ptBinEdges); | |
260 | } | |
a9e9e511 | 261 | |
6aa36480 | 262 | TH1D* hMixRawRatio =new TH1D("hMixRawRatio", "ratio of statistics in RAW and Mixed", nPtBins, ptBinEdges); |
a9e9e511 | 263 | |
47747ea2 | 264 | // Pt slice loop |
265 | for(Int_t ptBin=1; ptBin<=nPtBins; ptBin++){ | |
266 | //gSystem->Sleep(2000); | |
267 | canvas->Clear(); | |
268 | canvas->Divide(2,2); | |
269 | canvas->cd(1); | |
270 | printf("pt bin %i, %3.1f to %3.1f, ", ptBin, ptBinEdges[ptBin-1], ptBinEdges[ptBin]); | |
a9e9e511 | 271 | |
47747ea2 | 272 | //TList* ptBinOutputList = outputList; |
273 | TAxis * ptAxis=hPi0->GetYaxis() ; | |
274 | Int_t imin=ptAxis->FindBin(ptBinEdges[ptBin-1]+0.0001); | |
275 | Int_t imax=ptAxis->FindBin(ptBinEdges[ptBin]+0.0001) -1; | |
276 | if( TMath::Abs( ptAxis->GetBinLowEdge(imin) - ptBinEdges[ptBin -1] ) >0.0001 ) { | |
277 | Printf("\nError: hPi0 lower bin edge (%f) different from ptBinEdges lower (%f)", ptAxis->GetBinLowEdge(imin), ptBinEdges[ptBin -1]); | |
278 | continue; | |
279 | } | |
280 | if( TMath::Abs( ptAxis->GetBinLowEdge(imax+1) - ptBinEdges[ptBin] ) >0.0001 ) { | |
281 | Printf("\nError: hPi0 upper bin edge (%f) different from ptBinEdges upper (%f)", ptAxis->GetBinLowEdge(imax+1), ptBinEdges[ptBin]); | |
282 | continue; | |
283 | } | |
a9e9e511 | 284 | |
47747ea2 | 285 | Double_t pt=(ptBinEdges[ptBin]+ptBinEdges[ptBin-1])/2. ; |
286 | Double_t dpt= ptBinEdges[ptBin] - ptBinEdges[ptBin-1]; | |
287 | ||
288 | TH1D * hPi0Proj = hPi0->ProjectionX(Form("pt%03i_hPi0",ptBin), imin, imax); | |
289 | hPi0Proj->Sumw2(); | |
290 | TH1D * hPi0ProjMix= hPi0Mix->ProjectionX(Form("pt%03i_hPi0Mix",ptBin), imin, imax) ; | |
291 | hPi0ProjMix->Sumw2(); | |
a9e9e511 | 292 | |
47747ea2 | 293 | const Int_t pi0Entries = hPi0Proj->Integral(hPi0Proj->FindBin(lowerMass), hPi0Proj->FindBin(upperMass)); |
294 | const Int_t mixEntries = hPi0ProjMix->Integral(hPi0Proj->FindBin(lowerMass), hPi0Proj->FindBin(upperMass)); | |
6aa36480 | 295 | if( pi0Entries >0) { |
296 | hMixRawRatio->SetBinContent(ptBin, mixEntries/pi0Entries); | |
297 | hMixRawRatio->SetBinError(ptBin, TMath::Sqrt(mixEntries/pi0Entries/pi0Entries + mixEntries*mixEntries/pi0Entries/pi0Entries/pi0Entries)); | |
298 | } | |
47747ea2 | 299 | printf("statistics in bin is %i, mixed %i, ", pi0Entries, mixEntries); |
300 | if( pi0Entries < 10 ) { | |
301 | Printf("to few entries"); | |
302 | continue; | |
303 | } | |
a9e9e511 | 304 | |
47747ea2 | 305 | const bool rebin = pi0Entries < 1000; |
306 | if( rebin && pi0Entries >= 500 ) { | |
307 | printf("rebin by factor 2, "); | |
308 | hPi0Proj->Rebin(2); | |
309 | hPi0ProjMix->Rebin(2); | |
310 | } else if( rebin && pi0Entries >= 200) { | |
311 | printf("rebin by factor 4, "); | |
312 | hPi0Proj->Rebin(4); | |
313 | hPi0ProjMix->Rebin(4); | |
314 | } else if( rebin && pi0Entries >= 10) { | |
315 | printf("rebin by factor 5, "); | |
316 | hPi0Proj->Rebin(5); | |
317 | hPi0ProjMix->Rebin(5); | |
318 | } | |
319 | std::cout << std::endl; | |
a9e9e511 | 320 | |
321 | hPi0Proj->SetTitle(Form("M_{#gamma#gamma}, %.1f<p_{T}<%.1f, PID=%s, %s", ptBinEdges[ptBin-1], ptBinEdges[ptBin], outBin.PID().Data(), outBin.Trigger().Data())); | |
47747ea2 | 322 | hPi0Proj->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})"); |
a9e9e511 | 323 | hPi0ProjMix->SetTitle(Form("M_{#gamma#gamma}^{Mix}, %.1f<p_{T}<%.1f, PID=%s, %s",ptBinEdges[ptBin-1],ptBinEdges[ptBin],outBin.PID().Data(), outBin.Trigger().Data())); |
47747ea2 | 324 | hPi0ProjMix->SetXTitle("M_{#gamma#gamma}^{Mix} (GeV/c^{2})"); |
325 | ||
326 | ||
327 | // Error Fix | |
328 | for(Int_t ib=1; ib<=hPi0Proj->GetNbinsX();ib++){if(hPi0Proj ->GetBinContent(ib)==0)hPi0Proj ->SetBinError(ib,1.);} | |
329 | for(Int_t ib=1; ib<=hPi0Proj->GetNbinsX();ib++){if(hPi0ProjMix->GetBinContent(ib)==0)hPi0ProjMix->SetBinError(ib,1.);} | |
a9e9e511 | 330 | |
47747ea2 | 331 | // Signal-Mix Ratio |
332 | canvas->cd(2); | |
333 | TH1D * hPi0Ratio = (TH1D*)hPi0Proj->Clone( Form("pt%03i_hPi0Ratio",ptBin) ) ; | |
334 | hPi0Ratio->SetTitle(Form("#frac{M_{#gamma#gamma}}{M_{#gamma#gamma}^{Mix}}, %.1f<p_{T}<%.1f GeV/c", ptBinEdges[ptBin-1], ptBinEdges[ptBin])); | |
335 | hPi0Ratio->Divide(hPi0ProjMix) ; | |
336 | hPi0Ratio->SetMarkerStyle(20) ; | |
337 | hPi0Ratio->SetMarkerSize(0.7) ; | |
338 | ||
20fb7b55 | 339 | // if(outBin.Centrality()==0) rangeMax=0.4 ; |
47747ea2 | 340 | // if(ptBin==1){ |
341 | // rangeMin=0.06 ; | |
342 | // rangeMax=0.25 ; | |
343 | // } | |
a9e9e511 | 344 | |
345 | ||
346 | ||
47747ea2 | 347 | // ================================================ |
348 | // Fit Pol1 ratio | |
349 | // ================================================ | |
350 | printf("Pol1 ratio Fit, "); | |
351 | canvas->cd(2); | |
a9e9e511 | 352 | |
47747ea2 | 353 | funcRatioFit1->SetParameters(0.001,0.136,0.0055,0.0002,-0.002) ; |
354 | funcRatioFit1->SetParLimits(0,0.000,1.000) ; | |
355 | funcRatioFit1->SetParLimits(1,lowerMass,upperMass) ; | |
356 | funcRatioFit1->SetParLimits(2,lowerWidth,upperWidth) ; | |
357 | ||
a9e9e511 | 358 | |
d1f44d99 | 359 | TFitResultPtr ratioFitResultPtr1 = hPi0Ratio->Fit(funcRatioFit1,"MSQ" ,"",rangeMin,rangeMax) ; |
47747ea2 | 360 | if( int(ratioFitResultPtr1) % 4000 ) // "More" error is acceptable |
d1f44d99 | 361 | ratioFitResultPtr1 = hPi0Ratio->Fit(funcRatioFit1,"MSQ" ,"",rangeMin,rangeMax) ; |
47747ea2 | 362 | |
363 | Int_t ratioFitError1 = ratioFitResultPtr1; | |
364 | ratioFitError1 = ratioFitError1 % 4000; // "More" error is acceptable | |
6aa36480 | 365 | ratioFitError1 = ratioFitError1 || TMath::Abs( funcRatioFit1->GetParameter(1) - lowerMass) < 0.0001 || TMath::Abs( funcRatioFit1->GetParameter(1) - upperMass) < 0.0001;// "center" mass converged to limit |
366 | ratioFitError1 = ratioFitError1 || funcRatioFit1->GetParError(1) > (upperMass - lowerMass)/2; // to large of an error | |
367 | ratioFitError1 = ratioFitError1 || TMath::Abs( funcRatioFit1->GetParameter(2) - lowerWidth) < 0.0001 || TMath::Abs( funcRatioFit1->GetParameter(2) - upperWidth) < 0.0001;// st. error converged to limit | |
368 | ratioFitError1 = ratioFitError1 || funcRatioFit1->GetParError(2) > (upperWidth - lowerWidth)/2; // to large of an error | |
a9e9e511 | 369 | |
47747ea2 | 370 | if( ratioFitError1) { |
371 | Printf("in ERROR, %i", ratioFitError1); | |
6aa36480 | 372 | ((TH1D*) output.GetHistogram("mr1r_error", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
373 | ((TH1D*) output.GetHistogram("mr1r_error", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; | |
374 | ((TH1D*) output.GetHistogram("sr1r_error", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; | |
375 | ((TH1D*) output.GetHistogram("sr1r_error", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; | |
376 | ((TH1D*) output.GetHistogram("ar1_error", outBin))->SetBinContent(ptBin,funcRatioFit1->GetParameter(3)) ; | |
377 | ((TH1D*) output.GetHistogram("ar1_error", outBin))->SetBinError (ptBin,funcRatioFit1->GetParError(3)) ; | |
378 | ((TH1D*) output.GetHistogram("br1_error", outBin))->SetBinContent(ptBin,funcRatioFit1->GetParameter(4)) ; | |
379 | ((TH1D*) output.GetHistogram("br1_error", outBin))->SetBinError (ptBin,funcRatioFit1->GetParError(4)) ; | |
a9e9e511 | 380 | } |
6aa36480 | 381 | if( !ratioFitError1 || ignoreErrors ) { |
47747ea2 | 382 | Printf("converged, status:%i, covMatrixStatus: %i", ratioFitResultPtr1->Status(), ratioFitResultPtr1->CovMatrixStatus()); |
6aa36480 | 383 | ((TH1D*) output.GetHistogram("mr1r", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
384 | ((TH1D*) output.GetHistogram("mr1r", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; | |
385 | ((TH1D*) output.GetHistogram("sr1r", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; | |
386 | ((TH1D*) output.GetHistogram("sr1r", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; | |
20fb7b55 | 387 | ((TH1D*) output.GetHistogram("ar1", outBin))->SetBinContent(ptBin,funcRatioFit1->GetParameter(3)) ; |
388 | ((TH1D*) output.GetHistogram("ar1", outBin))->SetBinError (ptBin,funcRatioFit1->GetParError(3)) ; | |
389 | ((TH1D*) output.GetHistogram("br1", outBin))->SetBinContent(ptBin,funcRatioFit1->GetParameter(4)) ; | |
390 | ((TH1D*) output.GetHistogram("br1", outBin))->SetBinError (ptBin,funcRatioFit1->GetParError(4)) ; | |
47747ea2 | 391 | } |
392 | ||
a9e9e511 | 393 | |
394 | ||
47747ea2 | 395 | // ================================================ |
396 | // Fit Pol2 ratio | |
397 | // ================================================ | |
398 | printf("Pol2 ratio Fit, "); | |
399 | if( ratioFitError1 ) { | |
400 | funcRatioFit2->SetParameters(0.001,0.136,0.0055,0.0002,-0.002, 0) ; | |
401 | funcRatioFit2->SetParLimits(0,0.000,1.000) ; | |
402 | funcRatioFit2->SetParLimits(1,lowerMass,upperMass) ; | |
403 | funcRatioFit2->SetParLimits(2,lowerWidth,upperWidth) ; | |
404 | } else { | |
405 | funcRatioFit2->SetParameters(funcRatioFit1->GetParameters()) ; | |
406 | funcRatioFit2->SetParameter(5, 0); | |
407 | funcRatioFit2->SetParLimits(0,0.000,1.000) ; | |
408 | funcRatioFit2->SetParLimits(1,lowerMass,upperMass) ; | |
409 | funcRatioFit2->SetParLimits(2,lowerWidth,upperWidth) ; | |
410 | } | |
d1f44d99 | 411 | TFitResultPtr ratioFitResultPtr2 = hPi0Ratio->Fit(funcRatioFit2,"+MSQ" ,"",rangeMin,rangeMax) ; |
47747ea2 | 412 | if( int(ratioFitResultPtr2) != 4000 ) // if error, "More" error is acceptable |
d1f44d99 | 413 | ratioFitResultPtr2 = hPi0Ratio->Fit(funcRatioFit2,"MSQ" ,"",rangeMin,rangeMax) ; |
47747ea2 | 414 | |
415 | Int_t ratioFitError2 = ratioFitResultPtr2; | |
416 | ratioFitError2 = ratioFitError2 % 4000; // "More" error is acceptable | |
6aa36480 | 417 | ratioFitError2 = ratioFitError2 || TMath::Abs( funcRatioFit2->GetParameter(1) - lowerMass) < 0.0001 || TMath::Abs( funcRatioFit2->GetParameter(1) - upperMass) < 0.0001; // "center" mass converged to limit |
418 | //ratioFitError2 = ratioFitError2 || funcRatioFit2->GetParError(1) > (upperMass - lowerMass)/2; // to large of an error | |
419 | ratioFitError2 = ratioFitError2 || TMath::Abs( funcRatioFit2->GetParameter(2) - lowerWidth) < 0.0001 || TMath::Abs( funcRatioFit2->GetParameter(2) - upperWidth) < 0.0001; // st. error converged to limit | |
420 | //ratioFitError2 = ratioFitError2 || funcRatioFit2->GetParError(2) > (upperWidth - lowerWidth)/2; // to large of an error | |
47747ea2 | 421 | |
422 | if( ratioFitError2) { | |
423 | Printf("in ERROR, %i", ratioFitError2); | |
6aa36480 | 424 | ((TH1D*) output.GetHistogram("mr2r_error", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
425 | ((TH1D*) output.GetHistogram("mr2r_error", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; | |
426 | ((TH1D*) output.GetHistogram("sr2r_error", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; | |
427 | ((TH1D*) output.GetHistogram("sr2r_error", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; | |
428 | ((TH1D*) output.GetHistogram("ar2_error", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(3)) ; | |
429 | ((TH1D*) output.GetHistogram("ar2_error", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(3)) ; | |
430 | ((TH1D*) output.GetHistogram("br2_error", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(4)) ; | |
431 | ((TH1D*) output.GetHistogram("br2_error", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(4)) ; | |
432 | ((TH1D*) output.GetHistogram("cr2_error", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(5)) ; | |
433 | ((TH1D*) output.GetHistogram("cr2_error", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(5)) ; | |
a9e9e511 | 434 | } |
6aa36480 | 435 | if( !ratioFitError2 || ignoreErrors ) { |
47747ea2 | 436 | Printf("converged, status:%i, covMatrixStatus: %i", ratioFitResultPtr2->Status(), ratioFitResultPtr2->CovMatrixStatus()); |
6aa36480 | 437 | ((TH1D*) output.GetHistogram("mr2r", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
438 | ((TH1D*) output.GetHistogram("mr2r", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; | |
439 | ((TH1D*) output.GetHistogram("sr2r", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; | |
440 | ((TH1D*) output.GetHistogram("sr2r", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; | |
20fb7b55 | 441 | ((TH1D*) output.GetHistogram("ar2", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(3)) ; |
442 | ((TH1D*) output.GetHistogram("ar2", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(3)) ; | |
443 | ((TH1D*) output.GetHistogram("br2", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(4)) ; | |
444 | ((TH1D*) output.GetHistogram("br2", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(4)) ; | |
445 | ((TH1D*) output.GetHistogram("cr2", outBin))->SetBinContent(ptBin,funcRatioFit2->GetParameter(5)) ; | |
446 | ((TH1D*) output.GetHistogram("cr2", outBin))->SetBinError (ptBin,funcRatioFit2->GetParError(5)) ; | |
47747ea2 | 447 | } |
47747ea2 | 448 | |
a9e9e511 | 449 | |
450 | ||
47747ea2 | 451 | // ================================================ |
452 | // Plot Ratio Fits | |
453 | // ================================================ | |
454 | hPi0Ratio->GetXaxis()->SetRangeUser(rangeMin, rangeMax); | |
455 | hPi0Ratio->Draw(); | |
456 | canvas->Update(); | |
457 | ||
458 | ||
a9e9e511 | 459 | |
47747ea2 | 460 | // ================================================ |
461 | // Pol1 Scaled Background Subtraction | |
462 | // ================================================ | |
463 | canvas->cd(3); | |
464 | Double_t intRangeMin = PeakPosition(pt)-3.*PeakWidth(pt) ; | |
465 | Double_t intRangeMax = PeakPosition(pt)+3.*PeakWidth(pt) ; | |
466 | Int_t intBinMin = hPi0Proj->GetXaxis()->FindBin(intRangeMin) ; | |
467 | Int_t intBinMax = hPi0Proj->GetXaxis()->FindBin(intRangeMax) ; | |
468 | Double_t mixInt = hPi0ProjMix->Integral(intBinMin,intBinMax); | |
469 | ||
6aa36480 | 470 | if( ! ratioFitError1 || true) { |
47747ea2 | 471 | printf("Pol1 BS Fit, "); |
472 | TH1D * hPi0MixScaledPol1 = (TH1D*)hPi0ProjMix->Clone(Form("pt%03i_hPi0MixScaledPol1", ptBin)) ; | |
473 | TH1D * hPi0BSPol1 = (TH1D*)hPi0Proj->Clone(Form("pt%03i_hPi0BSPol1", ptBin)) ; | |
474 | ||
475 | // Scale Mix by linear part of ratio, yielding approx background | |
476 | fbg1->SetParameters(funcRatioFit1->GetParameter(3), funcRatioFit1->GetParameter(4)); | |
477 | hPi0MixScaledPol1 ->Multiply(fbg1) ; | |
478 | hPi0BSPol1->Add(hPi0MixScaledPol1 ,-1.) ; | |
479 | ||
480 | Int_t binPi0 = hPi0BSPol1->FindBin(funcRatioFit1->GetParameter(1)); | |
481 | Int_t nWidPi0 = 2 * (Int_t) (funcRatioFit1->GetParameter(2)/hPi0BSPol1->GetBinWidth(1)); | |
482 | Int_t integral = TMath::Abs(hPi0BSPol1->Integral(binPi0-nWidPi0,binPi0+nWidPi0)); | |
483 | fgs->SetParameters(integral/5., funcRatioFit1->GetParameter(1), funcRatioFit1->GetParameter(2)) ; | |
484 | fgs->SetParLimits(0,0.,pi0Entries) ; | |
485 | fgs->SetParLimits(1,lowerMass,upperMass) ; | |
486 | fgs->SetParLimits(2,lowerWidth,upperWidth) ; | |
487 | ||
488 | // Fit | |
d1f44d99 | 489 | TFitResultPtr bs1FitResultPtr = hPi0BSPol1->Fit(fgs,"MSQ","",rangeMin,rangeMax) ; |
47747ea2 | 490 | if( int(bs1FitResultPtr) != 4000 ) // if error, "More" error is acceptable |
d1f44d99 | 491 | bs1FitResultPtr = hPi0BSPol1->Fit(fgs,"MSQ","",rangeMin,rangeMax) ; |
a9e9e511 | 492 | |
47747ea2 | 493 | Int_t bs1FitError = bs1FitResultPtr; |
494 | bs1FitError = bs1FitError % 4000; // "More" error is acceptable | |
6aa36480 | 495 | bs1FitError = bs1FitError || TMath::Abs( fgs->GetParameter(1) - lowerMass) < 0.0001 || TMath::Abs( fgs->GetParameter(1) - upperMass) < 0.0001; // "center" mass converged to limit |
496 | //bs1FitError = bs1FitError || fgs->GetParError(1) > (upperMass - lowerMass)/2; // to large of an error | |
497 | bs1FitError = bs1FitError || TMath::Abs( fgs->GetParameter(2) - lowerWidth) < 0.0001 || TMath::Abs( fgs->GetParameter(2) - upperWidth) < 0.0001; // st. error converged to limit | |
498 | //bs1FitError = bs1FitError || fgs->GetParError(2) > (upperWidth - lowerWidth)/2; // to large of an error | |
47747ea2 | 499 | |
6aa36480 | 500 | Double_t y=fgs->GetParameter(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ; |
501 | Double_t ey=fgs->GetParError(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ; | |
502 | Double_t npiInt = hPi0BSPol1->Integral(intBinMin,intBinMax) ; | |
503 | Double_t norm = fbg1->GetParameter(0) ; | |
504 | Double_t normErr= fbg1->GetParError(0) ; | |
47747ea2 | 505 | if( bs1FitError) { |
506 | Printf("in ERROR, %i", bs1FitError); | |
6aa36480 | 507 | ((TH1D*) output.GetHistogram("mr1_error", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
508 | ((TH1D*) output.GetHistogram("mr1_error", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; | |
509 | ((TH1D*) output.GetHistogram("sr1_error", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; | |
510 | ((TH1D*) output.GetHistogram("sr1_error", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; | |
511 | ||
512 | ((TH1D*) output.GetHistogram("yr1_error", outBin))->SetBinContent(ptBin,y/dpt) ; | |
513 | ((TH1D*) output.GetHistogram("yr1_error", outBin))->SetBinError(ptBin,ey/dpt) ; | |
514 | if(npiInt>0.){ | |
515 | ((TH1D*) output.GetHistogram("yr1int_error", outBin))->SetBinContent(ptBin,npiInt/dpt) ; | |
516 | ((TH1D*) output.GetHistogram("yr1int_error", outBin))->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*mixInt + normErr*normErr*mixInt*mixInt + norm*norm*mixInt)/dpt) ; | |
517 | } | |
a9e9e511 | 518 | } |
6aa36480 | 519 | if( !bs1FitError || ignoreErrors ) { |
47747ea2 | 520 | Printf("converged, status:%i, covMatrixStatus: %i", bs1FitResultPtr->Status(), bs1FitResultPtr->CovMatrixStatus()); |
20fb7b55 | 521 | ((TH1D*) output.GetHistogram("mr1", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
522 | ((TH1D*) output.GetHistogram("mr1", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; | |
523 | ((TH1D*) output.GetHistogram("sr1", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; | |
524 | ((TH1D*) output.GetHistogram("sr1", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; | |
47747ea2 | 525 | |
20fb7b55 | 526 | ((TH1D*) output.GetHistogram("yr1", outBin))->SetBinContent(ptBin,y/dpt) ; |
20fb7b55 | 527 | ((TH1D*) output.GetHistogram("yr1", outBin))->SetBinError(ptBin,ey/dpt) ; |
47747ea2 | 528 | if(npiInt>0.){ |
20fb7b55 | 529 | ((TH1D*) output.GetHistogram("yr1int", outBin))->SetBinContent(ptBin,npiInt/dpt) ; |
530 | ((TH1D*) output.GetHistogram("yr1int", outBin))->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*mixInt + normErr*normErr*mixInt*mixInt + norm*norm*mixInt)/dpt) ; | |
47747ea2 | 531 | // maybe we should use TH1::IntegralAndError |
532 | } | |
533 | } | |
534 | // Ploting | |
535 | hPi0BSPol1->SetAxisRange(rangeMin, rangeMax); | |
536 | hPi0BSPol1->SetMaximum(hPi0BSPol1->GetMaximum()*1.5) ; | |
537 | hPi0BSPol1->SetMinimum(hPi0BSPol1->GetMinimum()*0.9) ; | |
538 | hPi0BSPol1->SetMarkerStyle(20) ; | |
539 | hPi0BSPol1->SetMarkerSize(0.7) ; | |
540 | hPi0BSPol1->Draw(); | |
541 | canvas->Update(); | |
542 | //printf("end pt"); | |
543 | } | |
a9e9e511 | 544 | |
545 | ||
47747ea2 | 546 | // ================================================ |
a9e9e511 | 547 | // Pol2 Scaled Background Subtraction |
47747ea2 | 548 | // ================================================ |
549 | canvas->cd(4); | |
550 | fbg2->SetParameters(funcRatioFit2->GetParameter(3),funcRatioFit2->GetParameter(4),funcRatioFit2->GetParameter(5)); | |
6aa36480 | 551 | if( ! ratioFitError2 || true) { |
47747ea2 | 552 | printf("Pol1 Scaled Background Subtraction, "); |
553 | TH1D * hPi0MixScaledPol2 = (TH1D*)hPi0ProjMix->Clone(Form("pt%03i_hPi0MixScaledPol2", ptBin)) ; | |
554 | TH1D * hPi0BSPol2 = (TH1D*)hPi0Proj ->Clone(Form("pt%03i_hPi0BSPol2", ptBin)) ; | |
a9e9e511 | 555 | |
47747ea2 | 556 | hPi0MixScaledPol2->Multiply(fbg2) ; |
557 | hPi0BSPol2 ->Add(hPi0MixScaledPol2,-1.) ; | |
6aa36480 | 558 | hPi0BSPol2->SetOption(); |
47747ea2 | 559 | |
560 | Int_t binPi0 = hPi0BSPol2->FindBin(funcRatioFit2->GetParameter(1)); | |
561 | Int_t nWidPi0 = 2 * (Int_t) (funcRatioFit2->GetParameter(2)/hPi0BSPol2->GetBinWidth(1)); | |
562 | Int_t integral = TMath::Abs(hPi0BSPol2->Integral(binPi0-nWidPi0,binPi0+nWidPi0)); | |
563 | fgs->SetParameters(integral/5., funcRatioFit2->GetParameter(1), funcRatioFit2->GetParameter(2)) ; | |
564 | fgs->SetParLimits(0,0.,pi0Entries) ; | |
565 | fgs->SetParLimits(1,lowerMass,upperMass) ; | |
566 | fgs->SetParLimits(2,lowerWidth,upperWidth) ; | |
567 | ||
568 | // Fit | |
d1f44d99 | 569 | TFitResultPtr bs2FitResultPtr = hPi0BSPol2->Fit(fgs,"MSQ","",rangeMin,rangeMax) ; |
47747ea2 | 570 | if( int(bs2FitResultPtr) != 4000 ) // if error, "More" error is acceptable |
d1f44d99 | 571 | bs2FitResultPtr = hPi0BSPol2->Fit(fgs,"MSQ","",rangeMin,rangeMax) ; |
a9e9e511 | 572 | |
47747ea2 | 573 | Int_t bs2FitError = bs2FitResultPtr; |
574 | bs2FitError = bs2FitError % 4000; // "More" error is acceptable | |
6aa36480 | 575 | bs2FitError = bs2FitError || TMath::Abs( fgs->GetParameter(1) - lowerMass) < 0.0001 || TMath::Abs( fgs->GetParameter(1) - upperMass) < 0.0001; // "center" mass converged to limit |
47747ea2 | 576 | // bs2FitError = bs2FitError || fgs->GetParError(1) > (upperMass - lowerMass)/2; // to large of an error |
6aa36480 | 577 | bs2FitError = bs2FitError || TMath::Abs( fgs->GetParameter(2) - lowerWidth) < 0.0001 || TMath::Abs( fgs->GetParameter(2) - upperWidth) < 0.0001; // st. error converged to limit |
47747ea2 | 578 | // bs2FitError = bs2FitError || fgs->GetParError(2) > (upperWidth - lowerWidth)/2; // to large of an error |
579 | ||
6aa36480 | 580 | Double_t y=fgs->GetParameter(0)/hPi0BSPol2->GetXaxis()->GetBinWidth(1) ; |
581 | Double_t ey=fgs->GetParError(0)/hPi0BSPol2->GetXaxis()->GetBinWidth(1) ; | |
582 | Double_t npiInt = hPi0BSPol2->Integral(intBinMin,intBinMax) ; | |
583 | Double_t norm = fbg2->GetParameter(0) ; | |
584 | Double_t normErr= fbg2->GetParError(0) ; | |
47747ea2 | 585 | if( bs2FitError ) { |
586 | Printf("in ERROR, %i", bs2FitError); | |
6aa36480 | 587 | ((TH1D*) output.GetHistogram("mr2_error", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
588 | ((TH1D*) output.GetHistogram("mr2_error", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; | |
589 | ((TH1D*) output.GetHistogram("sr2_error", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; | |
590 | ((TH1D*) output.GetHistogram("sr2_error", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; | |
591 | ||
592 | ((TH1D*) output.GetHistogram("yr2_error", outBin))->SetBinContent(ptBin,y/dpt) ; | |
593 | ((TH1D*) output.GetHistogram("yr2_error", outBin))->SetBinError(ptBin,ey/dpt) ; | |
594 | if(npiInt>0.){ | |
595 | ((TH1D*) output.GetHistogram("yr2int_error", outBin))->SetBinContent(ptBin,npiInt/dpt) ; | |
596 | ((TH1D*) output.GetHistogram("yr2int_error", outBin))->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*mixInt + normErr*normErr*mixInt*mixInt + norm*norm*mixInt)/dpt) ; | |
597 | // maybe we should use TH1::IntegralAndError | |
598 | } | |
a9e9e511 | 599 | } |
6aa36480 | 600 | if( !bs2FitError || ignoreErrors ) { |
47747ea2 | 601 | Printf("converged, status:%i, covMatrixStatus: %i", bs2FitResultPtr->Status(), bs2FitResultPtr->CovMatrixStatus()); |
20fb7b55 | 602 | ((TH1D*) output.GetHistogram("mr2", outBin))->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
603 | ((TH1D*) output.GetHistogram("mr2", outBin))->SetBinError (ptBin,fgs->GetParError(1) ) ; | |
604 | ((TH1D*) output.GetHistogram("sr2", outBin))->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; | |
605 | ((TH1D*) output.GetHistogram("sr2", outBin))->SetBinError (ptBin,fgs->GetParError(2) ) ; | |
47747ea2 | 606 | |
20fb7b55 | 607 | ((TH1D*) output.GetHistogram("yr2", outBin))->SetBinContent(ptBin,y/dpt) ; |
20fb7b55 | 608 | ((TH1D*) output.GetHistogram("yr2", outBin))->SetBinError(ptBin,ey/dpt) ; |
47747ea2 | 609 | if(npiInt>0.){ |
20fb7b55 | 610 | ((TH1D*) output.GetHistogram("yr2int", outBin))->SetBinContent(ptBin,npiInt/dpt) ; |
611 | ((TH1D*) output.GetHistogram("yr2int", outBin))->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*mixInt + normErr*normErr*mixInt*mixInt + norm*norm*mixInt)/dpt) ; | |
47747ea2 | 612 | // maybe we should use TH1::IntegralAndError |
613 | } | |
614 | } | |
615 | // Plotting | |
616 | hPi0BSPol2->SetAxisRange(rangeMin, rangeMax); | |
617 | hPi0BSPol2->SetMaximum(hPi0BSPol2->GetMaximum()*1.5) ; | |
618 | hPi0BSPol2->SetMinimum(hPi0BSPol2->GetMinimum()*0.9) ; | |
619 | hPi0BSPol2->SetMarkerStyle(20) ; | |
620 | hPi0BSPol2->SetMarkerSize(0.7) ; | |
621 | hPi0BSPol2->Draw(); | |
622 | canvas->Update(); | |
623 | ||
624 | } | |
a9e9e511 | 625 | |
47747ea2 | 626 | canvas->cd(1); |
627 | TH1D * hPi0MixScaled = (TH1D*)hPi0ProjMix ->Clone(Form("%sScaled", hPi0Proj->GetName())) ; | |
628 | //hPi0MixScaled->Scale(double(pi0Entries)/mixEntries); | |
629 | hPi0MixScaled->Scale(fbg2->Eval(0.134)); | |
630 | hPi0MixScaled->SetLineColor(kRed); | |
631 | hPi0MixScaled->SetTitle(Form("%s, Scaled", hPi0Proj->GetName())); | |
632 | hPi0Proj->SetAxisRange(rangeMin, 1.); | |
633 | //hPi0Proj->SetMaximum(TMath::Max(hPi0Proj->GetMaximum(), hPi0ProjMix->GetMaximum())*1.2); | |
634 | hPi0Proj->SetMinimum(0); | |
635 | hPi0Proj->Draw(); | |
636 | hPi0MixScaled->Draw("same"); | |
a9e9e511 | 637 | |
47747ea2 | 638 | canvas->Update(); |
a9e9e511 | 639 | |
20fb7b55 | 640 | canvas->Print(Form("imgs/%s_ptBin:%03i.pdf", outBin.Key().Data(), ptBin)); |
641 | canvas->Print(Form("imgs/%s_ptBin:%03i.png", outBin.Key().Data(), ptBin)); | |
a9e9e511 | 642 | |
47747ea2 | 643 | std::cout << std::endl; |
644 | }// end of Pt slice loop | |
645 | ||
646 | ||
647 | //Normalize by the number of events | |
648 | Int_t cMin=0, cMax=0; | |
d1f44d99 | 649 | if( input.Bin().Trigger().EqualTo("kCentral") ) |
a9e9e511 | 650 | switch(outBin.Centrality()) { |
47747ea2 | 651 | case 0: cMin = 1; cMax = 5; break; |
652 | case 1: cMin = 6; cMax = 10; break; | |
6aa36480 | 653 | case -1: cMin = 1; cMax = 10; break; |
47747ea2 | 654 | default: Printf("ERROR: cent bin not defined for trigger"); |
655 | } | |
d1f44d99 | 656 | else if( input.Bin().Trigger().EqualTo("kSemiCentral") ) |
a9e9e511 | 657 | switch(outBin.Centrality()) { |
47747ea2 | 658 | case 0: cMin = 11; cMax = 20; break; |
659 | case 1: cMin = 21; cMax = 30; break; | |
660 | case 2: cMin = 31; cMax = 40; break; | |
661 | case 3: cMin = 41; cMax = 50; break; | |
6aa36480 | 662 | case -2: cMin = 11; cMax = 20; break; |
663 | case -3: cMin = 21; cMax = 30; break; | |
664 | case -4: cMin = 31; cMax = 40; break; | |
665 | case -5: cMin = 41; cMax = 50; break; | |
47747ea2 | 666 | default: Printf("ERROR: cent bin not defined for trigger"); |
667 | } | |
d1f44d99 | 668 | else if ( input.Bin().Trigger().EqualTo("kMB") || input.Bin().Trigger().EqualTo("kPHOSPb") ) |
20fb7b55 | 669 | switch(outBin.Centrality()) { |
47747ea2 | 670 | case 0: cMin = 1; cMax = 5; break; |
671 | case 1: cMin = 6; cMax = 10; break; | |
672 | case 2: cMin = 11; cMax = 20; break; | |
673 | case 3: cMin = 21; cMax = 30; break; | |
674 | case 4: cMin = 31; cMax = 40; break; | |
675 | case 5: cMin = 41; cMax = 50; break; | |
676 | case 6: cMin = 51; cMax = 80; break; | |
677 | case -10: cMin=1; cMax = 80; break; | |
6aa36480 | 678 | case -1: cMin = 1; cMax = 10; break; |
679 | case -2: cMin = 11; cMax = 20; break; | |
680 | case -3: cMin = 21; cMax = 30; break; | |
681 | case -4: cMin = 31; cMax = 40; break; | |
682 | case -5: cMin = 41; cMax = 50; break; | |
683 | case -6: cMin = 51; cMax = 80; break; | |
47747ea2 | 684 | default: Printf("ERROR: cent bin not defined for trigger"); |
685 | } | |
686 | else | |
d1f44d99 | 687 | Printf("ERROR: cent bins not defined for trigger, %s", input.Bin().Trigger().Data()); |
47747ea2 | 688 | |
689 | Double_t nevents = hCentralityX->Integral(cMin,cMax); | |
690 | if ( nevents > 0.9 ) { | |
20fb7b55 | 691 | ((TH1D*) output.GetHistogram("yr1", outBin)) ->Scale(1./nevents) ; |
692 | ((TH1D*) output.GetHistogram("yr1int", outBin)) ->Scale(1./nevents) ; | |
693 | ((TH1D*) output.GetHistogram("yr2", outBin)) ->Scale(1./nevents) ; | |
694 | ((TH1D*) output.GetHistogram("yr2int", outBin)) ->Scale(1./nevents) ; | |
a9e9e511 | 695 | |
6aa36480 | 696 | ((TH1D*) output.GetHistogram("yr1_error", outBin)) ->Scale(1./nevents) ; |
697 | ((TH1D*) output.GetHistogram("yr1int_error", outBin)) ->Scale(1./nevents) ; | |
698 | ((TH1D*) output.GetHistogram("yr2_error", outBin)) ->Scale(1./nevents) ; | |
699 | ((TH1D*) output.GetHistogram("yr2int_error", outBin)) ->Scale(1./nevents) ; | |
47747ea2 | 700 | } else { |
701 | Printf("WARNING: non positive nEvents in centrality range, cMin:%d, cMax:%d, nEvents:%f", cMin, cMax, nevents ); | |
702 | ||
703 | } | |
704 | ||
a9e9e511 | 705 | |
47747ea2 | 706 | // store to file |
707 | // TFile* outputFile = TFile::Open(saveToFileName.Data(), "UPDATE"); | |
708 | //outputList->Write(input.KeySuggestion(), TObject::kSingleKey); | |
709 | //outputList->Write(); | |
20fb7b55 | 710 | // outputFile->Write(); |
711 | // outputFile->Close(); | |
712 | // delete outputFile; | |
47747ea2 | 713 | // delete list and content from memory |
714 | //outputList->SetOwner(kTRUE); | |
715 | //outputList->Delete("slow"); | |
716 | //delete outputList; | |
20fb7b55 | 717 | //output.Write(); |
6aa36480 | 718 | delete canvas; |
47747ea2 | 719 | } |
720 | ||
721 | ||
722 | //----------------------------------------------------------------------------- | |
723 | double GetPtBin(int bin){ | |
724 | // recusive function used by MakePtBins | |
725 | ||
726 | if( bin==0 ) | |
727 | return 1.; | |
728 | ||
729 | // return GetPtBin(bin-1) * 1.1; | |
730 | ||
731 | // if ( bin % 2 ) | |
732 | // return GetPtBin(bin-1) + 0.4; | |
733 | // else | |
734 | // return GetPtBin(bin-1) + 0.2; | |
735 | ||
736 | double previousBin = GetPtBin(bin-1); | |
737 | double linInc = 1.; | |
738 | double threshold = 5.; | |
739 | double logFact = 1 + linInc/threshold; | |
740 | if ( previousBin < threshold ) // linear | |
741 | return double(int((previousBin + linInc)*10))/10; | |
742 | else { // logarithmic | |
743 | return double(int((previousBin * logFact)*10))/10; | |
744 | } | |
745 | } | |
746 | ||
747 | //----------------------------------------------------------------------------- | |
748 | void MakePtBins() { | |
749 | // function for setting Pt bins | |
750 | ||
751 | int bin = -1; | |
752 | do { | |
753 | ++bin; | |
754 | ptBinEdges[bin] = GetPtBin(bin); | |
755 | } while(ptBinEdges[bin] < 40.); | |
756 | nPtBins = bin -2; | |
757 | ||
758 | printf("Making Pt Bins:\n"); | |
759 | for(int b=0; b < nPtBins+1; ++b) | |
760 | printf("%.1f, ", ptBinEdges[b]); | |
761 | printf("\n N. Bins: %d\n", nPtBins); | |
762 | ||
763 | ||
764 | // for(int bin = 0; bin <= nPtBins; ++bin){ | |
765 | // ptBinEdges[bin] = GetPtBin(bin); | |
766 | // printf("%.1f, ", ptBinEdges[bin]); | |
767 | // } | |
768 | // printf("\n"); | |
769 | } | |
770 | ||
771 | //----------------------------------------------------------------------------- | |
772 | Double_t CGausPol1(Double_t * x, Double_t * par){ | |
773 | //Parameterization of Real/Mixed ratio | |
774 | Double_t m=par[1] ; | |
775 | Double_t s=par[2] ; | |
776 | Double_t dx=(x[0]-m)/s ; | |
777 | return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean) ; | |
778 | } | |
779 | ||
780 | //----------------------------------------------------------------------------- | |
781 | Double_t CGausPol2(Double_t * x, Double_t * par){ | |
782 | //Another parameterization of Real/Mixed ratio7TeV/README | |
783 | Double_t m=par[1] ; | |
784 | Double_t s=par[2] ; | |
785 | Double_t dx=(x[0]-m)/s ; | |
786 | return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean)+par[5]*(x[0]-kMean)*(x[0]-kMean) ; | |
787 | } | |
788 | //----------------------------------------------------------------------------- | |
789 | Double_t CGausPol0(Double_t * x, Double_t * par){ | |
790 | //Parameterizatin of signal | |
791 | Double_t m=par[1] ; | |
792 | Double_t s=par[2] ; | |
793 | Double_t dx=(x[0]-m)/s ; | |
794 | return par[0]*exp(-dx*dx/2.)/TMath::Sqrt(TMath::TwoPi())/s+par[3] ; | |
795 | } | |
796 | //----------------------------------------------------------------------------- | |
797 | Double_t CPol1(Double_t * x, Double_t * par){ | |
798 | //Normalizatino of Mixed | |
799 | return par[0]+par[1]*(x[0]-kMean) ; | |
800 | } | |
801 | //----------------------------------------------------------------------------- | |
802 | Double_t CPol2(Double_t * x, Double_t * par){ | |
803 | //Another normalization of Mixed | |
804 | return par[0]+par[1]*(x[0]-kMean)+par[2]*(x[0]-kMean)*(x[0]-kMean) ; | |
805 | } | |
806 | ||
a9e9e511 | 807 | |
47747ea2 | 808 | // Input Definitions |
809 | TFile* Input::fFile = 0; | |
a9e9e511 | 810 | Input::Input(const TString& fileName, const RawProduction::TriggerBin& inputBin, TString listPath) |
d1f44d99 | 811 | : fList(0x0), fBin(inputBin.Trigger()) |
47747ea2 | 812 | { |
813 | // File | |
814 | if(fFile && !fileName.EqualTo(fFile->GetName())){ | |
815 | fFile->Close(); | |
816 | delete fFile; | |
817 | fFile = 0x0; | |
818 | } else if(! fFile) { | |
819 | Printf("Opening %s", fileName.Data()); | |
a9e9e511 | 820 | fFile = TFile::Open(fileName.Data(), "READ"); |
47747ea2 | 821 | } |
a9e9e511 | 822 | |
47747ea2 | 823 | if( listPath.EqualTo("") ) { |
824 | char cstr[256] = ""; | |
d1f44d99 | 825 | sprintf(cstr, "PHOSPi0Flow_%s/PHOSPi0Flow_%sCoutput1", fBin.Trigger().Data(), fBin.Trigger().Data()); |
47747ea2 | 826 | listPath = cstr; |
827 | } | |
a9e9e511 | 828 | |
47747ea2 | 829 | Printf("Getting list, %s", listPath.Data()); |
830 | fFile->GetObject(listPath.Data(), fList); | |
a9e9e511 | 831 | if( !fList ) |
47747ea2 | 832 | Printf("ERROR: list not found"); |
833 | } | |
a9e9e511 | 834 | |
47747ea2 | 835 | TH1* Input::GetHistogram(const char* name){ |
836 | TObject* obj = fList->FindObject(name); | |
837 | TH1* hist = dynamic_cast<TH1*> (obj); | |
838 | if( ! hist) | |
839 | Printf("MakePi0FitInput::GetHistogram: Error, could not find object of name: %s or cast to hist", name); | |
840 | return hist; | |
841 | } | |
842 | ||
d1f44d99 | 843 | //OutputBin Definitions |
a53ea5d3 | 844 | TriggerBin::TriggerBin(const TString& trigger) |
3374aa5b | 845 | : fTrigger(trigger), fDir(trigger), fKey(trigger) |
d1f44d99 | 846 | { } |
847 | ||
a53ea5d3 | 848 | TriggerBin::TriggerBin(const char* trigger) |
3374aa5b | 849 | : fTrigger(trigger), fDir(trigger), fKey(trigger) |
d1f44d99 | 850 | { } |
851 | ||
a53ea5d3 | 852 | TriCenPidBin::TriCenPidBin(Int_t centrality, const TString& pid, const TString& trigger) |
853 | : TriggerBin(trigger), fCentrality(centrality), fPID(pid) | |
d1f44d99 | 854 | { |
3374aa5b | 855 | fDir.Form("%s/c%03i/%s", trigger.Data(), centrality, pid.Data()); |
856 | fKey.Form("%s_c%03i_%s", trigger.Data(), centrality, pid.Data()); | |
d1f44d99 | 857 | } |
858 | ||
859 | Output::Output(const TString& fileName, const char* options) | |
6aa36480 | 860 | : fFile(0x0) |
d1f44d99 | 861 | { |
862 | fFile = TFile::Open(fileName.Data(), options); | |
d1f44d99 | 863 | } |
a9e9e511 | 864 | |
a53ea5d3 | 865 | void Output::SetDir(const TriggerBin& inBin) |
d1f44d99 | 866 | { |
3374aa5b | 867 | TDirectory* dir = (TDirectoryFile*) fFile; |
868 | TStringToken dirs(inBin.Dir(), "/"); | |
869 | while( dirs.NextToken() ) { | |
870 | Printf("%s", dirs.Data()); | |
871 | TDirectory* ndir = dir->GetDirectory(dirs.Data()); | |
872 | if(!ndir) | |
873 | dir = dir->mkdir(dirs.Data()); | |
874 | else | |
875 | dir = ndir; | |
d1f44d99 | 876 | } |
3374aa5b | 877 | dir->cd(); |
d1f44d99 | 878 | } |
a9e9e511 | 879 | |
a53ea5d3 | 880 | TH1* Output::GetHistogram(const TString& name, const RawProduction::TriggerBin& inBin) |
d1f44d99 | 881 | { |
3374aa5b | 882 | TDirectory* dir = fFile->GetDirectory(inBin.Dir().Data(), true); |
a9e9e511 | 883 | TH1* hist = dynamic_cast<TH1*>( dir->Get(name.Data()) ); |
d1f44d99 | 884 | if( hist ) |
885 | return hist; | |
886 | else { | |
a53ea5d3 | 887 | Printf("ERROR: Output::GetHistogram: %s could not be found", name.Data()); |
d1f44d99 | 888 | return 0x0; |
889 | } | |
890 | } | |
baacba22 | 891 | |
892 | TH1* Output::GetHistogram(const TString& name) { | |
893 | TH1* hist = dynamic_cast<TH1*>( fFile->Get(name.Data()) ); | |
894 | if( hist ) | |
895 | return hist; | |
896 | else { | |
897 | Printf("ERROR: Output::GetHistogram: %s could not be found", name.Data()); | |
898 | return 0x0; | |
899 | } | |
900 | } | |
901 | ||
6aa36480 | 902 | |
903 | ||
d1f44d99 | 904 | |
905 | void Output::Write() | |
906 | { | |
6aa36480 | 907 | fFile->Write(); |
d1f44d99 | 908 | } |
909 | ||
910 | ||
a9e9e511 | 911 | |
912 | ||
47747ea2 | 913 | TH1* GetHistogram_cent(Input& input, const TString& name, int centrality) |
914 | { | |
915 | // Getter (and merger) for histograms following the naming patern of %s_cen%i | |
916 | // | |
917 | // For certain negeative values, function is defined to merge centrality bins. | |
918 | // -1 0-10% | |
919 | // -2 10-20% | |
920 | // -3 20-30% | |
921 | // -4 30-40% | |
922 | // -5 40-50% | |
04f85d20 | 923 | // -6 50-80% |
924 | // -7 50-60% | |
925 | // -8 60-70% | |
926 | // -9 70-80% | |
a9e9e511 | 927 | // -10 0-80% |
04f85d20 | 928 | // -11 10-50% |
47747ea2 | 929 | if( centrality >= 0 ) { |
930 | char cname[256] = ""; | |
931 | sprintf(cname, "%s_cen%i", name.Data(), centrality); | |
932 | input.GetHistogram(cname); | |
933 | } | |
a9e9e511 | 934 | |
47747ea2 | 935 | TH1* hist = 0x0; |
04f85d20 | 936 | if( 1 == centBinVersion ) { |
937 | if( input.Bin().Trigger().EqualTo("kMB") || input.Bin().Trigger().EqualTo("kPHOSPb") ) { | |
938 | switch(centrality) { | |
939 | case -10: hist = MergeHistogram_cent(input, name, centrality, 0, 7); break; | |
940 | case -1: hist = MergeHistogram_cent(input, name, centrality, 0, 2); break; | |
941 | case -2: hist = MergeHistogram_cent(input, name, centrality, 2, 3); break; | |
942 | case -3: hist = MergeHistogram_cent(input, name, centrality, 3, 4); break; | |
943 | case -4: hist = MergeHistogram_cent(input, name, centrality, 4, 5); break; | |
944 | case -5: hist = MergeHistogram_cent(input, name, centrality, 5, 6); break; | |
945 | case -6: hist = MergeHistogram_cent(input, name, centrality, 6, 7); break; | |
946 | } | |
947 | } else if ( input.Bin().Trigger().EqualTo("kCentral") ) { | |
948 | switch( centrality ) { | |
949 | case -1: return MergeHistogram_cent(input, name, centrality, 0, 2); break; | |
950 | } | |
951 | } else if ( input.Bin().Trigger().EqualTo("kSemiCentral") ) { | |
952 | switch( centrality ) { | |
953 | case -2: return MergeHistogram_cent(input, name, centrality, 0, 1); break; | |
954 | case -3: return MergeHistogram_cent(input, name, centrality, 1, 2); break; | |
955 | case -4: return MergeHistogram_cent(input, name, centrality, 2, 3); break; | |
956 | case -5: return MergeHistogram_cent(input, name, centrality, 3, 4); break; | |
957 | } | |
47747ea2 | 958 | } |
04f85d20 | 959 | } |
960 | else if ( 2 == centBinVersion ) { | |
961 | if( input.Bin().Trigger().EqualTo("kMB") || input.Bin().Trigger().EqualTo("kPHOSPb") ) { | |
962 | switch(centrality) { | |
963 | case -10: hist = MergeHistogram_cent(input, name, centrality, 0, 8); break; | |
964 | case -1: hist = MergeHistogram_cent(input, name, centrality, 0, 1); break; | |
965 | case -2: hist = MergeHistogram_cent(input, name, centrality, 1, 2); break; | |
966 | case -3: hist = MergeHistogram_cent(input, name, centrality, 2, 3); break; | |
967 | case -4: hist = MergeHistogram_cent(input, name, centrality, 3, 4); break; | |
968 | case -5: hist = MergeHistogram_cent(input, name, centrality, 4, 5); break; | |
969 | case -6: hist = MergeHistogram_cent(input, name, centrality, 5, 8); break; | |
970 | case -7: hist = MergeHistogram_cent(input, name, centrality, 5, 6); break; | |
971 | case -8: hist = MergeHistogram_cent(input, name, centrality, 6, 7); break; | |
972 | case -9: hist = MergeHistogram_cent(input, name, centrality, 7, 8); break; | |
973 | } | |
974 | } else if ( input.Bin().Trigger().EqualTo("kCentral") ) { | |
975 | switch( centrality ) { | |
976 | case -1: return MergeHistogram_cent(input, name, centrality, 0, 4); break; | |
977 | } | |
978 | } else if ( input.Bin().Trigger().EqualTo("kSemiCentral") ) { | |
979 | switch( centrality ) { | |
980 | case -2: return MergeHistogram_cent(input, name, centrality, 0, 5); break; | |
981 | case -3: return MergeHistogram_cent(input, name, centrality, 5, 6); break; | |
982 | case -4: return MergeHistogram_cent(input, name, centrality, 6, 7); break; | |
983 | case -5: return MergeHistogram_cent(input, name, centrality, 7, 8); break; | |
984 | } | |
47747ea2 | 985 | } |
04f85d20 | 986 | |
47747ea2 | 987 | } |
988 | // in case not defined above | |
989 | if( ! hist ) { | |
d1f44d99 | 990 | Printf("ERROR:GetHistogram_cent: %i not possible for %s trigger", centrality, input.Bin().Trigger().Data()); |
47747ea2 | 991 | return 0x0; |
992 | } | |
a9e9e511 | 993 | |
47747ea2 | 994 | switch(centrality) { |
995 | case -10: hist->SetTitle( Form("%s, 0-80%% centrality", hist->GetTitle())); break; | |
04f85d20 | 996 | case -11: hist->SetTitle( Form("%s,10-50%% centrality", hist->GetTitle())); break; |
997 | case -1: hist->SetTitle(Form("%s, 0-10%% centrality", hist->GetTitle())); break; | |
47747ea2 | 998 | case -2: hist->SetTitle(Form("%s, 10-20%% centrality", hist->GetTitle())); break; |
999 | case -3: hist->SetTitle(Form("%s, 20-30%% centrality", hist->GetTitle())); break; | |
1000 | case -4: hist->SetTitle(Form("%s, 30-40%% centrality", hist->GetTitle())); break; | |
1001 | case -5: hist->SetTitle(Form("%s, 40-50%% centrality", hist->GetTitle())); break; | |
04f85d20 | 1002 | case -6: hist->SetTitle( Form("%s, 50-80%% centrality", hist->GetTitle())); break; |
1003 | case -7: hist->SetTitle(Form("%s, 50-60%% centrality", hist->GetTitle())); break; | |
1004 | case -8: hist->SetTitle(Form("%s, 60-70%% centrality", hist->GetTitle())); break; | |
1005 | case -9: hist->SetTitle(Form("%s, 70-80%% centrality", hist->GetTitle())); break; | |
47747ea2 | 1006 | } |
1007 | return hist; | |
1008 | } | |
a9e9e511 | 1009 | |
47747ea2 | 1010 | TH1* MergeHistogram_cent(Input& input, const TString& name, int newCentIndex, int fromCentIndex, int toCentIndex) |
1011 | { | |
1012 | // Merger (All cent) for histograms following the naming patern of %s_cen%i, from including to excluding | |
1013 | // | |
1014 | // Merges centralites bins into one histogram, from including to excluding, and names the histogram using the patern above. | |
1015 | // If an histogram with that name Allready exist in the current directory (gDirectory), then no merge | |
1016 | // occurs and this hist. is simply returned. | |
a9e9e511 | 1017 | |
47747ea2 | 1018 | char mergeHistName[256] = ""; |
1019 | sprintf(mergeHistName, "%s_cen%i", name.Data(), newCentIndex); | |
a9e9e511 | 1020 | |
47747ea2 | 1021 | // Check if histogram allready exists in current directory. |
1022 | TH1* mergeHist = dynamic_cast<TH1*>( gDirectory->FindObject(mergeHistName) ); | |
1023 | if ( mergeHist ) | |
1024 | return mergeHist; | |
a9e9e511 | 1025 | |
47747ea2 | 1026 | // If not so; get the first hist, clone, and add the others |
1027 | char cname[256] = ""; | |
1028 | sprintf(cname, "%s_cen%i", name.Data(), fromCentIndex); | |
1029 | TH1* hist0 = input.GetHistogram(cname); | |
1030 | sprintf(cname, "%s_cen%i", name.Data(), newCentIndex); | |
1031 | TH1 * histMerged = (TH1*) hist0->Clone(cname); | |
a9e9e511 | 1032 | |
47747ea2 | 1033 | for(int cent=fromCentIndex+1; cent < toCentIndex; ++cent) { |
1034 | sprintf(cname, "%s_cen%i", name.Data(), cent); | |
1035 | TH1* nextHist = input.GetHistogram(cname); | |
1036 | if( ! nextHist ) {Printf("ERROR: Merge of histograms failed"); delete histMerged; return 0x0; } | |
1037 | histMerged->Add(nextHist); | |
1038 | } | |
1039 | return histMerged; | |
1040 | } | |
1041 | ||
1042 | ||
a9e9e511 | 1043 | |
47747ea2 | 1044 | //----------------------------------------------------------------------------- |
1045 | void PPRstyle() | |
1046 | { | |
1047 | ||
1048 | ////////////////////////////////////////////////////////////////////// | |
1049 | // | |
1050 | // ROOT style macro for the TRD TDR | |
1051 | // | |
1052 | ////////////////////////////////////////////////////////////////////// | |
1053 | ||
1054 | gStyle->SetPalette(1); | |
1055 | gStyle->SetCanvasBorderMode(-1); | |
1056 | gStyle->SetCanvasBorderSize(1); | |
1057 | gStyle->SetCanvasColor(10); | |
1058 | ||
1059 | gStyle->SetFrameFillColor(10); | |
1060 | gStyle->SetFrameBorderSize(1); | |
1061 | gStyle->SetFrameBorderMode(-1); | |
1062 | gStyle->SetFrameLineWidth(1.2); | |
1063 | gStyle->SetFrameLineColor(1); | |
1064 | ||
1065 | gStyle->SetHistFillColor(0); | |
1066 | gStyle->SetHistLineWidth(1); | |
1067 | gStyle->SetHistLineColor(1); | |
1068 | ||
1069 | gStyle->SetPadColor(10); | |
1070 | gStyle->SetPadBorderSize(1); | |
1071 | gStyle->SetPadBorderMode(-1); | |
1072 | ||
1073 | gStyle->SetStatColor(10); | |
1074 | gStyle->SetTitleColor(kBlack,"X"); | |
1075 | gStyle->SetTitleColor(kBlack,"Y"); | |
1076 | ||
1077 | gStyle->SetLabelSize(0.04,"X"); | |
1078 | gStyle->SetLabelSize(0.04,"Y"); | |
1079 | gStyle->SetLabelSize(0.04,"Z"); | |
1080 | gStyle->SetTitleSize(0.04,"X"); | |
1081 | gStyle->SetTitleSize(0.04,"Y"); | |
1082 | gStyle->SetTitleSize(0.04,"Z"); | |
1083 | gStyle->SetTitleFont(42,"X"); | |
1084 | gStyle->SetTitleFont(42,"Y"); | |
1085 | gStyle->SetTitleFont(42,"X"); | |
1086 | gStyle->SetLabelFont(42,"X"); | |
1087 | gStyle->SetLabelFont(42,"Y"); | |
1088 | gStyle->SetLabelFont(42,"Z"); | |
1089 | gStyle->SetStatFont(42); | |
1090 | ||
1091 | gStyle->SetTitleOffset(1.0,"X"); | |
1092 | gStyle->SetTitleOffset(1.4,"Y"); | |
1093 | ||
1094 | gStyle->SetFillColor(kWhite); | |
1095 | gStyle->SetTitleFillColor(kWhite); | |
1096 | ||
1097 | gStyle->SetOptDate(0); | |
1098 | gStyle->SetOptTitle(1); | |
1099 | gStyle->SetOptStat(0); | |
1100 | gStyle->SetOptFit(111); | |
1101 | ||
1102 | } | |
1103 | } | |
1104 | ||
1105 | ||
1106 | void MakeRawProduction() | |
1107 | { | |
6aa36480 | 1108 | RawProduction::Output output; |
1109 | ||
1110 | ||
1111 | //TStringToken triggers("kMB kCentral kSemiCentral kPHOSPb", " "); | |
1112 | TStringToken triggers("kMB kPHOSPb", " "); | |
1113 | while(triggers.NextToken()) { | |
a53ea5d3 | 1114 | RawProduction::TriggerBin inBin(triggers); |
6aa36480 | 1115 | RawProduction::Input input("AnalysisResults.root", inBin); |
1116 | TStringToken pids("All Allcore Allwou Disp Disp2 Dispcore Dispwou CPV CPVcore CPV2 Both Bothcore", " "); | |
1117 | //TStringToken pids("Bothcore", " "); | |
1118 | while(pids.NextToken()) { | |
a53ea5d3 | 1119 | RawProduction::TriCenPidBin outBin(-10, pids, inBin.Trigger()); |
6aa36480 | 1120 | RawProduction::MakePi0Fit(input, outBin, output); |
1121 | } | |
1122 | } | |
1123 | output.Write(); | |
1124 | ||
1125 | } | |
1126 | ||
1127 | void MakeRawProductionAll() | |
1128 | { | |
1129 | //gStyle->SetOptStat(0); | |
1130 | gStyle->SetOptFit(1); | |
6aa36480 | 1131 | |
a9e9e511 | 1132 | RawProduction::Output output; |
6aa36480 | 1133 | |
a9e9e511 | 1134 | TStringToken triggers("kMB kCentral kSemiCentral kPHOSPb", " "); |
1135 | //TStringToken triggers("kCentral ", " "); | |
6aa36480 | 1136 | while(triggers.NextToken()) { |
a9e9e511 | 1137 | RawProduction::TriggerBin triggerBin(triggers); |
1138 | RawProduction::Input input("AnalysisResults.root", triggerBin); | |
1139 | ||
1140 | RawProduction::MakePi0Fit(input, triggerBin, output); | |
1141 | ||
6aa36480 | 1142 | //TStringToken pids("All Allcore Allwou Disp Disp2 Dispcore Dispwou CPV CPVcore CPV2 Both Bothcore", " "); |
1143 | //TStringToken pids("All Allcore Allwou Disp Disp2 Dispcore Dispwou", " "); | |
1144 | TStringToken pids("All", " "); | |
1145 | while(pids.NextToken()) { | |
1146 | for(int cent = -1; cent > -7; cent--) { | |
1147 | if(triggers.EqualTo("kCentral") && cent != -1) continue; | |
1148 | if(triggers.EqualTo("kSemiCentral") && !(-1 > cent && cent > -6 )) continue; | |
1149 | ||
a9e9e511 | 1150 | RawProduction::TriCenPidBin tcpBin(cent, pids, triggerBin.Trigger()); |
1151 | RawProduction::MakePi0FitTCP(input, tcpBin, output); | |
6aa36480 | 1152 | } |
a9e9e511 | 1153 | |
6aa36480 | 1154 | if( triggers.EqualTo("kCentral") || triggers.EqualTo("kSemiCentral") ) continue; |
a9e9e511 | 1155 | RawProduction::TriCenPidBin tcpBinc(-10, pids, triggerBin.Trigger()); |
1156 | RawProduction::MakePi0FitTCP(input, tcpBinc, output); | |
6aa36480 | 1157 | } |
20fb7b55 | 1158 | } |
1159 | output.Write(); | |
47747ea2 | 1160 | } |