Corrected spelling mistake Emperical -> Empirical
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / corrections / Emperical / GenerateDispVtxCorr.C
1 /** 
2  * Open a file 
3  * 
4  * @param filename Name of file 
5  * 
6  * @return Return value 
7  */
8 TFile* OpenFile(const char* filename)
9 {
10   TFile* file = TFile::Open(filename, "READ");
11   if (!file) 
12     Error("OpenFile", "Failed to open the file %s for reading", filename);
13   return file;
14 }
15 /** 
16  * Get an object from parent list 
17  * 
18  * @param l    Parent list 
19  * @param name Name of object to find 
20  * 
21  * @return 0 in case of problems, object pointer otherwise 
22  */
23 TObject* GetObject(const TList* l, const char* name)
24 {
25   if (!l) {
26     Warning("GetObject", "No parent list given");
27     return 0;
28   }
29   TObject* o = l->FindObject(name);
30   if (!o) { 
31     Warning("GetObject", "Object named %s not found in list %s", 
32             name, l->GetName());
33     return 0;
34   }
35   return o;
36 }
37
38 /** 
39  * Get a list from another list 
40  * 
41  * @param l    Parent list 
42  * @param name Name of list to find 
43  * 
44  * @return 0 in case of problems, object pointer otherwise 
45  */
46 TList* GetList(const TList* l, const char* name)
47 {
48   TObject* o = GetObject(l, name);
49   if (!o) return 0;
50   return static_cast<TList*>(o);
51 }
52 /** 
53  * Get a histogram from a list 
54  * 
55  * @param l    Parent list 
56  * @param name Name of histogram to find 
57  * 
58  * @return 0 in case of problems, object pointer otherwise 
59  */
60 TH1* GetHist(const TList* l, const char* name)
61 {
62   TObject* o = GetObject(l, name);
63   if (!o) return 0;
64   return static_cast<TH1*>(o);
65 }         
66 /** 
67  * Generate the emperical correction for unknown material in the
68  * AliROOT geometric description.  This is done by taking the ratio
69  * of the nominal vertex results to the results from satelitte events. 
70  *
71  * This can be done in two ways 
72  *
73  * - Either by using FMD results only 
74  * - Or by using the full combined FMD+V0+SPD results 
75  * 
76  * The correction is written to the file "EmpiricalCorrection.root",
77  * which will contain a TGraphErrors for each defined centrality
78  * class, plus a TGraphErrors object that contains the average over
79  * all centrality classes.
80  *
81  * Other results should be DIVIDED by the correction obtained from
82  * this script.
83  * 
84  * @param fmdfmd If true, use FMD results only 
85  */
86 void 
87 GenerateDispVtxCorr(Bool_t fmdfmd=true,
88                     const char* fmdDispVtx="forward_dndetaDispVtx.root",
89                     const char* fmdNomVtx="forward_dndetaNominalVtxZDC.root",
90                     const char* fullDispVtx="dndetaPreliminaryQM12.root")
91 {
92   
93   TString outName    = "EmpiricalCorrection.root";
94   TFile* fdispvtxFMD = OpenFile(fmdDispVtx);
95   TFile* fNominalFMD = OpenFile(fmdNomVtx);
96   TFile* fdispvtxAll = OpenFile(fullDispVtx);
97   if (!fdispvtxFMD || !fNominalFMD || !fdispvtxAll) return;
98  
99   TGraphErrors* corrcent[4];
100   TMultiGraph*  mg = new TMultiGraph();
101   mg->SetName("corrections");
102   
103   Int_t limits[] = {     0,        5,      10,         20,    30  };
104   Int_t colors[] = {kRed+2, kGreen+2, kBlue+1, kMagenta+1, kBlack };
105
106   TGraphErrors* allsym     = 0;
107   TGraphErrors* fmddispvtx = 0;
108   TGraphErrors* fmdnominal = 0;
109
110   TList* displist = static_cast<TList*>(fdispvtxFMD->Get("ForwardResults"));
111   TList* nomlist  = static_cast<TList*>(fNominalFMD->Get("ForwardResults"));
112     
113   Int_t nMin = 1000;
114   // Get each defined centrality and form the ratio 
115   for(Int_t i=0; i<4; i++) {
116     Int_t cl = limits[i];
117     Int_t ch = limits[i+1];
118     corrcent[i] = new TGraphErrors;
119     corrcent[i]->SetName(Form("correction_cent_%03d_%03d",cl, ch));
120     corrcent[i]->SetTitle(Form("%2d%% - %2d%%", cl, ch));
121     corrcent[i]->GetHistogram()->SetXTitle("#eta");
122     corrcent[i]->GetHistogram()->SetYTitle("#frac{dN_{ch}/d#eta|_{nominal}}{"
123                                            "dN_{ch}/d#eta|_{satelitte}}");  
124     corrcent[i]->SetMarkerColor(colors[i]);
125     corrcent[i]->SetLineColor(colors[i]);
126     corrcent[i]->SetMarkerStyle(8);
127     corrcent[i]->SetFillStyle(0);
128     corrcent[i]->SetFillColor(0);
129
130     mg->Add(corrcent[i]);
131     TGraphErrors* allsym = 
132       static_cast<TGraphErrors*>(fdispvtxAll->Get(Form("graphErrors_cent_%d_%d",
133                                                        cl, ch)));
134     
135     TString folderName   = Form("cent%03d_%03d",cl, ch);
136     TList*  dispcentlist = GetList(displist, folderName);
137     TList*  nomcentlist  = GetList(nomlist, folderName);
138     TString histName     = "dndetaForward_rebin05";
139     TH1D*   hDisp        = GetHist(dispcentlist, histName);
140     TH1D*   hNominal     = GetHist(nomcentlist, histName);
141     
142     // Make our temporary graph 
143     if   (fmdfmd) fmddispvtx = new TGraphErrors(hDisp);
144     else          fmddispvtx = static_cast<TGraphErrors*>(allsym->Clone());
145     fmdnominal               = new TGraphErrors(hNominal);
146     
147     Int_t nPoints = 0;
148     
149     for(Int_t n=0; n<fmdnominal->GetN(); n++) {
150       Double_t eta        = fmdnominal->GetX()[n];
151       Double_t nommult    = fmdnominal->GetY()[n];
152       Double_t nommulterr = fmdnominal->GetErrorY(n);
153       Double_t etaerr     = fmdnominal->GetErrorX(n);
154
155       // Ignore empty bins 
156       if(nommult < 0.0001) continue;
157
158       // Find the corresponding bin from the dispaclaced vertex analysis
159       for(Int_t m=0; m < fmddispvtx->GetN(); m++) {
160         Double_t eta1        = fmddispvtx->GetX()[m];
161         Double_t dispmult    = fmddispvtx->GetY()[m];
162         Double_t dispmulterr = fmddispvtx->GetErrorY(m);
163
164         if(TMath::Abs(eta-eta1) >= 0.001) continue;
165
166         Double_t corr  = nommult/dispmult;
167         Double_t rd    = dispmulterr/dispmult;
168         Double_t rn    = nommulterr/nommult;
169         Double_t error = (1/corr)*TMath::Sqrt(TMath::Power(rd,2) + 
170                                               TMath::Power(rn,2));
171         corrcent[i]->SetPoint(nPoints,eta,corr);
172         corrcent[i]->SetPointError(nPoints,etaerr,error);
173         nPoints++;
174       
175       }
176       //std::cout<<eta<<"  "<<nommult<<std::endl;
177     }
178     nMin = TMath::Min(nPoints, nMin);
179   }
180
181   // Calulate the average 
182   TGraphErrors* average = new TGraphErrors;
183   average->SetName(Form("correction_average"));
184   average->SetTitle(Form("%2d%% - %2d%%", limits[0], limits[4]));
185   average->SetMarkerColor(colors[4]);
186   average->SetLineColor(colors[4]);
187   average->SetMarkerStyle(8);
188   average->SetFillStyle(0);
189   average->SetFillColor(0);
190   mg->Add(average);
191
192   for(Int_t k = 0; k < nMin; k++) {
193     Double_t mean   = 0;
194     Double_t error2 = 0;
195     Double_t sumw2  = 0;
196     
197     // Loop over centralities 
198     for(Int_t l=0; l<4; l++) {
199       Double_t eta  = corrcent[l]->GetX()[k];
200       Double_t corr = corrcent[l]->GetY()[k];
201       Double_t err  = corrcent[l]->GetErrorY(k);
202       Double_t err2 = err * err;
203       Double_t w    = 1 / err2;
204       sumw2         += w;
205       mean          += w * corr;
206 #if 0
207       mean   += corr;
208       error2 += TMath::Power(err,2);
209 #endif
210     }
211     mean           /= sumw2;
212     Double_t error = TMath::Sqrt(1. / sumw2);
213 #if 0
214     mean           /= 4;
215     Double_t error =  TMath::Sqrt(error2) / 4;
216 #endif
217     
218     average->SetPoint(k,eta,mean);
219     average->SetPointError(k,0.125,error);    
220   }
221
222   Double_t min = +1000;
223   Double_t max = -1000;
224   TMultiGraph* ratios = new TMultiGraph;
225   for(Int_t l=0; l<4; l++) {
226     Int_t cl        = limits[l];
227     Int_t ch        = limits[l+1];
228     TGraphErrors* r = 
229       static_cast<TGraphErrors*>(corrcent[l]->Clone(Form("ratio%02d%02d",
230                                                         cl,ch)));
231     ratios->Add(r);
232     for (Int_t k = 0; k < r->GetN(); k++) { 
233       Double_t x = r->GetX()[k];
234       Double_t y = r->GetY()[k];
235       Double_t a = average->Eval(x);
236       Double_t s = (y-a)/a;
237       r->SetPoint(k, x, s);
238       min  = TMath::Min(s, min);
239       max  = TMath::Max(s, max);
240     }
241   }
242   Printf("Min=%f max=%f", min, max);
243   // Draw the results 
244   TCanvas* canvas = new TCanvas("overview","overview",800,600);
245   TPad* overview =  new TPad("overview", "Overview", 0, .3, 1, 1, 0, 0);
246   overview->SetBottomMargin(0);
247   overview->SetTopMargin(0.02);
248   overview->SetRightMargin(0.02);
249   overview->Draw();
250   overview->SetGridx();
251   overview->cd();
252
253   mg->Draw("APEL");
254   mg->GetXaxis()->SetTitle("#eta");
255   mg->GetYaxis()->SetTitle("#frac{dN_{ch}/d#eta|_{nominal}}{"
256                            "dN_{ch}/d#eta|_{satelitte}}");  
257   overview->Clear();
258   mg->DrawClone("APEC");
259   TLegend* leg = overview->BuildLegend(0.35,0.6,0.6,0.975, "Centrality");
260   leg->SetFillColor(0);
261   leg->SetFillStyle(0);
262   leg->SetBorderSize(0);
263
264   canvas->cd();
265   TPad* details =  new TPad("details", "Details", 0, 0, 1, .3, 0, 0);
266   details->SetBottomMargin(0.15);
267   details->SetTopMargin(0);
268   details->SetRightMargin(0.02);
269   details->Draw();
270   details->SetGridx();
271   details->cd();
272
273
274   ratios->Draw("APEL");
275   ratios->GetXaxis()->SetTitle("#eta");
276   ratios->GetYaxis()->SetTitle("#frac{c-#LTc#GT}{#LTc#GT}");  
277   ratios->GetYaxis()->SetTitleSize(0.09);
278   ratios->GetYaxis()->SetTitleOffset(0.5);
279   ratios->GetYaxis()->SetLabelSize(0.08);
280   ratios->GetXaxis()->SetTitleSize(0.09);
281   ratios->GetXaxis()->SetTitleOffset(0.5);
282   ratios->GetXaxis()->SetLabelSize(0.08);
283   ratios->GetYaxis()->SetNdivisions(10);
284
285   Double_t x1 = ratios->GetXaxis()->GetXmin();
286   Double_t x2 = ratios->GetXaxis()->GetXmax();
287   Double_t y1 = min; // ratios->GetHistogram()->GetMinimum();
288   Double_t y2 = max; // ratios->GetHistogram()->GetMaximum();
289   TGraphAsymmErrors* band = new TGraphAsymmErrors(2);
290   band->SetName("band");
291   band->SetTitle(Form("Error band min=%4.2f%% max=%4.2f%%", -100*min, 100*max));
292   band->SetPoint(0, x1, 0);
293   band->SetPoint(1, x2, 0);
294   band->SetPointError(0, 0, 0, -y1, y2);
295   band->SetPointError(1, 0, 0, -y1, y2);
296   band->SetFillColor(kYellow+2);
297   band->SetFillStyle(3001);
298   band->SetLineStyle(2);
299   band->SetLineColor(kBlack);
300   band->Draw("A E3 L");
301   band->GetXaxis()->SetTitle("#eta");
302   band->GetYaxis()->SetTitle("#frac{c-#LTc#GT}{#LTc#GT}");  
303   band->GetYaxis()->SetTitleSize(0.09);
304   band->GetYaxis()->SetTitleOffset(0.5);
305   band->GetYaxis()->SetLabelSize(0.08);
306   band->GetXaxis()->SetTitleSize(0.09);
307   band->GetXaxis()->SetTitleOffset(0.5);
308   band->GetXaxis()->SetLabelSize(0.08);
309   band->GetYaxis()->SetNdivisions(10);
310
311   details->Clear();
312   // band->DrawClone("A E3 L");
313   ratios->GetListOfGraphs()->AddFirst(band, "E3 L");
314   ratios->DrawClone("PEC same");
315
316   leg = details->BuildLegend(0.35,0.6,0.6,0.975, "Centrality");
317   leg->SetFillColor(0);
318   leg->SetFillStyle(0);
319   leg->SetBorderSize(0);
320   
321   // Store the results 
322   TFile* fout = TFile::Open(outName,"RECREATE");
323
324   for(Int_t p=0; p<4; p++) corrcent[p]->Write();
325   average->Write("average");
326   mg->Write("all");
327   ratios->Write("all");
328   fout->ls();
329   fout->Close();
330   Info("GenerateDispVtxCorr", "Corrections written to %s", outName.Data());
331 }
332 //
333 // EOF
334 //