commented define _ClusterTopology_ - to be used only for the special productions
[u/mrichter/AliRoot.git] / ITS / ParametrizeITSBetheBloch.C
CommitLineData
fba679b6 1#if !defined(__CINT__) || defined(__MAKECINT__)
2#include <TCanvas.h>
3#include <TGraph.h>
4#include <TGraphErrors.h>
5#include <TGraphAsymmErrors.h>
6#include <TH1.h>
7#include <TF1.h>
8#include <TH2.h>
9#include <TFile.h>
10#include <TList.h>
11#include <TString.h>
12#include <TDirectoryFile.h>
13#include <TLatex.h>
14#include <TStyle.h>
15#include <TMath.h>
16#include <TDatabasePDG.h>
17#include "AliITSPIDResponse.h"
18#endif
19
20/* $Id$ */
21
22
23// Macro to extract the parameterizations of dE/dx in ITS vs. p
24// using Phobos BB parameterization
25
26
27Double_t BetheBlochPhobos(const Double_t *x, const Double_t *par);
28Double_t BetheBlochPhobosBG(const Double_t *x, const Double_t *par);
29Double_t BetheBlochAleph(const Double_t *x, const Double_t *par);
30
31
32void ParametrizeITSBetheBloch(Bool_t optFromITSsa=kFALSE){
33
34 TFile *fil=new TFile("AnalysisResults.root");
35 TH2F *hsigvsp=0x0;
36 TString dfName="PWG2SpectraITSsa";
37 TString listName="clistITSsaMult-1to-1";
38 TString hisName="fHistDEDX";
39
40 if(optFromITSsa){
41 dfName="ITSsaTracks";
42 listName="clistITSsaTracks";
43 hisName="hdedxvsPt4clsITSpureSA";
44 }
45 TDirectoryFile* df=(TDirectoryFile*)fil->Get(dfName.Data());
46 if(!df){
47 printf("ERROR: directory file %s not found\n",dfName.Data());
48 return;
49 }
50 TList* l=(TList*)df->Get(listName.Data());
51 if(!l){
52 printf("ERROR: TList %s not found\n",listName.Data());
53 return;
54 }
55 hsigvsp=(TH2F*)l->FindObject(hisName.Data());
56 if(!hsigvsp){
57 printf("ERROR: historgam %s not found\n",hisName.Data());
58 return;
59 }
60 hsigvsp->GetXaxis()->SetTitle("momentum (GeV/c)");
61 hsigvsp->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
62
63 Double_t minpPionFit=0.1;
64 Double_t maxpPionFit=0.8;
65 Double_t minbgGlobFit=0.3;
66 Double_t maxbgGlobFit=2.5;
67 Double_t minbgPions=2.;
68 Double_t maxbgPions=10.;
69 Double_t minpKaons=0.2;
70 Double_t maxpKaons=0.4;
71 Double_t maxbgKaons=0.7;
72 Double_t minpProtons=0.35;
73 Double_t maxpProtons=0.57;
74 Double_t maxbgProtons=0.7;
75
76 Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
77 Double_t massKaon=TDatabasePDG::Instance()->GetParticle(321)->Mass();
78 Double_t massProton=TDatabasePDG::Instance()->GetParticle(2212)->Mass();
79 Double_t massDeuteron=1.8756;
80
81 Int_t iPad=1;
82 Int_t iCanvas=0;
83 Int_t nbinsX=hsigvsp->GetNbinsX();
84 printf("Total Number of bins =%d\n",nbinsX);
85 Int_t group=2;
86
87 Int_t totPads=(Int_t)((Float_t)nbinsX/(Float_t)group+0.5);
88 Int_t totCanvas=totPads/16;
89 TCanvas** c0=new TCanvas*[totCanvas];
90 c0[iCanvas]=new TCanvas(Form("c0_%d",iCanvas),Form("c0_%d",iCanvas));
91 c0[iCanvas]->Divide(4,4);
92 TCanvas** c0k=new TCanvas*[totCanvas];
93 c0k[iCanvas]=new TCanvas(Form("c0k_%d",iCanvas),Form("c0k_%d",iCanvas));
94 c0k[iCanvas]->Divide(4,4);
95 TCanvas** c0p=new TCanvas*[totCanvas];
96 c0p[iCanvas]=new TCanvas(Form("c0p_%d",iCanvas),Form("c0p_%d",iCanvas));
97 c0p[iCanvas]->Divide(4,4);
98
99 TGraphErrors* gdedxpi=new TGraphErrors(0);
100 TGraphErrors* gdedxk=new TGraphErrors(0);
101 TGraphErrors* gdedxp=new TGraphErrors(0);
102 TGraphErrors* gdedxbg=new TGraphErrors(0);
103 gdedxbg->SetTitle("");
104
105 Int_t nPpi=0;
106 Int_t nPk=0;
107 Int_t nPp=0;
108 Int_t nPbg=0;
109 for(Int_t iBin=1; iBin<=nbinsX; iBin+=group){
110 if(iPad>16){
111 iCanvas++;
112 c0[iCanvas]=new TCanvas(Form("c0_%d",iCanvas),Form("c0_%d",iCanvas));
113 c0[iCanvas]->Divide(4,4);
114 c0k[iCanvas]=new TCanvas(Form("c0k_%d",iCanvas),Form("c0k_%d",iCanvas));
115 c0k[iCanvas]->Divide(4,4);
116 c0p[iCanvas]=new TCanvas(Form("c0p_%d",iCanvas),Form("c0p_%d",iCanvas));
117 c0p[iCanvas]->Divide(4,4);
118 iPad=1;
119 }
120
121 Float_t sum=0;
122 Float_t sumwei=0;
123 Double_t pmed=0.;
124 for(Int_t ibx=iBin; ibx<iBin+group;ibx++){
125 for(Int_t iby=1; iby<=hsigvsp->GetNbinsY();iby++){
126 sum+=hsigvsp->GetXaxis()->GetBinCenter(ibx)*hsigvsp->GetBinContent(ibx,iby);
127 sumwei+=hsigvsp->GetBinContent(ibx,iby);
128 }
129 }
130 if(sumwei>0.) pmed=sum/sumwei;
131 if(pmed<0.1) continue;
132 Int_t iPmed=(Int_t)(pmed*1000);
133 TH1D* hProj=hsigvsp->TH2F::ProjectionY(Form("proj%d",iPmed),iBin,iBin+group-1);
134 TH1D* hKaon=(TH1D*)hProj->Clone(Form("projk%d",iPmed));
135 c0[iCanvas]->cd(iPad);
136 hProj->SetLineColor(2);
137 hProj->Draw();
138
139 Float_t peakpos=hProj->GetBinCenter(hProj->GetMaximumBin());
140 Float_t peaksig=0.1*peakpos;
141 hProj->Fit("gaus","LL","",peakpos-3*peaksig,peakpos+3.*peaksig);
142 TF1* fgaus=(TF1*)hProj->GetListOfFunctions()->FindObject("gaus");
143 gdedxpi->SetPoint(nPpi,sum/sumwei,fgaus->GetParameter(1));
144 gdedxpi->SetPointError(nPpi,0.,fgaus->GetParError(1));
145 nPpi++;
146 if(pmed/massPion>minbgPions && pmed/massPion<maxbgPions){
147 gdedxbg->SetPoint(nPbg,pmed/massPion,fgaus->GetParameter(1));
148 gdedxbg->SetPointError(nPbg,0.,fgaus->GetParError(1));
149 nPbg++;
150 }
151 for(Int_t iBinK=1; iBinK<=hKaon->GetNbinsX(); iBinK++){
152 Double_t cont=hKaon->GetBinContent(iBinK)-fgaus->Eval(hKaon->GetBinCenter(iBinK));
153 if(cont<0) cont=0.;
154 if(hKaon->GetBinCenter(iBinK)<fgaus->GetParameter(1)) cont=0;
155 hKaon->SetBinContent(iBinK,cont);
156 }
157 c0k[iCanvas]->cd(iPad);
158 hKaon->SetLineColor(4);
159 hKaon->Draw();
160 TH1D* hProton=(TH1D*)hKaon->Clone(Form("projp%d",iPmed));
161 TF1* fgausk=0x0;
162 Double_t mink=250-600*(pmed-0.2);
163 Double_t maxk=mink+mink;
164 hKaon->Fit("gaus","LL","",mink,maxk);
165 fgausk=(TF1*)hKaon->GetListOfFunctions()->FindObject("gaus");
166 if(pmed>minpKaons && pmed<maxpKaons){
167 gdedxk->SetPoint(nPk,pmed,fgausk->GetParameter(1));
168 gdedxk->SetPointError(nPk,0.,fgausk->GetParError(1));
169 nPk++;
170 if(pmed/massKaon<maxbgKaons){
171 gdedxbg->SetPoint(nPbg,pmed/massKaon,fgausk->GetParameter(1));
172 gdedxbg->SetPointError(nPbg,0.,fgausk->GetParError(1));
173 nPbg++;
174 }
175 }
176 if(fgausk){
177 for(Int_t iBinP=1; iBinP<=hProton->GetNbinsX(); iBinP++){
178 Double_t cont=hKaon->GetBinContent(iBinP)-fgausk->Eval(hKaon->GetBinCenter(iBinP));
179 if(cont<0) cont=0.;
180 if(hProton->GetBinCenter(iBinP)<fgausk->GetParameter(1)) cont=0;
181 hProton->SetBinContent(iBinP,cont);
182 }
183 }
184 c0p[iCanvas]->cd(iPad);
185 hProton->SetLineColor(kGreen+1);
186 hProton->Draw();
187 if(pmed>minpProtons && pmed<maxpProtons){
188 Double_t minP=300-600*(pmed-0.35);
189 Double_t maxP=minP+minP;
190 hProton->Fit("gaus","LL","",minP,maxP);
191 TF1* fgausp=(TF1*)hProton->GetListOfFunctions()->FindObject("gaus");
192 gdedxp->SetPoint(nPp,pmed,fgausp->GetParameter(1));
193 gdedxp->SetPointError(nPp,0.,fgausp->GetParError(1));
194 nPp++;
195 if(pmed/massProton<maxbgProtons){
196 gdedxbg->SetPoint(nPbg,pmed/massProton,fgausp->GetParameter(1));
197 gdedxbg->SetPointError(nPbg,0.,fgausp->GetParError(1));
198 nPbg++;
199 }
200 }
201 ++iPad;
202 }
203
204 gStyle->SetPalette(1);
205 hsigvsp->GetXaxis()->SetRangeUser(0.07,1.8);
206 hsigvsp->SetStats(0);
207 TCanvas* c1=new TCanvas("c1","Pion Fit");
208 c1->SetLogz();
209 hsigvsp->Draw("col");
210 gdedxpi->Draw("Psame");
211 gdedxk->Draw("Psame");
212 gdedxp->Draw("Psame");
213
214 TF1 *fPhobos=new TF1("fPhobos",BetheBlochPhobos,minpPionFit,maxpPionFit,5);
215 fPhobos->SetParLimits(0,0.,1.E9);
216 fPhobos->SetParLimits(1,0.,100.);
217 fPhobos->SetParLimits(2,0.,1.);
218 fPhobos->SetParLimits(3,0.,1.);
219 fPhobos->SetParLimits(4,0.,1.);
220 gdedxpi->Fit("fPhobos","R");
221
222 TCanvas* c1g=new TCanvas("c1g","Global Fit");
223 TF1 *fPhobosbg=new TF1("fPhobosbg",BetheBlochPhobosBG,minbgGlobFit,maxbgGlobFit,5);
224 fPhobosbg->SetParLimits(0,0.,1.E9);
225 fPhobosbg->SetParLimits(1,0.,100.);
226 fPhobosbg->SetParLimits(2,0.,1.);
227 fPhobosbg->SetParLimits(3,0.,1.);
228 fPhobosbg->SetParLimits(4,0.,1.);
229 gdedxbg->Draw("AP");
230 gdedxbg->GetXaxis()->SetTitle("#beta#gamma");
231 gdedxbg->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
232 gdedxbg->Fit("fPhobosbg","R");
233 c1g->Update();
234
235
236 AliITSPIDResponse* itsresGlobFit=new AliITSPIDResponse(kFALSE);
237 AliITSPIDResponse* itsres=new AliITSPIDResponse(kFALSE);
238 Double_t parameters[5],parametersGlob[5];
239
240 printf("-------- Pion Fit Parameters -----------\n");
241 for(Int_t iPar=0; iPar<5; iPar++){
242 parameters[iPar]=fPhobos->GetParameter(iPar);
243 printf("parameters[%d]=%g;\n",iPar,parameters[iPar]);
244 }
245 printf("-------- Global Fit Parameters -----------\n");
246 for(Int_t iPar=0; iPar<5; iPar++){
247 parametersGlob[iPar]=fPhobosbg->GetParameter(iPar);
248 printf("parametersGlob[%d]=%g;\n",iPar,parametersGlob[iPar]);
249 }
250
251 itsresGlobFit->SetBetheBlochParamsITSsa(parametersGlob);
252 itsres->SetBetheBlochParamsITSsa(parameters);
253
254 TGraph* gPion=new TGraph(0);
255 TGraph* gKaon=new TGraph(0);
256 TGraph* gProton=new TGraph(0);
257 TGraph* gPionGlobFit=new TGraph(0);
258 TGraph* gKaonGlobFit=new TGraph(0);
259 TGraph* gProtonGlobFit=new TGraph(0);
260
261 TGraphAsymmErrors* gBBpions=new TGraphAsymmErrors(0);
262 TGraphAsymmErrors* gBBkaons=new TGraphAsymmErrors(0);
263 TGraphAsymmErrors* gBBprotons=new TGraphAsymmErrors(0);
264 TGraphAsymmErrors* gBBpionsGlobFit=new TGraphAsymmErrors(0);
265 TGraphAsymmErrors* gBBkaonsGlobFit=new TGraphAsymmErrors(0);
266 TGraphAsymmErrors* gBBprotonsGlobFit=new TGraphAsymmErrors(0);
267
268 Double_t numberofsigmapions=2.;
269 Float_t resodedx[4];
270 resodedx[0]=0.0001; //default value
271 resodedx[1]=0.0001; //default value
272 resodedx[2]=0.116;
273 resodedx[3]=0.103;
274
275
276 for(Int_t ip=0; ip<400; ip++){
277 Double_t mom=0.1+0.02*ip;
278 Double_t dedxpions=itsres->Bethe(mom,massPion,kTRUE);
279 Double_t dedxkaons=itsres->Bethe(mom,massKaon,kTRUE);
280 Double_t dedxprotons=itsres->Bethe(mom,massProton,kTRUE);
281 Double_t dedxdeutons=itsres->Bethe(mom,massDeuteron,kTRUE);
282 Double_t dedxpionsGlobFit=itsresGlobFit->Bethe(mom,massPion,kTRUE);
283 Double_t dedxkaonsGlobFit=itsresGlobFit->Bethe(mom,massKaon,kTRUE);
284 Double_t dedxprotonsGlobFit=itsresGlobFit->Bethe(mom,massProton,kTRUE);
285 Double_t dedxdeutonsGlobFit=itsresGlobFit->Bethe(mom,massDeuteron,kTRUE);
286
287 gPion->SetPoint(ip,mom,dedxpions);
288 gKaon->SetPoint(ip,mom,dedxkaons);
289 gProton->SetPoint(ip,mom,dedxprotons);
290 gPionGlobFit->SetPoint(ip,mom,dedxpionsGlobFit);
291 gKaonGlobFit->SetPoint(ip,mom,dedxkaonsGlobFit);
292 gProtonGlobFit->SetPoint(ip,mom,dedxprotonsGlobFit);
293
294 gBBpions->SetPoint(ip,mom,dedxpions);
295 gBBpions->SetPointError(ip,0,0,(dedxpions*resodedx[3]*numberofsigmapions),TMath::Abs((dedxpions-dedxkaons)/2));
296 gBBkaons->SetPoint(ip,mom,dedxkaons);
297 gBBkaons->SetPointError(ip,0,0,TMath::Abs((dedxpions-dedxkaons)/2),TMath::Abs((dedxprotons-dedxkaons)/2));
298 gBBprotons->SetPoint(ip,mom,dedxprotons);
299 gBBprotons->SetPointError(ip,0,0,TMath::Abs((dedxprotons-dedxkaons)/2),TMath::Abs((dedxprotons-dedxdeutons)/2));
300
301 gBBpionsGlobFit->SetPoint(ip,mom,dedxpionsGlobFit);
302 gBBpionsGlobFit->SetPointError(ip,0,0,(dedxpionsGlobFit*resodedx[3]*numberofsigmapions),TMath::Abs((dedxpions-dedxkaonsGlobFit)/2));
303 gBBkaonsGlobFit->SetPoint(ip,mom,dedxkaonsGlobFit);
304 gBBkaonsGlobFit->SetPointError(ip,0,0,TMath::Abs((dedxpionsGlobFit-dedxkaonsGlobFit)/2),TMath::Abs((dedxprotonsGlobFit-dedxkaonsGlobFit)/2));
305 gBBprotonsGlobFit->SetPoint(ip,mom,dedxprotonsGlobFit);
306 gBBprotonsGlobFit->SetPointError(ip,0,0,TMath::Abs((dedxprotonsGlobFit-dedxkaonsGlobFit)/2),TMath::Abs((dedxprotonsGlobFit-dedxdeutonsGlobFit)/2));
307 }
308
309
310 TCanvas* c2=new TCanvas("c2","Parameterizations");
311 c2->SetLogz();
312 hsigvsp->SetMinimum(20);
313 hsigvsp->Draw("col");
314 gPion->Draw("LSAME");
315 gKaon->Draw("LSAME");
316 gProton->Draw("LSAME");
317 gPionGlobFit->SetLineColor(kRed+1);
318 gKaonGlobFit->SetLineColor(kRed+1);
319 gProtonGlobFit->SetLineColor(kRed+1);
320 gPionGlobFit->Draw("LSAME");
321 gKaonGlobFit->Draw("LSAME");
322 gProtonGlobFit->Draw("LSAME");
323 TLatex* t1=new TLatex(0.7,0.85,"Pion Fit");
324 t1->SetNDC();
325 t1->Draw();
326 TLatex* tG=new TLatex(0.7,0.75,"Global Fit");
327 tG->SetNDC();
328 tG->SetTextColor(kRed+1);
329 tG->Draw();
330
331
332 gBBpions->SetLineColor(2);
333 gBBpions->SetLineWidth(2);
334 gBBpions->SetFillColor(2);
335 gBBpions->SetFillStyle(3001);
336 gBBkaons->SetLineColor(3);
337 gBBkaons->SetLineWidth(2);
338 gBBkaons->SetFillColor(3);
339 gBBkaons->SetFillStyle(3001);
340 gBBprotons->SetLineColor(4);
341 gBBprotons->SetLineWidth(2);
342 gBBprotons->SetFillColor(4);
343 gBBprotons->SetFillStyle(3001);
344
345 TCanvas* c3=new TCanvas("c3","Asymm Bands, PionFit");
346 c3->SetLogz();
347 c3->SetLogx();
348 hsigvsp->SetMinimum(20);
349 hsigvsp->Draw("col");
350 gBBpions->Draw("L4same");
351 gBBkaons->Draw("L4same");
352 gBBprotons->Draw("L4same");
353
354 gBBpionsGlobFit->SetLineColor(2);
355 gBBpionsGlobFit->SetLineWidth(2);
356 gBBpionsGlobFit->SetFillColor(2);
357 gBBpionsGlobFit->SetFillStyle(3001);
358 gBBkaonsGlobFit->SetLineColor(3);
359 gBBkaonsGlobFit->SetLineWidth(2);
360 gBBkaonsGlobFit->SetFillColor(3);
361 gBBkaonsGlobFit->SetFillStyle(3001);
362 gBBprotonsGlobFit->SetLineColor(4);
363 gBBprotonsGlobFit->SetLineWidth(2);
364 gBBprotonsGlobFit->SetFillColor(4);
365 gBBprotonsGlobFit->SetFillStyle(3001);
366
367 TCanvas* c3g=new TCanvas("c3g","Asymm Bands, GlobFit");
368 c3g->SetLogz();
369 c3g->SetLogx();
370 hsigvsp->SetMinimum(20);
371 hsigvsp->Draw("col");
372 gBBpionsGlobFit->Draw("L4same");
373 gBBkaonsGlobFit->Draw("L4same");
374 gBBprotonsGlobFit->Draw("L4same");
375
376}
377
378
379
380
381
382
383Double_t BetheBlochPhobos(const Double_t *x, const Double_t *par){
384
385 //
386 // returns AliExternalTrackParam::BetheBloch normalized to
387 // fgMIP at the minimum
388 //
389
390 Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
391 Double_t bg=x[0]/massPion;
392 Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
393 Double_t gamma=bg/beta;
394 Double_t eff=1.0;
395 if(bg<par[2])
396 eff=(bg-par[3])*(bg-par[3])+par[4];
397 else
398 eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
399
400 Double_t bb=0.;
401 if(gamma>=0. && beta>0.){
402 bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
403 }
404 return bb;
405}
406
407Double_t BetheBlochPhobosBG(const Double_t *x, const Double_t *par){
408
409 //
410 // returns AliExternalTrackParam::BetheBloch normalized to
411 // fgMIP at the minimum
412 //
413
414 Double_t bg=x[0];
415 Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
416 Double_t gamma=bg/beta;
417 Double_t eff=1.0;
418 if(bg<par[2])
419 eff=(bg-par[3])*(bg-par[3])+par[4];
420 else
421 eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
422
423 Double_t bb=0.;
424 if(gamma>=0. && beta>0.){
425 bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
426 }
427 return bb;
428}
429
430Double_t BetheBlochAleph(const Double_t *x, const Double_t *par){
431
432 //
433 // This is the empirical ALEPH parameterization of the Bethe-Bloch formula.
434 // It is normalized to 1 at the minimum.
435 //
436 // bg - beta*gamma
437 //
438 // The default values for the kp* parameters are for ALICE TPC.
439 // The returned value is in MIP units
440 //
441
442 Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
443 Double_t bg=x[0]/massPion;
444 Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
445
446 Double_t aa = TMath::Power(beta,par[3]);
447 Double_t bb = TMath::Power(1./bg,par[4]);
448
449 bb=TMath::Log(par[2]+bb);
450
451 return (par[1]-aa-bb)*par[0]/aa;
452}