Fixes
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / DrawRes.C
1 #include <TH1D.h>
2 #include <TH2D.h>
3 #include <TCanvas.h>
4 #include <TPad.h>
5 #include <TFile.h>
6 #include <TTree.h>
7 #include <TError.h>
8 #include <TStyle.h>
9 #include <THStack.h>
10 #include <TLegend.h>
11 #include <TMath.h>
12 #include <TMultiGraph.h>
13 #include <TGraph.h>
14 #include <TGraphErrors.h>
15 #include <TLatex.h>
16 #include "AliAODForwardMult.h"
17 #include "OtherData.C"
18
19 /** 
20  * Draw the data stored in the AOD 
21  *
22  * To use this, do 
23  * 
24  * @code 
25  *   Root> .L $ALICE_ROOT/PWG2/FORWARD/analysis2/Compile.C
26  *   Root> Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/DrawRes.C")
27  *   Root> DrawRes dr
28  *   Root> dr.Run("AliAODs.root",-10,10,5,AliAODForwardMult::kInel,900)
29  * @endcode 
30  * 
31  * The output is stored in a ROOT file 
32  * 
33  * See also the script Pass2.C 
34  * 
35  * @ingroup pwg2_forward_analysis_scripts
36  */
37 struct DrawRes 
38 {
39 public: 
40   /** AOD tree */
41   TTree*             fTree;
42   /** AOD object */
43   AliAODForwardMult* fAOD;
44   /** Output file */
45   TFile*             fOut;
46   /** Summed histogram */
47   TH2D*              fSum;
48   /** Vertex efficiency */
49   Double_t           fVtxEff;
50   /** Title to put on the plot */
51   TString            fTitle;
52
53   //__________________________________________________________________
54   /** 
55    * Constructor 
56    * 
57    */
58   DrawRes()
59     : fTree(0), 
60       fAOD(0),
61       fOut(0), 
62       fSum(0),
63       fVtxEff(0),
64       fTitle("")
65   {}
66   //__________________________________________________________________
67   /** 
68    * Reset internal structures 
69    * 
70    */
71   void Clear(Option_t* )
72   {
73     if (fTree && fTree->GetCurrentFile()) { 
74       fTree->GetCurrentFile()->Close();
75       delete fTree;
76     }
77     if (fOut) {
78       fOut->Close();
79       delete fOut;
80     }
81     if (fSum) delete fSum;
82     fTree   = 0;
83     fOut    = 0;
84     fSum    = 0;
85     fVtxEff = 0;
86     
87   }
88   //__________________________________________________________________
89   /** 
90    * Run 
91    * 
92    * @param file   Input file with AOD tree
93    * @param vzMin  Minimum interaction point z coordinate 
94    * @param vzMax  Maximum interaction point z coordinate 
95    * @param rebin  How many times to re-bin the @f$ dN_{ch}/d\eta@f$
96    * @param mask   Trigger mask 
97    * @param energy Collision energy @f$ \sqrt{s}@f$ or @f$\sqrt{s_{NN}}@f$ 
98    * @param title  Title to put on the plot
99    * 
100    * @return True on success, false otherwise 
101    */
102   Bool_t Run(const char* file="AliAODs.root", 
103              Double_t vzMin=-10, Double_t vzMax=10, Int_t rebin=1,
104              Int_t mask=AliAODForwardMult::kInel, Int_t energy=900,
105              const char* title="")
106   {
107     TString trgName;
108     if    (mask & AliAODForwardMult::kInel)    trgName.Append("inel-");
109     if    (mask & AliAODForwardMult::kInelGt0) trgName.Append("inelgt0-");
110     if    (mask & AliAODForwardMult::kNSD)     trgName.Append("nsd-");
111     if (trgName.EndsWith("-")) trgName.Remove(trgName.Length()-1);
112     if (trgName.IsNull()) trgName = "unknown";
113     TString outName = 
114       TString::Format("dndeta_%04dGeV_%c%02d-%c%02dcm_rb%02d_%s.root",
115                       energy, 
116                       vzMin < 0 ? 'm' : 'p', Int_t(TMath::Abs(vzMin)+.5), 
117                       vzMax < 0 ? 'm' : 'p', Int_t(TMath::Abs(vzMax)+.5), 
118                       rebin, trgName.Data());
119     fTitle = title;
120     if (!Open(file, outName)) return kFALSE;
121     if (!Process(vzMin,vzMax,mask)) return kFALSE;
122     if (!Finish(rebin, mask, energy)) return kFALSE;
123
124     return kTRUE;
125   }
126   //__________________________________________________________________
127   /** 
128    * Open the files @a fname containg the tree with AliAODEvent objects, 
129    * and also open the output file @a outname for writting. 
130    * 
131    * @param fname   Input file name 
132    * @param outname Output file name 
133    * 
134    * @return true on success, false otherwise 
135    */
136   Bool_t Open(const char* fname, const char* outname) 
137   {
138     Clear("");
139     
140     TFile* file = TFile::Open(fname, "READ");
141     if (!file) { 
142       Error("Open", "Cannot open file %s", fname);
143       return kFALSE;
144     }
145
146     // Get the AOD tree 
147     fTree = static_cast<TTree*>(file->Get("aodTree"));
148     if (!fTree) {
149       Error("Init", "Couldn't get the tree");
150       return kFALSE;
151     }
152
153     // Set the branch pointer 
154     fTree->SetBranchAddress("Forward", &fAOD);
155
156     fOut = TFile::Open(outname, "RECREATE");
157     if (!fOut) { 
158       Error("Open", "Couldn't open %s", outname);
159       return kFALSE;
160     }
161     return kTRUE;
162   }
163   //__________________________________________________________________
164   /** 
165    * Process the events 
166    * 
167    * @param vzMin  Minimum interaction point z coordinate 
168    * @param vzMax  Maximum interaction point z coordinate 
169    * @param mask   Trigger mask 
170    *
171    * @return true on success, false otherwise 
172    */
173   Bool_t Process(Double_t vzMin, Double_t vzMax, Int_t mask) 
174   {
175     Int_t nTriggered  = 0;                    // # of triggered ev.
176     Int_t nWithVertex = 0;                    // # of ev. w/vertex
177     Int_t nAccepted   = 0;                    // # of ev. used
178     Int_t nAvailable  = fTree->GetEntries();  // How many entries
179  
180     for (int i = 0; i < nAvailable; i++) { 
181       fTree->GetEntry(i);
182
183       // Create sum histogram on first event - to match binning to input
184       if (!fSum) {
185         fSum = static_cast<TH2D*>(fAOD->GetHistogram().Clone("d2ndetadphi"));
186         Info("Process", "Created sum histogram (%d,%f,%f)x(%d,%f,%f)", 
187              fSum->GetNbinsX(), 
188              fSum->GetXaxis()->GetXmin(),
189              fSum->GetXaxis()->GetXmax(),
190              fSum->GetNbinsY(), 
191              fSum->GetYaxis()->GetXmin(),
192              fSum->GetYaxis()->GetXmax());
193       }
194  
195       // fAOD->Print();
196       // Other trigger/event requirements could be defined 
197       if (!fAOD->IsTriggerBits(mask)) continue; 
198       nTriggered++;
199
200       // Check if there's a vertex 
201       if (!fAOD->HasIpZ()) continue; 
202       nWithVertex++;
203
204       // Select vertex range (in centimeters) 
205       if (!fAOD->InRange(vzMin, vzMax)) continue; 
206       nAccepted++;
207  
208       // Add contribution from this event
209       fSum->Add(&(fAOD->GetHistogram()));
210     }
211     fVtxEff = Double_t(nWithVertex)/nTriggered;
212
213     Info("Process", "Total of %9d events\n"
214          "               %9d has a trigger\n" 
215          "               %9d has a vertex\n" 
216          "               %9d was used\n",
217          nAvailable, nTriggered, nWithVertex, nAccepted);
218
219     return kTRUE;
220   }
221   //__________________________________________________________________
222   /** 
223    * Finish the stuff and draw 
224    * 
225    * @param rebin  How many times to rebin
226    * @param energy Collision energy 
227    * @param mask   Trigger mask 
228    * 
229    * @return true on success, false otherwise 
230    */
231   Bool_t Finish(Int_t rebin, Int_t mask, Int_t energy)
232   {
233     fOut->cd();
234
235     // Get acceptance normalisation from underflow bins 
236     TH1D* norm   = fSum->ProjectionX("norm", 0, 1, "");
237     // Project onto eta axis - _ignoring_underflow_bins_!
238     TH1D* dndeta = fSum->ProjectionX("dndeta", 1, -1, "e");
239     dndeta->SetTitle("1/N dN_{ch}/d#eta|_{forward}");
240     // Normalize to the acceptance 
241     dndeta->Divide(norm);
242     // Scale by the vertex efficiency 
243     dndeta->Scale(fVtxEff, "width");
244     dndeta->SetMarkerColor(kRed+1);
245     dndeta->SetMarkerStyle(20);
246     dndeta->SetMarkerSize(1);
247     dndeta->SetFillStyle(0);
248
249     TH1* sym = Symmetrice(dndeta);
250     
251     Rebin(dndeta, rebin);
252     Rebin(sym, rebin);
253
254     THStack* stack = new THStack("results", "Results");
255     stack->Add(dndeta);
256     stack->Add(sym);
257
258     TString hhdName(fOut->GetName());
259     hhdName.ReplaceAll("dndeta", "hhd");    
260     TH1* hhd    = GetHHD(hhdName.Data(), mask & AliAODForwardMult::kNSD);
261     TH1* hhdsym = 0;
262     if (hhd) { 
263       hhdsym = Symmetrice(hhd);
264       stack->Add(hhd);
265       stack->Add(hhdsym);
266     }
267
268     TMultiGraph* mg     = GetData(energy, mask);
269
270     gStyle->SetOptTitle(0);
271     Double_t yd = (mg->GetListOfGraphs()->GetEntries() ? 0.3 : 0);
272     TCanvas* c = new TCanvas("c", "C", 900, 800);
273     c->SetFillColor(0);
274     c->SetBorderMode(0);
275     c->SetBorderSize(0);
276     c->cd();
277
278     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0);
279     p1->SetTopMargin(0.05);
280     p1->SetBottomMargin(0.001);
281     p1->SetRightMargin(0.05);
282     p1->SetGridx();
283     p1->SetTicks(1,1);
284     p1->Draw();
285     p1->cd();
286     stack->SetMinimum(-0.1);
287     FixAxis(stack, 1/(1-yd)/1.5, "#frac{1}{N} #frac{dN_{ch}}{#eta}");
288     p1->Clear();
289     stack->DrawClone("nostack e1");
290
291
292     THStack* ratios = new THStack("ratios", "Ratios");
293     TGraphAsymmErrors* o = 0;
294     TIter nextG(mg->GetListOfGraphs());
295     while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
296       o->Draw("same p");
297       ratios->Add(Ratio(dndeta, o));
298       ratios->Add(Ratio(sym,    o));
299       ratios->Add(Ratio(hhd,    o));
300       ratios->Add(Ratio(hhdsym, o));
301     }
302     if (hhd && hhdsym) { 
303       TH1F* t1 = static_cast<TH1F*>(dndeta->Clone(Form("%s_%s", 
304                                                        dndeta->GetName(), 
305                                                        hhd->GetName())));
306       t1->SetTitle(Form("%s / %s", dndeta->GetTitle(), hhd->GetTitle()));
307       TH1F* t2 = static_cast<TH1F*>(sym->Clone(Form("%s_%s", 
308                                                     sym->GetName(), 
309                                                     hhdsym->GetName())));
310       t2->SetTitle(Form("%s / %s", sym->GetTitle(), hhdsym->GetTitle()));
311       t1->Divide(hhd);
312       t2->Divide(hhdsym);
313       ratios->Add(t1);
314       ratios->Add(t2);
315     }
316
317     // Make a legend
318     TString trigs(AliAODForwardMult::GetTriggerString(mask));
319     TLegend* l = p1->BuildLegend(.3,p1->GetBottomMargin()+.01,.7,.4,
320                                  Form("#sqrt{s}=%dGeV, %s", 
321                                       energy, trigs.Data()));
322     l->SetFillColor(0);
323     l->SetFillStyle(0);
324     l->SetBorderSize(0);
325     p1->cd();
326     c->cd();
327
328     if (!ratios->GetHists() || ratios->GetHists()->GetEntries() <= 0) {
329       p1->SetPad(0, 0, 1, 1);
330       p1->SetBottomMargin(0.1);
331       l->SetY1(0.11);
332       FixAxis(stack, (1-yd)/1,  "#frac{1}{N} #frac{dN_{ch}}{#eta}",false);
333       p1->cd();
334       c->cd();
335     }
336     else {
337       TPad* p2 = new TPad("p2", "p2", 0, 0.0, 1.0, yd, 0, 0);
338       p2->SetTopMargin(0.001);
339       p2->SetRightMargin(0.05);
340       p2->SetBottomMargin(1/yd * 0.07);
341       p2->SetGridx();
342       p2->SetTicks(1,1);
343       p2->Draw();
344       p2->cd();
345       FixAxis(ratios, 1/yd/1.5, "Ratios");
346       p2->Clear();
347       ratios->DrawClone("nostack e1");
348       
349       TLegend* l2 = p2->BuildLegend(.3,p2->GetBottomMargin()+.01,.7,.6);
350       l2->SetFillColor(0);
351       l2->SetFillStyle(0);
352       l2->SetBorderSize(0);
353       
354       p2->cd();
355       TGraphErrors* band = new TGraphErrors(2);
356       band->SetPoint(0, sym->GetXaxis()->GetXmin(), 1);
357       band->SetPoint(1, dndeta->GetXaxis()->GetXmax(), 1);
358       band->SetPointError(0, 0, .1);
359       band->SetPointError(1, 0, .1);
360       band->SetFillColor(kYellow+2);
361       band->SetFillStyle(3002);
362       band->Draw("3 same");
363       ratios->DrawClone("nostack e1 same");
364
365       c->cd();
366     }
367     p1->cd();
368     TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
369     tit->SetNDC();
370     tit->SetTextSize(0.05);
371     tit->Draw();
372     
373     TString imgName(fOut->GetName());
374     imgName.ReplaceAll(".root", ".png");
375     c->SaveAs(imgName.Data());
376
377     stack->Write();
378     mg->Write();
379     ratios->Write();
380
381     fOut->Close();
382
383     return kTRUE;
384   }
385   //__________________________________________________________________
386   /** 
387    * Get the result from previous analysis code 
388    * 
389    * @param fn File to open 
390    * 
391    * @return null or result of previous analysis code 
392    */
393   TH1* GetHHD(const char* fn="fmd_dNdeta_mult.root", bool nsd=false) 
394   {
395     TDirectory* savdir = gDirectory;
396     TFile* file = TFile::Open(fn, "READ");
397     if (!file) { 
398       Warning("GetHHD", "couldn't open HHD file %s", fn);
399       savdir->cd();
400       return 0;
401     }
402     TString hist(Form("dNdeta_dNdeta%s", (nsd ? "NSD" : "")));
403     TH1* h = static_cast<TH1*>(file->Get(hist.Data()));
404     if (!h) { 
405       Warning("GetHHD", "Couldn't find HHD histogram %s in %s", 
406               hist.Data(), fn);
407       file->ls();
408       file->Close();
409       savdir->cd();
410       return 0;
411     }
412     TH1* r = static_cast<TH1*>(h->Clone("dndeta_hhd"));
413     r->SetTitle("1/N dN_{ch}/d#eta (HHD)");
414     r->SetFillStyle(0);
415     r->SetFillColor(0);
416     r->SetMarkerStyle(21);
417     r->SetMarkerColor(kPink+1);
418     r->SetDirectory(savdir);
419
420     file->Close();
421     savdir->cd();
422     return r;
423   }
424   //__________________________________________________________________
425   /** 
426    * Fix the apperance of the axis in a stack 
427    * 
428    * @param stack  stack of histogram
429    * @param s      Scaling factor 
430    * @param ytitle Y axis title 
431    * @param force  Whether to draw the stack first or not 
432    */
433   void FixAxis(THStack* stack, Double_t s, const char* ytitle, 
434                Bool_t force=true) 
435   {
436     if (force) stack->Draw("nostack e1");
437
438     TH1* h = stack->GetHistogram();
439     if (!h) return;
440
441     h->SetXTitle("#eta");
442     h->SetYTitle(ytitle);
443     TAxis* xa = h->GetXaxis();
444     TAxis* ya = h->GetYaxis();
445     if (xa) { 
446       xa->SetTitle("#eta");
447       // xa->SetTicks("+-");
448       xa->SetTitleSize(s*xa->GetTitleSize());
449       xa->SetLabelSize(s*xa->GetLabelSize());
450       xa->SetTickLength(s*xa->GetTickLength());
451     }
452     if (ya) { 
453       ya->SetTitle(ytitle);
454       // ya->SetTicks("+-");
455       ya->SetNdivisions(10);
456       ya->SetTitleSize(s*ya->GetTitleSize());
457       ya->SetLabelSize(s*ya->GetLabelSize());
458     }      
459   }
460   //__________________________________________________________________
461   /** 
462    * Compute the ratio of @a h to @a g.  @a g is evaluated at the bin
463    * centers of @a h
464    * 
465    * @param h  Numerator 
466    * @param g  Divisor 
467    * 
468    * @return h/g 
469    */
470   TH1* Ratio(const TH1* h, const TGraph* g) const 
471   {
472     if (!h || !g) return 0;
473
474     TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
475     ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
476     ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
477     ret->Reset();
478     ret->SetMarkerStyle(h->GetMarkerStyle());
479     ret->SetMarkerColor(g->GetMarkerColor());
480     ret->SetMarkerSize(0.7*g->GetMarkerSize());
481     Double_t xlow  = g->GetX()[0];
482     Double_t xhigh = g->GetX()[g->GetN()-1];
483     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
484
485     for (Int_t i = 1; i <= h->GetNbinsX(); i++) { 
486       Double_t c = h->GetBinContent(i);
487       if (c <= 0) continue;
488
489       Double_t x = h->GetBinCenter(i);
490       if (x < xlow || x > xhigh) continue; 
491
492       Double_t f = g->Eval(x);
493       if (f <= 0) continue; 
494
495       ret->SetBinContent(i, c / f);
496       ret->SetBinError(i, h->GetBinError(i) / f);
497     }
498     if (ret->GetEntries() <= 0) { delete ret; ret = 0; }
499     return ret;
500   }
501
502   //__________________________________________________________________
503   /** 
504    * Make an extension of @a h to make it symmetric about 0 
505    * 
506    * @param h Histogram to symmertrice 
507    * 
508    * @return Symmetric extension of @a h 
509    */
510   TH1* Symmetrice(const TH1* h) const
511   {
512     Int_t nBins = h->GetNbinsX();
513     TH1*  s     = new TH1D(Form("%s_mirror", h->GetName()),
514                            Form("%s (mirrored)", h->GetTitle()), 
515                            nBins, 
516                            -h->GetXaxis()->GetXmax(), 
517                            -h->GetXaxis()->GetXmin());
518     s->SetMarkerColor(h->GetMarkerColor());
519     s->SetMarkerSize(h->GetMarkerSize());
520     s->SetMarkerStyle(h->GetMarkerStyle()+4);
521     s->SetFillColor(h->GetFillColor());
522     s->SetFillStyle(h->GetFillStyle());
523
524     // Find the first and last bin with data 
525     Int_t first = nBins+1;
526     Int_t last  = 0;
527     for (Int_t i = 1; i <= nBins; i++) { 
528       if (h->GetBinContent(i) <= 0) continue; 
529       first = TMath::Min(first, i);
530       last  = TMath::Max(last,  i);
531     }
532     
533     Double_t xfirst = h->GetBinCenter(first-1);
534     Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
535     Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
536     for (Int_t i = f1, j=l2; i <= last; i++,j--) { 
537       s->SetBinContent(j, h->GetBinContent(i));
538       s->SetBinError(j, h->GetBinError(i));
539     }
540     return s;
541   }
542   //__________________________________________________________________
543   /** 
544    * Rebin a histogram 
545    * 
546    * @param h     Histogram to rebin
547    * @param rebin Rebinning factor 
548    * 
549    * @return 
550    */
551   virtual void Rebin(TH1* h, Int_t rebin) const
552   { 
553     if (rebin <= 1) return;
554
555     Int_t nBins = h->GetNbinsX();
556     if(nBins % rebin != 0) {
557       Warning("Rebin", "Rebin factor %d is not a devisor of current number "
558               "of bins %d in the histogram %s", rebin, nBins, h->GetName());
559       return;
560     }
561     
562     // Make a copy 
563     TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
564     tmp->Rebin(rebin);
565     tmp->SetDirectory(0);
566
567     // The new number of bins 
568     Int_t nBinsNew = nBins / rebin;
569     for(Int_t i = 1;i<= nBinsNew; i++) {
570       Double_t content = 0;
571       Double_t sumw    = 0;
572       Double_t wsum    = 0;
573       Int_t    nbins   = 0;
574       for(Int_t j = 1; j<=rebin;j++) {
575         Int_t bin = (i-1)*rebin + j;
576         if(h->GetBinContent(bin) <= 0) continue;
577         Double_t c =  h->GetBinContent(bin);
578         Double_t w = 1 / TMath::Power(c,2);
579         content    += c;
580         sumw       += w;
581         wsum       += w * c;
582         nbins++;
583       }
584       
585       if(content > 0 ) {
586         tmp->SetBinContent(i, wsum / sumw);
587         tmp->SetBinError(i,TMath::Sqrt(sumw));
588       }
589     }
590
591     // Finally, rebin the histogram, and set new content
592     h->Rebin(rebin);
593     for(Int_t i =1;i<=nBinsNew; i++) {
594       h->SetBinContent(i,tmp->GetBinContent(i));
595       // h->SetBinError(i,tmp->GetBinError(i));
596     }
597     
598     delete tmp;
599   }
600 };
601
602 //____________________________________________________________________
603 //
604 // EOF
605 //
606