]>
Commit | Line | Data |
---|---|---|
b8bea077 | 1 | /* $Id$ */ |
2 | ||
7195e4ab | 3 | #include "TFile.h" |
4 | #include "TF1.h" | |
5 | #include "TH1F.h" | |
6 | #include "TH2F.h" | |
7 | #include "TSystem.h" | |
8 | #include "TStyle.h" | |
9 | #include "TCanvas.h" | |
10 | #include "TMath.h" | |
11 | ||
12 | void PPRstyle(); | |
13 | ||
14 | Double_t PeakPosition(Double_t pt); | |
15 | Double_t PeakWidth(Double_t pt); | |
16 | Double_t CB(Double_t * x, Double_t * par); | |
17 | Double_t CB2(Double_t * x, Double_t * par); | |
18 | Double_t CBs(Double_t * x, Double_t * par); | |
19 | Double_t BG1(Double_t * x, Double_t * par); | |
20 | Double_t BG2(Double_t * x, Double_t * par); | |
21 | ||
22 | ||
2eb453ea | 23 | // const Int_t nPadX = 3, nPadY = 2; |
24 | // const Int_t nPtBins=6 ; | |
25 | // const Double_t ptBinEdges[21]={1., 2., 3., 4., 5., 7., 10.} ; | |
26 | ||
5fb9fcd6 | 27 | const Int_t nPadX = 5, nPadY = 4; |
7195e4ab | 28 | Int_t nPtBins=0; |
29 | Double_t ptBinEdges[1000] = {0}; | |
30 | double GetPtBin(int bin){ | |
31 | if( bin==0 ) | |
32 | return 1.; | |
33 | ||
34 | // return GetPtBin(bin-1) * 1.1; | |
5fb9fcd6 | 35 | |
7195e4ab | 36 | // if ( bin % 2 ) |
37 | // return GetPtBin(bin-1) + 0.4; | |
38 | // else | |
39 | // return GetPtBin(bin-1) + 0.2; | |
5fb9fcd6 | 40 | |
7195e4ab | 41 | double previousBin = GetPtBin(bin-1); |
42 | double linInc = 0.2; | |
43 | double threshold = 5.; | |
44 | double logFact = 1 + linInc/threshold; | |
45 | if ( previousBin < threshold ) // linear | |
46 | return previousBin + linInc; | |
47 | else { // logarithmic | |
48 | return previousBin * logFact; | |
49 | } | |
50 | } | |
51 | void MakePtBins() { | |
52 | int bin = -1; | |
53 | do { | |
54 | ++bin; | |
55 | ptBinEdges[bin] = GetPtBin(bin); | |
56 | } while(ptBinEdges[bin] < 40.); | |
57 | nPtBins = bin -2; | |
58 | ||
5fb9fcd6 | 59 | printf("Making Pt Bins:\n"); |
7195e4ab | 60 | for(int b=0; b < nPtBins+1; ++b) |
61 | printf("%.1f, ", ptBinEdges[b]); | |
62 | printf("\n N. Bins: %d\n", nPtBins); | |
5fb9fcd6 | 63 | |
7195e4ab | 64 | |
65 | // for(int bin = 0; bin <= nPtBins; ++bin){ | |
66 | // ptBinEdges[bin] = GetPtBin(bin); | |
67 | // printf("%.1f, ", ptBinEdges[bin]); | |
68 | // } | |
69 | // printf("\n"); | |
70 | } | |
71 | ||
2eb453ea | 72 | |
73 | const int kNCentralityBins = 3; | |
efb6c264 | 74 | |
8e6cc59c | 75 | const Double_t kMean=0.136 ; //Approximate peak position to facilitate error estimate |
76 | ||
5fb9fcd6 | 77 | const char format[] = ".pdf"; // say if you want .pdf |
efb6c264 | 78 | |
2eb453ea | 79 | TH2F* FindPi0(TList *histoList, bool mix, const char* pid, int centrality); |
80 | ||
8e6cc59c | 81 | //----------------------------------------------------------------------------- |
efb6c264 | 82 | void MakeMmixPi0(const TString filename, |
83 | const TString listPath = "PHOSPi0Flow/PHOSPi0FlowCoutput1", // lego train | |
5fb9fcd6 | 84 | const Int_t centrality=0, |
efb6c264 | 85 | const char* pid="CPV", |
5fb9fcd6 | 86 | const char* trigger="kCentral", |
efb6c264 | 87 | const char* saveToDir="") |
8e6cc59c | 88 | { |
7195e4ab | 89 | MakePtBins(); |
efb6c264 | 90 | Printf("\nMakeMmixPi0(%s, %s, %i, %s, %s)", filename.Data(), listPath.Data(), centrality, pid, saveToDir); |
91 | ||
2eb453ea | 92 | if( TString(saveToDir).Length() ) |
93 | gSystem->mkdir(saveToDir, true); | |
94 | ||
8e6cc59c | 95 | //Fit Real/Mixed ratio, normalize Mixed and subtract it from Real |
96 | ||
efb6c264 | 97 | TFile * file = new TFile(filename) ; |
98 | TList *histoList = (TList*)file->Get(listPath); | |
8b4ce769 | 99 | |
8e6cc59c | 100 | char key[125] ; |
101 | ||
102 | TH1F * hev = (TH1F*)histoList->FindObject("hTotSelEvents") ; | |
103 | TH2F * hCentrality = (TH2F*)histoList->FindObject("hCenPHOSCells") ; | |
efb6c264 | 104 | TH1D * hCentralityX = hCentrality->ProjectionX(); |
5fb9fcd6 | 105 | |
106 | printf("TotSelEvents (4): %.0f \n",hev->GetBinContent(4)) ; | |
107 | printf("Centrality: %.0f \n",hCentralityX->Integral()) ; | |
108 | ||
2eb453ea | 109 | TH2F *hPi0 = FindPi0(histoList, false, pid, centrality); |
110 | TH2F *hPi0Mix = FindPi0(histoList, true, pid, centrality); | |
efb6c264 | 111 | if( !hPi0 || !hPi0Mix || hPi0->GetEntries() < 10000) { |
ebfd4ce0 | 112 | Printf(Form("no histogram(0x%p, 0x%p) or to low number of entries(%.1f) in hPi0, skipping this cent", hPi0, hPi0Mix, hPi0->GetEntries())); |
efb6c264 | 113 | file->Close(); |
114 | delete file; | |
8e6cc59c | 115 | return; |
116 | } | |
117 | ||
5fb9fcd6 | 118 | |
119 | TFile* outFile = new TFile(Form("%sPi0_FitResult.root", saveToDir, centrality),"update"); | |
8e6cc59c | 120 | |
121 | PPRstyle(); | |
122 | gStyle->SetPadLeftMargin(0.14); | |
123 | gStyle->SetPadRightMargin(0.01); | |
bd741163 | 124 | //gStyle->SetPadTopMargin(0.01); |
8e6cc59c | 125 | gStyle->SetPadBottomMargin(0.08); |
126 | ||
5fb9fcd6 | 127 | //Fit real only |
8e6cc59c | 128 | //Linear Bg |
129 | char kkey[55]; | |
5fb9fcd6 | 130 | sprintf(kkey, Form("%s_cent%d_%s",pid, centrality, trigger)) ; |
8e6cc59c | 131 | char key2[155]; |
132 | sprintf(key,"Mix%s",kkey) ; | |
133 | sprintf(key2,"%s_mr1",key) ; | |
25a78f54 | 134 | TH1D * mr1 = new TH1D(key2,"Mass",nPtBins,ptBinEdges) ; |
8e6cc59c | 135 | sprintf(key2,"%s_sr1",key) ; |
25a78f54 | 136 | TH1D * sr1 = new TH1D(key2,"Width",nPtBins,ptBinEdges) ; |
8e6cc59c | 137 | sprintf(key2,"%s_ar1",key) ; |
25a78f54 | 138 | TH1D * ar1 = new TH1D(key2,"a",nPtBins,ptBinEdges) ; |
8e6cc59c | 139 | sprintf(key2,"%s_br1",key) ; |
25a78f54 | 140 | TH1D * br1 = new TH1D(key2,"a",nPtBins,ptBinEdges) ; |
8e6cc59c | 141 | sprintf(key2,"%s_yr1",key) ; |
25a78f54 | 142 | TH1D * nr1 = new TH1D(key2,"Raw yield",nPtBins,ptBinEdges) ; |
8e6cc59c | 143 | sprintf(key2,"%s_yr1int",key) ; |
25a78f54 | 144 | TH1D * nr1int = new TH1D(key2,"Raw yield, integrated",nPtBins,ptBinEdges) ; |
8e6cc59c | 145 | |
146 | //Quadratic Bg | |
147 | sprintf(key2,"%s_mr2",key) ; | |
25a78f54 | 148 | TH1D * mr2 = new TH1D(key2,"Mass",nPtBins,ptBinEdges) ; |
8e6cc59c | 149 | sprintf(key2,"%s_sr2",key) ; |
25a78f54 | 150 | TH1D * sr2 = new TH1D(key2,"Width",nPtBins,ptBinEdges) ; |
8e6cc59c | 151 | sprintf(key2,"%s_ar2",key) ; |
25a78f54 | 152 | TH1D * ar2 = new TH1D(key2,"a",nPtBins,ptBinEdges) ; |
8e6cc59c | 153 | sprintf(key2,"%s_br2",key) ; |
25a78f54 | 154 | TH1D * br2 = new TH1D(key2,"a",nPtBins,ptBinEdges) ; |
8e6cc59c | 155 | sprintf(key2,"%s_cr2",key) ; |
25a78f54 | 156 | TH1D * cr2 = new TH1D(key2,"a",nPtBins,ptBinEdges) ; |
8e6cc59c | 157 | sprintf(key2,"%s_yr2",key) ; |
25a78f54 | 158 | TH1D * nr2 = new TH1D(key2,"Raw yield",nPtBins,ptBinEdges) ; |
8e6cc59c | 159 | sprintf(key2,"%s_yr2int",key) ; |
25a78f54 | 160 | TH1D * nr2int = new TH1D(key2,"Raw yield, integrated",nPtBins,ptBinEdges) ; |
8e6cc59c | 161 | |
5fb9fcd6 | 162 | // Functions: |
163 | TF1 * funcRatioFit1 = new TF1("funcRatioFit1",CB,0.,1.,6) ; | |
164 | funcRatioFit1->SetParName(0,"A") ; | |
165 | funcRatioFit1->SetParName(1,"m_{0}") ; | |
166 | funcRatioFit1->SetParName(2,"#sigma") ; | |
167 | funcRatioFit1->SetParName(3,"a_{0}") ; | |
168 | funcRatioFit1->SetParName(4,"a_{1}") ; | |
169 | funcRatioFit1->SetParName(5,"a_{2}") ; | |
170 | funcRatioFit1->SetLineWidth(2) ; | |
171 | funcRatioFit1->SetLineColor(2) ; | |
8e6cc59c | 172 | TF1 * fgs = new TF1("gs",CBs,0.,1.,4) ; |
173 | fgs->SetParName(0,"A") ; | |
174 | fgs->SetParName(1,"m_{0}") ; | |
175 | fgs->SetParName(2,"#sigma") ; | |
176 | fgs->SetParName(3,"B") ; | |
177 | fgs->SetLineColor(2) ; | |
178 | fgs->SetLineWidth(1) ; | |
5fb9fcd6 | 179 | TF1 * funcRatioFit2 = new TF1("funcRatioFit2",CB2,0.,1.,7) ; |
180 | funcRatioFit2->SetParName(0,"A") ; | |
181 | funcRatioFit2->SetParName(1,"m_{0}") ; | |
182 | funcRatioFit2->SetParName(2,"#sigma") ; | |
183 | funcRatioFit2->SetParName(3,"a_{0}") ; | |
184 | funcRatioFit2->SetParName(4,"a_{1}") ; | |
185 | funcRatioFit2->SetParName(5,"a_{2}") ; | |
186 | funcRatioFit2->SetParName(6,"a_{3}") ; | |
187 | funcRatioFit2->SetLineWidth(2) ; | |
188 | funcRatioFit2->SetLineColor(4) ; | |
189 | funcRatioFit2->SetLineStyle(2) ; | |
8e6cc59c | 190 | TF1 * fbg1 = new TF1("bg1",BG1,0.,1.,3) ; |
191 | TF1 * fbg2 = new TF1("bg2",BG2,0.,1.,4) ; | |
192 | ||
8e6cc59c | 193 | |
5fb9fcd6 | 194 | // Canvases |
25a78f54 | 195 | TCanvas * rawCanvas = new TCanvas("rawCanvas","rawCanvas",10,10,1200,800); |
196 | rawCanvas->Divide(nPadX, nPadY); | |
5fb9fcd6 | 197 | TCanvas * ratioCanvas = new TCanvas("mggFit1","mggFit1",10,10,1200,800) ; |
198 | ratioCanvas->Divide(nPadX,nPadY) ; | |
199 | TCanvas * sbsCanvas = new TCanvas("mggFit1_Signal","mggFitCB",10,10,1200,800) ; | |
200 | sbsCanvas->Divide(nPadX,nPadY) ; | |
201 | ||
25a78f54 | 202 | |
25a78f54 | 203 | for(Int_t ptBin=1; ptBin<=nPtBins; ptBin++){ |
5fb9fcd6 | 204 | ratioCanvas->cd(ptBin) ; |
205 | TAxis * pta=hPi0->GetYaxis() ; | |
25a78f54 | 206 | Int_t imin=pta->FindBin(ptBinEdges[ptBin-1]+0.0001); |
207 | Int_t imax=pta->FindBin(ptBinEdges[ptBin]-0.0001) ; | |
208 | Double_t pt=(ptBinEdges[ptBin]+ptBinEdges[ptBin-1])/2. ; | |
2eb453ea | 209 | |
210 | TH1D * hPi0Proj = hPi0->ProjectionX(Form("hPi0_ptBin%d",ptBin), imin, imax) ; | |
211 | TH1D * hPi0ProjMix= hPi0Mix->ProjectionX(Form("hPi0Mix_ptBin%d",ptBin), imin, imax) ; | |
212 | ||
213 | hPi0Proj->SetTitle(Form("M_{#gamma#gamma}, PID=%s, %.1f<p_{T}<%.1f GeV/c",pid,ptBinEdges[ptBin-1],ptBinEdges[ptBin])); | |
214 | hPi0Proj->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})"); | |
215 | hPi0ProjMix->SetTitle(Form("M_{#gamma#gamma}^{Mix}, PID=%s, %.1f<p_{T}<%.1f GeV/c",pid,ptBinEdges[ptBin-1],ptBinEdges[ptBin])); | |
216 | hPi0ProjMix->SetXTitle("M_{#gamma#gamma}^{Mix} (GeV/c^{2})"); | |
ebfd4ce0 | 217 | |
ebfd4ce0 | 218 | Printf("\n\t%.1f<pt<%.1f GeV/c, entries: %.0f",ptBinEdges[ptBin-1],ptBinEdges[ptBin], hPi0Proj->GetEntries()); |
219 | if( hPi0Proj->GetEntries() < 100. ) { | |
5fb23089 | 220 | Printf("skipping this bin as n. entries is to low"); |
5fb23089 | 221 | continue; |
222 | } | |
5fb23089 | 223 | |
ebfd4ce0 | 224 | hPi0Proj ->Rebin(4) ; |
225 | hPi0ProjMix->Rebin(4) ; | |
ebfd4ce0 | 226 | hPi0Proj->SetName(Form("%s_rebin", hPi0Proj->GetName())); |
227 | hPi0ProjMix->SetName(Form("%s_rebin", hPi0ProjMix->GetName())); | |
228 | ||
5fb9fcd6 | 229 | // Error Fix |
2eb453ea | 230 | for(Int_t ib=1; ib<=hPi0Proj->GetNbinsX();ib++){if(hPi0Proj ->GetBinContent(ib)==0)hPi0Proj ->SetBinError(ib,1.);} |
231 | for(Int_t ib=1; ib<=hPi0Proj->GetNbinsX();ib++){if(hPi0ProjMix->GetBinContent(ib)==0)hPi0ProjMix->SetBinError(ib,1.);} | |
5fb9fcd6 | 232 | |
233 | // Signal-Mix Ratio | |
234 | TH1D * hPi0Ratio = (TH1D*)hPi0Proj->Clone( Form("hPi0Ratio_ptBin%d",ptBin) ) ; | |
235 | hPi0Ratio->SetTitle(Form("#frac{M_{#gamma#gamma}}{M_{#gamma#gamma}^{Mix}}, %.1f<p_{T}<%.1f GeV/c", ptBinEdges[ptBin-1], ptBinEdges[ptBin])); | |
236 | hPi0Ratio->Divide(hPi0ProjMix) ; | |
237 | hPi0Ratio->SetMarkerStyle(20) ; | |
238 | hPi0Ratio->SetMarkerSize(0.7) ; | |
239 | ||
240 | funcRatioFit1->SetParameters(0.001,0.136,0.0055,0.0002,-0.002,0.0) ; | |
241 | funcRatioFit1->SetParLimits(0,0.000,1.000) ; | |
242 | funcRatioFit1->SetParLimits(1,0.120,0.145) ; | |
243 | funcRatioFit1->SetParLimits(2,0.004,0.012) ; | |
8e6cc59c | 244 | |
245 | Double_t rangeMin=0.05 ; | |
246 | Double_t rangeMax=0.3 ; | |
247 | if(centrality==0) rangeMax=0.4 ; | |
25a78f54 | 248 | if(ptBin==1){ |
8e6cc59c | 249 | rangeMin=0.06 ; |
250 | rangeMax=0.25 ; | |
251 | } | |
5fb9fcd6 | 252 | int error = hPi0Ratio->Fit(funcRatioFit1,"Q" ,"",rangeMin,rangeMax) ; |
2eb453ea | 253 | if( error ) { |
254 | Printf("fit (ratio) error: %d ", error); | |
255 | continue; | |
256 | } | |
5fb9fcd6 | 257 | error = hPi0Ratio->Fit(funcRatioFit1,"MQ","",rangeMin,rangeMax) ; |
2eb453ea | 258 | if( error % 4000) { |
259 | Printf("fit (ratio more) error: %d ", error); | |
260 | continue; | |
261 | } | |
262 | ||
8e6cc59c | 263 | |
5fb9fcd6 | 264 | ar1->SetBinContent(ptBin,funcRatioFit1->GetParameter(3)) ; |
265 | ar1->SetBinError (ptBin,funcRatioFit1->GetParError(3)) ; | |
266 | br1->SetBinContent(ptBin,funcRatioFit1->GetParameter(4)) ; | |
267 | br1->SetBinError (ptBin,funcRatioFit1->GetParError(4)) ; | |
8e6cc59c | 268 | |
5fb9fcd6 | 269 | funcRatioFit2->SetParameters(funcRatioFit1->GetParameters()) ; |
270 | funcRatioFit2->SetParLimits(0,0.000,1.000) ; | |
271 | funcRatioFit2->SetParLimits(1,0.110,0.145) ; | |
272 | funcRatioFit2->SetParLimits(2,0.004,0.012) ; | |
8e6cc59c | 273 | |
5fb9fcd6 | 274 | error = hPi0Ratio->Fit(funcRatioFit2,"+NQ","",rangeMin,rangeMax) ; |
2eb453ea | 275 | if( error ) { |
5fb9fcd6 | 276 | Printf("fit (funcRatioFit2) error: %d ", error); |
2eb453ea | 277 | continue; |
278 | } | |
5fb9fcd6 | 279 | error =hPi0Ratio->Fit(funcRatioFit2,"+MQ","",rangeMin,rangeMax) ; |
2eb453ea | 280 | if( error % 4000) { |
5fb9fcd6 | 281 | Printf("fit (funcRatioFit2 more) error: %d ", error); |
2eb453ea | 282 | continue; |
283 | } | |
8e6cc59c | 284 | |
5fb9fcd6 | 285 | ar2->SetBinContent(ptBin,funcRatioFit2->GetParameter(3)) ; |
286 | ar2->SetBinError (ptBin,funcRatioFit2->GetParError(3)) ; | |
287 | br2->SetBinContent(ptBin,funcRatioFit2->GetParameter(4)) ; | |
288 | br2->SetBinError (ptBin,funcRatioFit2->GetParError(4)) ; | |
289 | cr2->SetBinContent(ptBin,funcRatioFit2->GetParameter(5)) ; | |
290 | cr2->SetBinError (ptBin,funcRatioFit2->GetParError(5)) ; | |
291 | hPi0Ratio->GetXaxis()->SetRangeUser(0.05,0.30) ; | |
292 | hPi0Ratio->DrawCopy() ; | |
293 | hPi0Ratio->Write(0,TObject::kOverwrite) ; | |
294 | ratioCanvas->Update() ; | |
295 | ||
296 | fbg1->SetParameters(funcRatioFit1->GetParameter(3),funcRatioFit1->GetParameter(4),funcRatioFit1->GetParameter(5)); | |
297 | fbg2->SetParameters(funcRatioFit2->GetParameter(3),funcRatioFit2->GetParameter(4),funcRatioFit2->GetParameter(5), | |
298 | funcRatioFit2->GetParameter(6)); | |
8e6cc59c | 299 | |
300 | Double_t intRangeMin = PeakPosition(pt)-3.*PeakWidth(pt) ; | |
301 | Double_t intRangeMax = PeakPosition(pt)+3.*PeakWidth(pt) ; | |
2eb453ea | 302 | Int_t intBinMin = hPi0Proj->GetXaxis()->FindBin(intRangeMin) ; |
303 | Int_t intBinMax = hPi0Proj->GetXaxis()->FindBin(intRangeMax) ; | |
5fb9fcd6 | 304 | Double_t errStat = hPi0ProjMix->Integral(intBinMin,intBinMax); |
305 | ||
bd741163 | 306 | rawCanvas->cd(ptBin); |
5fb9fcd6 | 307 | TH1D * hPi0MixScaled = (TH1D*)hPi0ProjMix ->Clone(Form("%s_clone2", hPi0Proj->GetName())) ; |
308 | hPi0MixScaled->Scale(fbg1->Eval(0.1349)); | |
2eb453ea | 309 | hPi0Proj->SetLineColor(kBlack); |
310 | hPi0Proj->SetAxisRange(0.0, 0.3); | |
311 | hPi0Proj->DrawCopy(); | |
312 | hPi0Proj->Write(0,TObject::kOverwrite) ; | |
5fb9fcd6 | 313 | hPi0MixScaled->SetLineColor(kRed); |
314 | hPi0MixScaled->SetTitle(Form("M_{#gamma#gamma}^{Mix,scaled}, PID=%s, %.1f<p_{T}<%.1f GeV/c",pid,ptBinEdges[ptBin-1],ptBinEdges[ptBin])); | |
315 | hPi0MixScaled->DrawCopy("same"); | |
bd741163 | 316 | rawCanvas->Update(); |
5fb9fcd6 | 317 | hPi0MixScaled->Write(0,TObject::kOverwrite) ; |
318 | ||
319 | // Linear background subtraction | |
320 | TH1D * hPi0MixScaledPol1 = (TH1D*)hPi0ProjMix->Clone(Form("hPi0MixScaledPol1_ptBin%d", ptBin)) ; | |
321 | hPi0MixScaledPol1 ->Multiply(fbg1) ; | |
322 | TH1D * hPi0BSPol1 = (TH1D*)hPi0Proj->Clone(Form("hPi0BSPol1_ptBin%d", ptBin)) ; | |
323 | hPi0BSPol1->Add(hPi0MixScaledPol1 ,-1.) ; | |
324 | ||
325 | // Quadratic background | |
326 | TH1D * hPi0MixScaledPol2 = (TH1D*)hPi0ProjMix->Clone(Form("hPi0MixScaledPol2_ptBin%d", ptBin)) ; | |
327 | hPi0MixScaledPol2->Multiply(fbg2) ; | |
328 | TH1D * hPi0BSPol2 = (TH1D*)hPi0Proj ->Clone(Form("hPi0BSPol2_ptBin%d", ptBin)) ; | |
329 | hPi0BSPol2 ->Add(hPi0MixScaledPol2,-1.) ; | |
330 | ||
331 | sbsCanvas->cd(ptBin) ; | |
332 | ||
333 | Int_t binPi0 = hPi0BSPol1->FindBin(funcRatioFit1->GetParameter(1)); | |
334 | Int_t nWidPi0 = 2 * (Int_t) (funcRatioFit1->GetParameter(2)/hPi0BSPol1->GetBinWidth(1)); | |
335 | fgs->SetParameters(hPi0BSPol1->Integral(binPi0-nWidPi0,binPi0+nWidPi0)/5., funcRatioFit1->GetParameter(1), funcRatioFit1->GetParameter(2)) ; | |
336 | fgs->SetParLimits(0,0.,1.e+5) ; | |
033398eb | 337 | fgs->SetParLimits(1,0.110,0.145) ; |
8e6cc59c | 338 | fgs->SetParLimits(2,0.004,0.02) ; |
5fb9fcd6 | 339 | error = hPi0BSPol1->Fit(fgs,"Q","",rangeMin,rangeMax) ; |
2eb453ea | 340 | if( error ) { |
5fb9fcd6 | 341 | Printf("fit (hPi0BSPol1 fgs) error: %d ", error); |
2eb453ea | 342 | continue; |
343 | } | |
5fb9fcd6 | 344 | hPi0BSPol1->SetMaximum(hPi0BSPol2->GetMaximum()*1.5) ; |
345 | hPi0BSPol1->SetMinimum(hPi0BSPol2->GetMinimum()*1.1) ; | |
346 | hPi0BSPol1->SetMarkerStyle(20) ; | |
347 | hPi0BSPol1->SetMarkerSize(0.7) ; | |
25a78f54 | 348 | mr1->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
349 | mr1->SetBinError (ptBin,fgs->GetParError(1) ) ; | |
350 | sr1->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; | |
351 | sr1->SetBinError (ptBin,fgs->GetParError(2) ) ; | |
8e6cc59c | 352 | |
5fb9fcd6 | 353 | Double_t y=fgs->GetParameter(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ; |
25a78f54 | 354 | nr1->SetBinContent(ptBin,y) ; |
5fb9fcd6 | 355 | Double_t ey=fgs->GetParError(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ; |
25a78f54 | 356 | nr1->SetBinError(ptBin,ey) ; |
8e6cc59c | 357 | |
5fb9fcd6 | 358 | Double_t npiInt = hPi0BSPol1->Integral(intBinMin,intBinMax) ; |
8e6cc59c | 359 | Double_t norm = fbg1->GetParameter(0) ; |
360 | Double_t normErr= fbg1->GetParError(0) ; | |
361 | if(npiInt>0.){ | |
25a78f54 | 362 | nr1int->SetBinContent(ptBin,npiInt) ; |
363 | nr1int->SetBinError(ptBin,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ; | |
8e6cc59c | 364 | } |
5fb9fcd6 | 365 | hPi0BSPol2->GetXaxis()->SetRangeUser(0.05,0.3) ; |
366 | hPi0BSPol2->SetMaximum(hPi0BSPol2->GetMaximum()*1.5) ; | |
367 | hPi0BSPol2->SetMinimum(hPi0BSPol2->GetMinimum()*1.1) ; | |
368 | hPi0BSPol2->SetMarkerStyle(20) ; | |
369 | hPi0BSPol2->SetMarkerSize(0.7) ; | |
370 | hPi0BSPol2->Fit(fgs,"Q","",rangeMin,rangeMax) ; | |
2eb453ea | 371 | if( error ) { |
5fb9fcd6 | 372 | Printf("fit (hPi0BSPol2 fgs) error: %d ", error); |
2eb453ea | 373 | continue; |
374 | } | |
25a78f54 | 375 | mr2->SetBinContent(ptBin,fgs->GetParameter(1)) ; |
376 | mr2->SetBinError (ptBin,fgs->GetParError(1)) ; | |
377 | sr2->SetBinContent(ptBin,TMath::Abs(fgs->GetParameter(2))) ; | |
378 | sr2->SetBinError (ptBin,fgs->GetParError(2)) ; | |
5fb9fcd6 | 379 | y=fgs->GetParameter(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ; |
25a78f54 | 380 | nr2->SetBinContent(ptBin,y) ; |
5fb9fcd6 | 381 | ey= fgs->GetParError(0)/hPi0BSPol1->GetXaxis()->GetBinWidth(1) ; |
25a78f54 | 382 | nr2->SetBinError(ptBin,ey) ; |
5fb9fcd6 | 383 | Double_t npiInt2=hPi0BSPol2->Integral(intBinMin,intBinMax) ; |
8e6cc59c | 384 | norm=fbg2->GetParameter(0) ; |
385 | normErr=fbg2->GetParError(0) ; | |
5fb9fcd6 | 386 | if(npiInt2>0.){ |
387 | nr2int->SetBinContent(ptBin,npiInt2) ; | |
388 | nr2int->SetBinError(ptBin,TMath::Sqrt(npiInt2 + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ; | |
389 | } | |
390 | hPi0BSPol2->SetTitle(Form("M_{#gamma#gamma}^{BS_{2}}, PID=%s, %.1f<p_{T}<%.1f GeV/c",pid,ptBinEdges[ptBin-1],ptBinEdges[ptBin])); | |
391 | hPi0BSPol2->DrawCopy() ; | |
392 | sbsCanvas->Update() ; | |
393 | hPi0BSPol2->Write(0,TObject::kOverwrite) ; | |
8e6cc59c | 394 | } |
5fb9fcd6 | 395 | |
efb6c264 | 396 | char name[55] ; |
397 | sprintf(name,"%sPi0_ratio_%s%s", saveToDir, kkey, format) ; | |
5fb9fcd6 | 398 | ratioCanvas->Print(name) ; |
efb6c264 | 399 | sprintf(name,"%sPi0_signal_%s%s", saveToDir, kkey, format) ; |
5fb9fcd6 | 400 | sbsCanvas->Print(name) ; |
efb6c264 | 401 | sprintf(name,"%sPi0_raw_%s%s", saveToDir, kkey, format) ; |
402 | rawCanvas->Print(name); | |
8e6cc59c | 403 | |
404 | //Normalize by the number of events | |
efb6c264 | 405 | Int_t cMin=0, cMax=100; |
8e6cc59c | 406 | if (centrality == 0) { |
407 | cMin=1; | |
408 | cMax=10; | |
409 | } | |
410 | else if (centrality == 1) { | |
efb6c264 | 411 | cMin=11; |
8e6cc59c | 412 | cMax=40; |
413 | } | |
414 | else if (centrality == 2) { | |
415 | cMin=41; | |
8e6cc59c | 416 | cMax=80; |
417 | } | |
2eb453ea | 418 | else if (centrality == -1) { |
419 | cMin=1; | |
420 | cMax=80; | |
421 | } | |
efb6c264 | 422 | else { |
423 | Printf("ERROR: Centrality: %d not defined !!! ERROR", centrality); | |
424 | return; | |
425 | } | |
5fb9fcd6 | 426 | |
efb6c264 | 427 | Double_t nevents = hCentralityX->Integral(cMin,cMax); |
2eb453ea | 428 | if ( nevents > 0.9 ) { |
429 | nr1 ->Scale(1./nevents) ; | |
430 | nr1int->Scale(1./nevents) ; | |
431 | nr2 ->Scale(1./nevents) ; | |
432 | nr2int->Scale(1./nevents) ; | |
433 | } else { | |
434 | Printf("WARNING: non positive nEvents in centrality range, cMin:%d, cMax:%d, nEvents:%f", cMin, cMax, nevents ); | |
5fb9fcd6 | 435 | |
2eb453ea | 436 | } |
2eb453ea | 437 | |
8e6cc59c | 438 | mr1->Write(0,TObject::kOverwrite) ; |
439 | sr1->Write(0,TObject::kOverwrite) ; | |
440 | ar1->Write(0,TObject::kOverwrite) ; | |
441 | br1->Write(0,TObject::kOverwrite) ; | |
442 | nr1->Write(0,TObject::kOverwrite) ; | |
443 | nr1int->Write(0,TObject::kOverwrite) ; | |
444 | ar2->Write(0,TObject::kOverwrite) ; | |
445 | br2->Write(0,TObject::kOverwrite) ; | |
446 | cr2->Write(0,TObject::kOverwrite) ; | |
447 | mr2->Write(0,TObject::kOverwrite) ; | |
448 | sr2->Write(0,TObject::kOverwrite) ; | |
449 | nr2->Write(0,TObject::kOverwrite) ; | |
450 | nr2int->Write(0,TObject::kOverwrite) ; | |
8e6cc59c | 451 | |
efb6c264 | 452 | |
453 | outFile->Close(); | |
454 | delete outFile; | |
455 | ||
456 | file->Close(); | |
457 | delete file; | |
2eb453ea | 458 | |
5fb9fcd6 | 459 | delete ratioCanvas; |
460 | delete sbsCanvas; | |
461 | delete rawCanvas; | |
8e6cc59c | 462 | } |
463 | ||
464 | //----------------------------------------------------------------------------- | |
465 | Double_t PeakPosition(Double_t pt){ | |
466 | //Fit to the measured peak position | |
467 | return 4.99292e-003*exp(-pt*9.32300e-001)+1.34944e-001 ; | |
468 | } | |
469 | //----------------------------------------------------------------------------- | |
470 | Double_t PeakWidth(Double_t pt){ | |
471 | //fit to the measured peak width | |
472 | Double_t a=0.0068 ; | |
473 | Double_t b=0.0025 ; | |
474 | Double_t c=0.000319 ; | |
475 | return TMath::Sqrt(a*a+b*b/pt/pt+c*c*pt*pt) ; | |
476 | } | |
5fb9fcd6 | 477 | |
8e6cc59c | 478 | //----------------------------------------------------------------------------- |
479 | Double_t CB(Double_t * x, Double_t * par){ | |
480 | //Parameterization of Real/Mixed ratio | |
481 | Double_t m=par[1] ; | |
482 | Double_t s=par[2] ; | |
483 | Double_t dx=(x[0]-m)/s ; | |
484 | return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean) ; | |
485 | } | |
486 | ||
487 | //----------------------------------------------------------------------------- | |
488 | Double_t CB2(Double_t * x, Double_t * par){ | |
489 | //Another parameterization of Real/Mixed ratio7TeV/README | |
490 | Double_t m=par[1] ; | |
491 | Double_t s=par[2] ; | |
492 | Double_t dx=(x[0]-m)/s ; | |
493 | return par[0]*exp(-dx*dx/2.)+par[3]+par[4]*(x[0]-kMean)+par[5]*(x[0]-kMean)*(x[0]-kMean) ; | |
494 | } | |
495 | //----------------------------------------------------------------------------- | |
496 | Double_t CBs(Double_t * x, Double_t * par){ | |
497 | //Parameterizatin of signal | |
498 | Double_t m=par[1] ; | |
499 | Double_t s=par[2] ; | |
500 | Double_t dx=(x[0]-m)/s ; | |
501 | return par[0]*exp(-dx*dx/2.)/TMath::Sqrt(TMath::TwoPi())/s+par[3] ; | |
502 | } | |
503 | //----------------------------------------------------------------------------- | |
504 | Double_t BG1(Double_t * x, Double_t * par){ | |
505 | //Normalizatino of Mixed | |
506 | return par[0]+par[1]*(x[0]-kMean) ; | |
507 | } | |
508 | //----------------------------------------------------------------------------- | |
509 | Double_t BG2(Double_t * x, Double_t * par){ | |
510 | //Another normalization of Mixed | |
511 | return par[0]+par[1]*(x[0]-kMean)+par[2]*(x[0]-kMean)*(x[0]-kMean) ; | |
512 | } | |
513 | ||
514 | ||
515 | //----------------------------------------------------------------------------- | |
ebfd4ce0 | 516 | void PPRstyle() |
8e6cc59c | 517 | { |
518 | ||
519 | ////////////////////////////////////////////////////////////////////// | |
520 | // | |
521 | // ROOT style macro for the TRD TDR | |
522 | // | |
523 | ////////////////////////////////////////////////////////////////////// | |
524 | ||
525 | gStyle->SetPalette(1); | |
526 | gStyle->SetCanvasBorderMode(-1); | |
527 | gStyle->SetCanvasBorderSize(1); | |
528 | gStyle->SetCanvasColor(10); | |
529 | ||
530 | gStyle->SetFrameFillColor(10); | |
531 | gStyle->SetFrameBorderSize(1); | |
532 | gStyle->SetFrameBorderMode(-1); | |
533 | gStyle->SetFrameLineWidth(1.2); | |
534 | gStyle->SetFrameLineColor(1); | |
535 | ||
536 | gStyle->SetHistFillColor(0); | |
537 | gStyle->SetHistLineWidth(1); | |
538 | gStyle->SetHistLineColor(1); | |
539 | ||
540 | gStyle->SetPadColor(10); | |
541 | gStyle->SetPadBorderSize(1); | |
542 | gStyle->SetPadBorderMode(-1); | |
543 | ||
544 | gStyle->SetStatColor(10); | |
545 | gStyle->SetTitleColor(kBlack,"X"); | |
546 | gStyle->SetTitleColor(kBlack,"Y"); | |
547 | ||
548 | gStyle->SetLabelSize(0.04,"X"); | |
549 | gStyle->SetLabelSize(0.04,"Y"); | |
550 | gStyle->SetLabelSize(0.04,"Z"); | |
551 | gStyle->SetTitleSize(0.04,"X"); | |
552 | gStyle->SetTitleSize(0.04,"Y"); | |
553 | gStyle->SetTitleSize(0.04,"Z"); | |
554 | gStyle->SetTitleFont(42,"X"); | |
555 | gStyle->SetTitleFont(42,"Y"); | |
556 | gStyle->SetTitleFont(42,"X"); | |
557 | gStyle->SetLabelFont(42,"X"); | |
558 | gStyle->SetLabelFont(42,"Y"); | |
559 | gStyle->SetLabelFont(42,"Z"); | |
560 | gStyle->SetStatFont(42); | |
561 | ||
562 | gStyle->SetTitleOffset(1.0,"X"); | |
563 | gStyle->SetTitleOffset(1.4,"Y"); | |
564 | ||
565 | gStyle->SetFillColor(kWhite); | |
566 | gStyle->SetTitleFillColor(kWhite); | |
567 | ||
568 | gStyle->SetOptDate(0); | |
569 | gStyle->SetOptTitle(1); | |
570 | gStyle->SetOptStat(0); | |
571 | gStyle->SetOptFit(111); | |
572 | ||
573 | } | |
574 | ||
efb6c264 | 575 | // For different tasks (e.g. triggers) |
576 | void MakeMmixPi0() | |
577 | { | |
5fb23089 | 578 | // Take care when uncommenting, seems to leak memory consumption. |
579 | ||
2eb453ea | 580 | // // All Centrality |
5fb23089 | 581 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow/PHOSPi0FlowCoutput1", -1, "CPV", "out/kCentral/"); |
2eb453ea | 582 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_SemiCentral/PHOSPi0Flow_SemiCentralCoutput1", -1, "CPV", "out/kSemiCentral/"); |
583 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_MB/PHOSPi0Flow_MBCoutput1", -1, "CPV", "out/kMB/"); | |
584 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_PHOSPb/PHOSPi0Flow_PHOSPbCoutput1", -1, "CPV", "out/kPHOSPb/"); | |
5fb9fcd6 | 585 | |
2eb453ea | 586 | // 0-10% |
ebfd4ce0 | 587 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow/PHOSPi0FlowCoutput1", 0, "CPV", "out/kCentral/"); |
2eb453ea | 588 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_MB/PHOSPi0Flow_MBCoutput1", 0, "CPV", "out/kMB/"); |
589 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_PHOSPb/PHOSPi0Flow_PHOSPbCoutput1", 0, "CPV", "out/kPHOSPb/"); | |
590 | ||
591 | // // 10-40% | |
5fb23089 | 592 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_SemiCentral/PHOSPi0Flow_SemiCentralCoutput1", 1, "CPV", "out/kSemiCentral/"); |
2eb453ea | 593 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_MB/PHOSPi0Flow_MBCoutput1", 1, "CPV", "out/kMB/"); |
594 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_PHOSPb/PHOSPi0Flow_PHOSPbCoutput1", 1, "CPV", "out/kPHOSPb/"); | |
595 | ||
596 | // // 40-80% | |
597 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_SemiCentral/PHOSPi0Flow_SemiCentralCoutput1", 2, "CPV", "out/kSemiCentral/"); | |
598 | // MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_MB/PHOSPi0Flow_MBCoutput1", 2, "CPV", "out/kMB/"); | |
ebfd4ce0 | 599 | MakeMmixPi0("AnalysisResults.root", "PHOSPi0Flow_PHOSPb/PHOSPi0Flow_PHOSPbCoutput1", 2, "CPV", "out/kPHOSPb/"); |
2eb453ea | 600 | } |
601 | ||
602 | ||
603 | ||
604 | TH2F* FindPi0(TList *histoList, bool mix, const char* pid, int centrality) | |
605 | { | |
606 | // consider mixed or not | |
607 | TString mixStr(""); | |
608 | if( mix ) | |
609 | mixStr = "Mi"; | |
5fb9fcd6 | 610 | |
2eb453ea | 611 | // If centrality is some integer in range [0, kNCentralityBins] return Pi0 plot for that cent. |
612 | if( centrality >= 0 && centrality < kNCentralityBins ) { | |
613 | TString name = Form("h%sPi0%s_cen%d", mixStr.Data(), pid, centrality); | |
614 | TH2F* hist = (TH2F*) histoList->FindObject(name.Data()); | |
615 | hist->Sumw2(); | |
616 | return hist; | |
5fb9fcd6 | 617 | } |
2eb453ea | 618 | // If centrality is '-1' Merge [0, kNCentralityBins) |
619 | else if ( centrality == -1 ) { | |
620 | TString name = Form("h%sPi0%s_cen%d", mixStr.Data(), pid, 0); | |
621 | TH2F* histMerge = (TH2F*) histoList->FindObject(name.Data()) -> Clone(Form("h%sPi0%s_merged", mixStr.Data(), pid)); | |
622 | histMerge->Sumw2(); | |
5fb9fcd6 | 623 | |
2eb453ea | 624 | for( int cent = 1; cent < kNCentralityBins; ++cent ) { |
625 | name = Form("h%sPi0%s_cen%d", mixStr.Data(), pid, cent); | |
626 | TH2F* other = (TH2F*) histoList->FindObject(name.Data()); | |
627 | other->Sumw2(); | |
628 | histMerge->Add( other ); | |
629 | } | |
630 | return histMerge; | |
efb6c264 | 631 | } |
ebfd4ce0 | 632 | else { // case not defined |
2eb453ea | 633 | Printf("ERROR: Centrality must be in range: [%d,%d]", -1, kNCentralityBins - 1 ); |
ebfd4ce0 | 634 | return 0; |
635 | } | |
efb6c264 | 636 | } |