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