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