PWGGA/PHOSTasks/PHOS_PbPb/macros/production/MakeRawProduction.C

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