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