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