B2 analysis code
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Nuclei / B2 / macros / DrawSec.C
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // Draw corrections of secondaries
17 // author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
18
19 #include <Riostream.h>
20 #include <TROOT.h>
21 #include <TFile.h>
22 #include <TSystem.h>
23 #include <TH1D.h>
24 #include <TString.h>
25 #include <TMath.h>
26 #include <TCanvas.h>
27
28 #include "B2.h"
29
30
31 Int_t GetNumberCols(Int_t lowbin, Int_t hibin);
32
33 TCanvas* DrawMC(TH1D** hData, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax);
34
35 TCanvas* DrawFit(TH1D** hData, TH1D** hDataFit, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax);
36
37 TCanvas* DrawGoF(TH1D** hData, TH1D** hDataFit, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax);
38
39 void DrawSec(const TString& inputFile="debug.root", const TString& tag="test", const TString& particle="Proton", Int_t lowbin=1, Int_t hibin=10, Double_t dcaxyMin=-1.5, Double_t dcaxyMax=1.5)
40 {
41 //
42 // draw corrections of secondaries
43 //
44         const Int_t kNBin = hibin;
45         
46         TFile* finput = new TFile(inputFile.Data());
47         if (finput->IsZombie()) exit(1);
48         
49         // data and templates
50         
51         TH1D* hData[kNBin];
52         TH1D* hDataMC[kNBin];
53         TH1D* hPrim[kNBin];
54         TH1D* hMat[kNBin];
55         TH1D* hFdwn[kNBin];
56         
57         // fitted data
58         
59         TH1D* hDataFit[kNBin];
60         TH1D* hPrimFit[kNBin];
61         TH1D* hMatFit[kNBin];
62         TH1D* hFdwnFit[kNBin];
63         
64         for(Int_t i=lowbin; i<hibin && i<kNBin; ++i)
65         {
66                 // MC
67                 hData[i]   = (TH1D*)FindObj(finput, tag, particle + Form("_Data_DCAxy_%02d",i));
68                 hDataMC[i] = (TH1D*)FindObj(finput, tag, particle + Form("_SimData_DCAxy_%02d",i));
69                 hPrim[i]   = (TH1D*)FindObj(finput, tag, particle + Form("_Prim_DCAxy_%02d",i));
70                 hMat[i]    = (TH1D*)FindObj(finput, tag, particle + Form("_Mat_DCAxy_%02d",i));
71                 
72                 if(particle.Contains("Proton"))
73                 {
74                         hFdwn[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Fdwn_DCAxy_%02d",i));
75                 }
76                 else
77                 {
78                         hFdwn[i] = 0;
79                 }
80                 
81                 // fitted data
82                 hDataFit[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Fit_Data_DCAxy_%02d",i));
83                 hPrimFit[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Fit_Prim_DCAxy_%02d",i));
84                 hMatFit[i]  = (TH1D*)FindObj(finput, tag, particle + Form("_Fit_Mat_DCAxy_%02d",i));
85                 
86                 if(particle.Contains("Proton"))
87                 {
88                         hFdwnFit[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Fit_Fdwn_DCAxy_%02d",i));
89                 }
90                 else
91                 {
92                         hFdwnFit[i] = 0;
93                 }
94         }
95         
96         // pt bin width
97         TH1D* hPidPt = (TH1D*)FindObj(finput, tag, particle + "_PID_Pt");
98         Double_t binwidth = hPidPt->GetBinWidth(0);
99         
100         // draw
101         
102         TCanvas* c0 = 0;
103         
104         if(hDataMC[lowbin]!=0)
105         {
106                 c0 = DrawMC(hDataMC, hPrim, hMat, hFdwn, lowbin, hibin, particle + ".MC", Form("MC simulation for %s",particle.Data()), binwidth, dcaxyMin, dcaxyMax);
107         }
108         
109         TCanvas* c1;
110         
111         if(hDataFit[lowbin]!=0)
112         {
113                 c1 = DrawFit(hData, hDataFit, hPrimFit, hMatFit, hFdwnFit, lowbin, hibin, particle + ".Fit", Form("Fitted DCA for %s",particle.Data()), binwidth, dcaxyMin, dcaxyMax);
114         }
115         
116         TCanvas* c2;
117         
118         if(hDataFit[lowbin]!=0)
119         {
120                 c2 = DrawGoF(hData, hDataFit, lowbin, hibin, particle + ".GoF", Form("GoF for %s",particle.Data()), binwidth, dcaxyMin, dcaxyMax);
121         }
122 }
123
124 Int_t GetNumberCols(Int_t lowbin, Int_t hibin)
125 {
126 //
127 // number of rows and cols for plotting the slices
128 // square matrix (cols = rows)
129 //
130         return TMath::CeilNint(TMath::Sqrt(TMath::Abs(hibin - lowbin)));
131 }
132
133 TCanvas* DrawMC(TH1D** hData, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax)
134 {
135 //
136 // Draw MC DCA distributions for each pt bin
137 //
138         Int_t ncol = GetNumberCols(lowbin, hibin);
139         
140         TCanvas* c0 = new TCanvas(name.Data(), title.Data());
141         c0->Divide(ncol, ncol);
142         
143         for(Int_t i=lowbin; i < hibin; ++i)
144         {
145                 c0->cd(i+1-lowbin);
146                 
147                 gPad->SetGrid();
148                 gPad->SetLogy();
149                 
150                 hData[i]->SetTitle(Form("%0.2f < p_{T} < %0.2f GeV/c", (i-1)*binwidth, i*binwidth));
151                 hData[i]->SetMinimum(0.2);
152                 
153                 hData[i]->SetAxisRange(dcaxyMin, dcaxyMax, "X");
154                 
155                 hData[i]->SetMarkerStyle(kFullCircle);
156                 hData[i]->Draw("E");
157                 
158                 if(hPrim[i] != 0)
159                 {
160                         hPrim[i]->SetLineColor(9);
161                         hPrim[i]->Draw("sameHist");
162                 }
163                 
164                 if(hMat[i] != 0)
165                 {
166                         hMat[i]->SetLineColor(46);
167                         hMat[i]->Draw("sameHist");
168                 }
169                 
170                 if(hFdwn != 0)
171                 if(hFdwn[i] != 0)
172                 {
173                         hFdwn[i]->SetLineColor(8);
174                         hFdwn[i]->Draw("sameHist");
175                 }
176         }
177         
178         return c0;
179 }
180
181 TCanvas* DrawFit(TH1D** hData, TH1D** hDataFit, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax)
182 {
183 //
184 // Draw DCA fitted model for each pt bin in current canvas
185 //
186         Int_t ncol = GetNumberCols(lowbin, hibin);
187         
188         TCanvas* c0 = new TCanvas(name.Data(), title.Data());
189         c0->Divide(ncol, ncol);
190         
191         for(Int_t i=lowbin; i < hibin; ++i)
192         {
193                 c0->cd(i+1-lowbin);
194                 
195                 gPad->SetGridx();
196                 gPad->SetGridy();
197                 gPad->SetLogy();
198                 
199                 hData[i]->SetTitle(Form("%0.2f < p_{T} < %0.2f GeV/c", (i-1)*binwidth, i*binwidth));
200                 hData[i]->SetMinimum(0.2);
201                 
202                 hData[i]->SetAxisRange(dcaxyMin, dcaxyMax, "X");
203                 
204                 hData[i]->SetMarkerStyle(kFullCircle);
205                 hData[i]->Draw("E");
206                 
207                 hDataFit[i]->SetLineColor(2);
208                 hDataFit[i]->Draw("sameHist");
209                 
210                 if(hPrim[i] != 0)
211                 {
212                         hPrim[i]->SetLineColor(9);
213                         hPrim[i]->Draw("sameHist");
214                 }
215                 
216                 if(hMat[i] != 0)
217                 {
218                         hMat[i]->SetLineColor(46);
219                         hMat[i]->Draw("sameHist");
220                 }
221                 
222                 if(hFdwn != 0)
223                 if(hFdwn[i] != 0)
224                 {
225                         hFdwn[i]->SetLineColor(8);
226                         hFdwn[i]->Draw("sameHist");
227                 }
228         }
229         
230         return c0;
231 }
232
233 TCanvas* DrawGoF(TH1D** hData, TH1D** hDataFit, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax)
234 {
235 //
236 // Draw a goodness of fit for each pt bin consisting of data/fit
237 //
238         Int_t ncol = GetNumberCols(lowbin, hibin);
239         
240         TCanvas* c0 = new TCanvas(name.Data(), title.Data());
241         c0->Divide(ncol, ncol);
242         
243         for(Int_t i=lowbin; i < hibin; ++i)
244         {
245                 c0->cd(i+1-lowbin);
246                 
247                 gPad->SetGrid();
248                 
249                 TH1D* hDiv = (TH1D*)hData[i]->Clone("_Div_Data_Fit_");
250                 hDiv->Sumw2();
251                 hDiv->Divide(hDataFit[i]);
252                 
253                 hDiv->SetTitle(Form("%0.2f < p_{T} < %0.2f GeV/c", (i-1)*binwidth, i*binwidth));
254                 hDiv->SetAxisRange(dcaxyMin, dcaxyMax, "X");
255                 hDiv->SetAxisRange(0,3,"Y");
256                 hDiv->SetMarkerStyle(kFullCircle);
257                 hDiv->DrawCopy("E");
258                 
259                 delete hDiv;
260         }
261         
262         return c0;
263 }