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