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