Corrected spelling mistake Emperical -> Empirical
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / corrections / Empirical / GenerateEmpirical.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  * Get an object from parent list 
39  * 
40  * @param l    Parent list 
41  * @param name Name of object to find 
42  * 
43  * @return 0 in case of problems, object pointer otherwise 
44  */
45 TObject* GetObject(const TDirectory* l, const char* name)
46 {
47   if (!l) {
48     Warning("GetObject", "No parent directory given");
49     return 0;
50   }
51   TObject* o = l->Get(name);
52   if (!o) { 
53     Warning("GetObject", "Object named %s not found in directory %s", 
54             name, l->GetName());
55     return 0;
56   }
57   return o;
58 }
59 Bool_t CheckType(const TObject* o, const TClass* cl)
60 {
61   if (!o) return false;
62   if (!o->IsA()->InheritsFrom(cl)) { 
63     Warning("CheckType", "Object %s is not a %s, but a %s", 
64             o->GetName(), cl->GetName(), o->ClassName());
65     return false;
66   }
67   return true;
68 }
69 /** 
70  * Get a list from another list 
71  * 
72  * @param l    Parent list 
73  * @param name Name of list to find 
74  * 
75  * @return 0 in case of problems, object pointer otherwise 
76  */
77 TList* GetList(const TList* l, const char* name)
78 {
79   TObject* o = GetObject(l, name);
80   if (!CheckType(o, TList::Class())) return 0;
81   return static_cast<TList*>(o);
82 }
83 /** 
84  * Get a list from another list 
85  * 
86  * @param l    Parent list 
87  * @param name Name of list to find 
88  * 
89  * @return 0 in case of problems, object pointer otherwise 
90  */
91 TList* GetList(const TDirectory* l, const char* name)
92 {
93   TObject* o = GetObject(l, name);
94   if (!CheckType(o, TList::Class())) return 0;
95   return static_cast<TList*>(o);
96 }
97 /** 
98  * Get a histogram from a list 
99  * 
100  * @param l    Parent list 
101  * @param name Name of histogram to find 
102  * 
103  * @return 0 in case of problems, object pointer otherwise 
104  */
105 TH1* GetHist(const TList* l, const char* name)
106 {
107   TObject* o = GetObject(l, name);
108   if (!CheckType(o, TH1::Class())) return 0;
109   return static_cast<TH1*>(o);
110 }         
111 /** 
112  * Get a histogram from a list 
113  * 
114  * @param l    Parent list 
115  * @param name Name of histogram to find 
116  * 
117  * @return 0 in case of problems, object pointer otherwise 
118  */
119 TH1* GetHist(const TDirectory* l, const char* name)
120 {
121   TObject* o = GetObject(l, name);
122   if (!CheckType(o, TH1::Class())) return 0;
123   return static_cast<TH1*>(o);
124 }         
125 /** 
126  * Get a histogram from a list 
127  * 
128  * @param l    Parent list 
129  * @param name Name of histogram to find 
130  * 
131  * @return 0 in case of problems, object pointer otherwise 
132  */
133 TGraphErrors* GetGraph(const TList* l, const char* name)
134 {
135   TObject* o = GetObject(l, name);
136   if (!CheckType(o, TGraphErrors::Class())) return 0;
137   return static_cast<TGraphErrors*>(o);
138 }         
139 /** 
140  * Get a histogram from a list 
141  * 
142  * @param l    Parent list 
143  * @param name Name of histogram to find 
144  * 
145  * @return 0 in case of problems, object pointer otherwise 
146  */
147 TGraphErrors* GetGraph(const TDirectory* l, const char* name)
148 {
149   TObject* o = GetObject(l, name);
150   if (!CheckType(o, TGraphErrors::Class())) return 0;
151   return static_cast<TGraphErrors*>(o);
152 }         
153
154 TGraphErrors* Ratio(const TGraphErrors* num, const TGraphErrors* denom)
155 {
156   TGraphErrors* r = static_cast<TGraphErrors*>(num->Clone("tmp"));
157   for (Int_t k = 0; k < r->GetN(); k++) { 
158     Double_t xn = r->GetX()[k];
159     Double_t yn = r->GetY()[k];
160     Double_t en = r->GetErrorY(k);
161     // Double_t xd = denom->GetX()[k];
162     Double_t yd = denom->GetY()[k];
163     Double_t ed = denom->GetErrorY(k);
164     Double_t s  = (yn-yd)/yd;
165     Double_t e  = 0;
166     r->SetPoint(k, xn, s);
167     r->SetPointError(k, 0, e);
168   }
169   return r;
170 }
171
172 /** 
173  * Generate the emperical correction for unknown material in the
174  * AliROOT geometric description.  This is done by taking the ratio
175  * of the nominal vertex results to the results from satelitte events. 
176  *
177  * This can be done in two ways 
178  *
179  * - Either by using FMD results only 
180  * - Or by using the full combined FMD+V0+SPD results 
181  * 
182  * The correction is written to the file "EmpiricalCorrection.root",
183  * which will contain a TGraphErrors for each defined centrality
184  * class, plus a TGraphErrors object that contains the average over
185  * all centrality classes.
186  *
187  * Other results should be DIVIDED by the correction obtained from
188  * this script.
189  * 
190  * @param fmdfmd If true, use FMD results only 
191  */
192 void 
193 GenerateEmpirical(const char* fmdDispVtx="forward_dndetaDispVtx.root",
194                   const char* fmdNomVtx="forward_dndetaNominalVtxZDC.root",
195                   const char* fullDispVtx="dndetaPreliminaryQM12.root")
196 {
197   
198   TString outName    = "EmpiricalCorrection.root";
199   TFile* fdispvtxFMD = OpenFile(fmdDispVtx);
200   TFile* fNominalFMD = OpenFile(fmdNomVtx);
201   TFile* fdispvtxAll = OpenFile(fullDispVtx);
202   if (!fdispvtxFMD || !fNominalFMD || !fdispvtxAll) return;
203
204   TFile* out = TFile::Open(outName, "RECREATE");
205   if (!out) { 
206     Error("GenerateEmpirical", "Failed to open output file %s", 
207           outName.Data());
208     return;
209   }
210
211   Int_t limits[] = {     0,        5,      10,         20,    30  };
212   Int_t colors[] = {kRed+2, kGreen+2, kBlue+1, kMagenta+1, kBlack };
213
214   TGraphErrors* afmdfmd  = 0;
215   TGraphErrors* afmdfull = 0;
216
217   // --- Do two things fmd-to-fmd and fmd-to-full
218   for (Int_t iM = 0; iM < 2; iM++) {
219     Bool_t       fmdfmd   = (iM == 0);
220     TMultiGraph* mg       = new TMultiGraph();
221     TMultiGraph* mgfmdnom = (iM == 0 ? new TMultiGraph() : 0);
222     TMultiGraph* mgref    = new TMultiGraph();
223     mgref->SetTitle(Form("Satellite, FMD%s", fmdfmd ? "" : "+V0+Tracklets"));
224     if (mgfmdnom) mgfmdnom->SetTitle("Nominal FMD");
225
226     Info("", "Doing %s/FMD", fmdfmd ? "FMD" : "Full");
227     TGraphErrors* corrcent[4];    
228     TGraphErrors* allsym     = 0;
229     TGraphErrors* fmddispvtx = 0;
230     TGraphErrors* fmdnominal = 0;
231
232     TList* displist = GetList(fdispvtxFMD,"ForwardResults");
233     TList* nomlist  = GetList(fNominalFMD,"ForwardResults");
234     
235     Int_t nMin = 1000;
236     // --- Get each defined centrality and form the ratio ------------
237     for(Int_t iC=0; iC<4; iC++) {
238       Int_t cl = limits[iC];
239       Int_t ch = limits[iC+1];
240       TString base; base.Form("_cent_%03d_%03d", cl, ch);
241       corrcent[iC] = new TGraphErrors;
242       corrcent[iC]->SetName(Form("correction%s",base.Data()));
243       corrcent[iC]->SetTitle(Form("%2d%% - %2d%%", cl, ch));
244       corrcent[iC]->GetHistogram()->SetXTitle("#eta");
245       corrcent[iC]->GetHistogram()->SetYTitle("#frac{dN_{ch}/d#eta|_{nominal}}{"
246                                              "dN_{ch}/d#eta|_{satelitte}}");  
247       corrcent[iC]->SetMarkerColor(colors[iC]);
248       corrcent[iC]->SetLineColor(colors[iC]);
249       corrcent[iC]->SetMarkerStyle(fmdfmd ? 20 : 24);
250       corrcent[iC]->SetMarkerSize(fmdfmd ? 1 : 1.2);
251       corrcent[iC]->SetFillStyle(0);
252       corrcent[iC]->SetFillColor(0);
253       
254       mg->Add(corrcent[iC]);
255       TGraphErrors* allsym = GetGraph(fdispvtxAll,Form("graphErrors_cent_%d_%d",
256                                                        cl, ch));
257     
258       TString folderName   = Form("cent%03d_%03d",cl, ch);
259       TList*  dispcentlist = GetList(displist, folderName);
260       TList*  nomcentlist  = GetList(nomlist, folderName);
261       TString histName     = "dndetaForward_rebin05";
262       TH1D*   hDisp        = GetHist(dispcentlist, histName);
263       TH1D*   hNominal     = GetHist(nomcentlist, histName);
264     
265       // Make our temporary graph 
266       if   (fmdfmd) fmddispvtx = new TGraphErrors(hDisp);
267       else          fmddispvtx = static_cast<TGraphErrors*>(allsym->Clone());
268       fmdnominal               = new TGraphErrors(hNominal);
269       if (mgfmdnom) {
270         TGraph* g = static_cast<TGraph*>(fmdnominal->Clone(Form("nominal%s", 
271                                                                 base.Data())));
272         g->SetMarkerColor(colors[iC]);
273         g->SetLineColor(colors[iC]);
274         g->SetMarkerStyle(21);
275         g->SetTitle("Nominal FMD");
276         mgfmdnom->Add(g);
277       }
278       TGraph* ref = static_cast<TGraph*>(fmddispvtx->Clone(Form("satellite%s",
279                                                                 base.Data())));
280       ref->SetTitle(Form("%2d - %2d", cl, ch));
281       ref->SetMarkerColor(colors[iC]);
282       ref->SetLineColor(colors[iC]);
283       ref->SetMarkerStyle(fmdfmd ? 22 : 20);
284       ref->SetMarkerSize(1);
285       mgref->Add(ref);
286     
287       Int_t nPoints = 0;
288     
289       for(Int_t iN=0; iN<fmdnominal->GetN(); iN++) {
290         Double_t eta        = fmdnominal->GetX()[iN];
291         Double_t nommult    = fmdnominal->GetY()[iN];
292         Double_t nommulterr = fmdnominal->GetErrorY(iN);
293         Double_t etaerr     = fmdnominal->GetErrorX(iN);
294
295         // Ignore empty bins 
296         if(nommult < 0.0001) continue;
297
298         // Find the corresponding bin from the dispaclaced vertex analysis
299         for(Int_t iS=0; iS < fmddispvtx->GetN(); iS++) {
300           Double_t eta1        = fmddispvtx->GetX()[iS];
301           Double_t dispmult    = fmddispvtx->GetY()[iS];
302           Double_t dispmulterr = fmddispvtx->GetErrorY(iS);
303
304           if(TMath::Abs(eta-eta1) >= 0.001) continue;
305           
306           Double_t corr  = nommult/dispmult;
307           Double_t rd    = dispmulterr/dispmult;
308           Double_t rn    = nommulterr/nommult;
309           Double_t error = (1/corr)*TMath::Sqrt(TMath::Power(rd,2) + 
310                                               TMath::Power(rn,2));
311           corrcent[iC]->SetPoint(nPoints,eta,corr);
312           corrcent[iC]->SetPointError(nPoints,etaerr,error);
313           nPoints++;
314           
315         }
316         //std::cout<<eta<<"  "<<nommult<<std::endl;
317       }
318       nMin = TMath::Min(nPoints, nMin);
319     }
320
321     // --- Calulate the average --------------------------------------
322     TGraphErrors* average = new TGraphErrors;
323     average->SetName(Form("average"));
324     average->SetTitle(Form("%2d%% - %2d%%", limits[0], limits[4]));
325     average->SetMarkerColor(colors[4]);
326     average->SetLineColor(colors[4]);
327     average->SetMarkerStyle(fmdfmd ? 20 : 24);
328     average->SetMarkerSize(fmdfmd ? 1 : 1.2);
329     average->SetFillStyle(0);
330     average->SetFillColor(0);
331     mg->Add(average);
332     
333     for(Int_t iA = 0; iA < nMin; iA++) {
334       Double_t mean   = 0;
335       Double_t error2 = 0;
336       Double_t sumw2  = 0;
337       
338       // Loop over centralities 
339       for(Int_t iC=0; iC<4; iC++) {
340         Double_t eta  = corrcent[iC]->GetX()[iA];
341         Double_t corr = corrcent[iC]->GetY()[iA];
342         Double_t err  = corrcent[iC]->GetErrorY(iA);
343         Double_t err2 = err * err;
344         Double_t w    = 1 / err2;
345         sumw2         += w;
346         mean          += w * corr;
347       }
348       mean           /= sumw2;
349       Double_t error = TMath::Sqrt(1. / sumw2);
350           average->SetPoint(iA,eta,mean);
351       average->SetPointError(iA,0.125,error);    
352     }
353     if (fmdfmd) afmdfmd  = average;
354     else        afmdfull = average;
355
356     // --- Calculate ratios ------------------------------------------
357     TMultiGraph* ratios = new TMultiGraph;
358     ratios->SetName("ratios");
359     for(Int_t iC=0; iC<4; iC++) {
360       Int_t cl        = limits[iC];
361       Int_t ch        = limits[iC+1];
362       TGraphErrors* r = Ratio(corrcent[iC], average);
363       r->SetName(Form("ratio%s",base.Data()));
364       ratios->Add(r);
365     }
366
367     // --- Store the result ------------------------------------------
368     TDirectory* d = out->mkdir(fmdfmd ? "fmdfmd" : "fmdfull", 
369                                Form("Empirical corrections %s", 
370                                     fmdfmd ? "FMD/FMD" : "FULL/FMD"));
371     d->cd();
372     for(Int_t p=0; p<4; p++) corrcent[p]->Write();
373     average->Write("average");
374     mg->Write("all");
375     ratios->Write("ratios");
376     mgref->Write("satellite");
377     out->cd();
378     if (mgfmdnom) mgfmdnom->Write("nominal");
379     out->ls();
380   } // End method loop
381
382   TGraphErrors* r = Ratio(afmdfull, afmdfmd);
383   r->Write("ratio");
384
385   out->Close();
386 }
387 //
388 // EOF
389 //
390