New implementation of the forward multiplicity analysis.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / DrawRes.C
CommitLineData
7e4038b5 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 */
37struct DrawRes
38{
39public:
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 nAccepted = 0; // # of ev. w/vertex
177 Int_t nAvailable = fTree->GetEntries(); // How many entries
178
179 for (int i = 0; i < nAvailable; i++) {
180 fTree->GetEntry(i);
181
182 // Create sum histogram on first event - to match binning to input
183 if (!fSum) {
184 fSum = static_cast<TH2D*>(fAOD->GetHistogram().Clone("d2ndetadphi"));
185 Info("Process", "Created sum histogram (%d,%f,%f)x(%d,%f,%f)",
186 fSum->GetNbinsX(),
187 fSum->GetXaxis()->GetXmin(),
188 fSum->GetXaxis()->GetXmax(),
189 fSum->GetNbinsY(),
190 fSum->GetYaxis()->GetXmin(),
191 fSum->GetYaxis()->GetXmax());
192 }
193
194 // fAOD->Print();
195 // Other trigger/event requirements could be defined
196 if (!fAOD->IsTriggerBits(mask)) continue;
197 nTriggered++;
198
199 // Select vertex range (in centimeters)
200 if (!fAOD->InRange(vzMin, vzMax)) continue;
201 nAccepted++;
202
203 // Add contribution from this event
204 fSum->Add(&(fAOD->GetHistogram()));
205 }
206 fVtxEff = Double_t(nAccepted)/nTriggered;
207
208 Info("Process", "Got %6d triggers, accepted %6d out of %6d events",
209 nTriggered, nAccepted, nAvailable);
210
211 return kTRUE;
212 }
213 //__________________________________________________________________
214 /**
215 * Finish the stuff and draw
216 *
217 * @param rebin How many times to rebin
218 * @param energy Collision energy
219 * @param mask Trigger mask
220 *
221 * @return true on success, false otherwise
222 */
223 Bool_t Finish(Int_t rebin, Int_t mask, Int_t energy)
224 {
225 fOut->cd();
226
227 // Get acceptance normalisation from underflow bins
228 TH1D* norm = fSum->ProjectionX("norm", 0, 1, "");
229 // Project onto eta axis - _ignoring_underflow_bins_!
230 TH1D* dndeta = fSum->ProjectionX("dndeta", 1, -1, "e");
231 dndeta->SetTitle("1/N dN_{ch}/d#eta|_{forward}");
232 // Normalize to the acceptance
233 dndeta->Divide(norm);
234 // Scale by the vertex efficiency
235 dndeta->Scale(fVtxEff, "width");
236 dndeta->SetMarkerColor(kRed+1);
237 dndeta->SetMarkerStyle(20);
238 dndeta->SetMarkerSize(1);
239 dndeta->SetFillStyle(0);
240
241 TH1* sym = Symmetrice(dndeta);
242
243 Rebin(dndeta, rebin);
244 Rebin(sym, rebin);
245
246 THStack* stack = new THStack("results", "Results");
247 stack->Add(dndeta);
248 stack->Add(sym);
249
250 TString hhdName(fOut->GetName());
251 hhdName.ReplaceAll("dndeta", "hhd");
252 TH1* hhd = GetHHD(hhdName.Data(), mask & AliAODForwardMult::kNSD);
253 TH1* hhdsym = 0;
254 if (hhd) {
255 hhdsym = Symmetrice(hhd);
256 stack->Add(hhd);
257 stack->Add(hhdsym);
258 }
259
260 TMultiGraph* mg = GetData(energy, mask);
261
262 gStyle->SetOptTitle(0);
263 Double_t yd = (mg->GetListOfGraphs()->GetEntries() ? 0.3 : 0);
264 TCanvas* c = new TCanvas("c", "C", 900, 800);
265 c->SetFillColor(0);
266 c->SetBorderMode(0);
267 c->SetBorderSize(0);
268 c->cd();
269
270 TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0);
271 p1->SetTopMargin(0.05);
272 p1->SetBottomMargin(0.001);
273 p1->SetRightMargin(0.05);
274 p1->SetGridx();
275 p1->SetTicks(1,1);
276 p1->Draw();
277 p1->cd();
278 stack->SetMinimum(-0.1);
279 FixAxis(stack, 1/(1-yd)/1.5, "#frac{1}{N} #frac{dN_{ch}}{#eta}");
280 p1->Clear();
281 stack->DrawClone("nostack e1");
282
283
284 THStack* ratios = new THStack("ratios", "Ratios");
285 TGraphAsymmErrors* o = 0;
286 TIter nextG(mg->GetListOfGraphs());
287 while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
288 o->Draw("same p");
289 ratios->Add(Ratio(dndeta, o));
290 ratios->Add(Ratio(sym, o));
291 ratios->Add(Ratio(hhd, o));
292 ratios->Add(Ratio(hhdsym, o));
293 }
294 // Make a legend
295 TString trigs(AliAODForwardMult::GetTriggerString(mask));
296 TLegend* l = p1->BuildLegend(.3,p1->GetBottomMargin()+.01,.7,.4,
297 Form("#sqrt{s}=%dGeV, %s",
298 energy, trigs.Data()));
299 l->SetFillColor(0);
300 l->SetFillStyle(0);
301 l->SetBorderSize(0);
302 p1->cd();
303 c->cd();
304
305 if (!ratios->GetHists() || ratios->GetHists()->GetEntries() <= 0) {
306 p1->SetPad(0, 0, 1, 1);
307 p1->SetBottomMargin(0.1);
308 l->SetY1(0.11);
309 FixAxis(stack, (1-yd)/1, "#frac{1}{N} #frac{dN_{ch}}{#eta}",false);
310 p1->cd();
311 c->cd();
312 }
313 else {
314 TPad* p2 = new TPad("p2", "p2", 0, 0.0, 1.0, yd, 0, 0);
315 p2->SetTopMargin(0.001);
316 p2->SetRightMargin(0.05);
317 p2->SetBottomMargin(1/yd * 0.07);
318 p2->SetGridx();
319 p2->SetTicks(1,1);
320 p2->Draw();
321 p2->cd();
322 FixAxis(ratios, 1/yd/1.5, "Ratios");
323 p2->Clear();
324 ratios->DrawClone("nostack e1");
325
326 TLegend* l2 = p2->BuildLegend(.3,p2->GetBottomMargin()+.01,.7,.6);
327 l2->SetFillColor(0);
328 l2->SetFillStyle(0);
329 l2->SetBorderSize(0);
330
331 p2->cd();
332 TGraphErrors* band = new TGraphErrors(2);
333 band->SetPoint(0, sym->GetXaxis()->GetXmin(), 1);
334 band->SetPoint(1, dndeta->GetXaxis()->GetXmax(), 1);
335 band->SetPointError(0, 0, .1);
336 band->SetPointError(1, 0, .1);
337 band->SetFillColor(kYellow+2);
338 band->SetFillStyle(3002);
339 band->Draw("3 same");
340 ratios->DrawClone("nostack e1 same");
341
342 c->cd();
343 }
344 p1->cd();
345 TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
346 tit->SetNDC();
347 tit->SetTextSize(0.05);
348 tit->Draw();
349
350 TString imgName(fOut->GetName());
351 imgName.ReplaceAll(".root", ".png");
352 c->SaveAs(imgName.Data());
353
354 stack->Write();
355 mg->Write();
356 ratios->Write();
357
358 fOut->Close();
359
360 return kTRUE;
361 }
362 //__________________________________________________________________
363 /**
364 * Get the result from previous analysis code
365 *
366 * @param fn File to open
367 *
368 * @return null or result of previous analysis code
369 */
370 TH1* GetHHD(const char* fn="fmd_dNdeta_mult.root", bool nsd=false)
371 {
372 TDirectory* savdir = gDirectory;
373 TFile* file = TFile::Open(fn, "READ");
374 if (!file) {
375 Warning("GetHHD", "couldn't open HHD file %s", fn);
376 savdir->cd();
377 return 0;
378 }
379 TString hist(Form("dNdeta_dNdeta%s", (nsd ? "NSD" : "")));
380 TH1* h = static_cast<TH1*>(file->Get(hist.Data()));
381 if (!h) {
382 Warning("GetHHD", "Couldn't find HHD histogram %s in %s",
383 hist.Data(), fn);
384 file->ls();
385 file->Close();
386 savdir->cd();
387 return 0;
388 }
389 TH1* r = static_cast<TH1*>(h->Clone("dndeta_hhd"));
390 r->SetTitle("1/N dN_{ch}/d#eta (HHD)");
391 r->SetFillStyle(0);
392 r->SetFillColor(0);
393 r->SetMarkerStyle(21);
394 r->SetMarkerColor(kPink+1);
395 r->SetDirectory(savdir);
396
397 file->Close();
398 savdir->cd();
399 return r;
400 }
401 //__________________________________________________________________
402 /**
403 * Fix the apperance of the axis in a stack
404 *
405 * @param stack stack of histogram
406 * @param s Scaling factor
407 * @param ytitle Y axis title
408 * @param force Whether to draw the stack first or not
409 */
410 void FixAxis(THStack* stack, Double_t s, const char* ytitle,
411 Bool_t force=true)
412 {
413 if (force) stack->Draw("nostack e1");
414
415 TH1* h = stack->GetHistogram();
416 if (!h) return;
417
418 h->SetXTitle("#eta");
419 h->SetYTitle(ytitle);
420 TAxis* xa = h->GetXaxis();
421 TAxis* ya = h->GetYaxis();
422 if (xa) {
423 xa->SetTitle("#eta");
424 // xa->SetTicks("+-");
425 xa->SetTitleSize(s*xa->GetTitleSize());
426 xa->SetLabelSize(s*xa->GetLabelSize());
427 xa->SetTickLength(s*xa->GetTickLength());
428 }
429 if (ya) {
430 ya->SetTitle(ytitle);
431 // ya->SetTicks("+-");
432 ya->SetNdivisions(10);
433 ya->SetTitleSize(s*ya->GetTitleSize());
434 ya->SetLabelSize(s*ya->GetLabelSize());
435 }
436 }
437 //__________________________________________________________________
438 /**
439 * Compute the ratio of @a h to @a g. @a g is evaluated at the bin
440 * centers of @a h
441 *
442 * @param h Numerator
443 * @param g Divisor
444 *
445 * @return h/g
446 */
447 TH1* Ratio(const TH1* h, const TGraph* g) const
448 {
449 if (!h || !g) return 0;
450
451 TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
452 ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
453 ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
454 ret->Reset();
455 ret->SetMarkerStyle(h->GetMarkerStyle());
456 ret->SetMarkerColor(g->GetMarkerColor());
457 ret->SetMarkerSize(0.7*g->GetMarkerSize());
458 Double_t xlow = g->GetX()[0];
459 Double_t xhigh = g->GetX()[g->GetN()-1];
460 if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
461
462 for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
463 Double_t c = h->GetBinContent(i);
464 if (c <= 0) continue;
465
466 Double_t x = h->GetBinCenter(i);
467 if (x < xlow || x > xhigh) continue;
468
469 Double_t f = g->Eval(x);
470 if (f <= 0) continue;
471
472 ret->SetBinContent(i, c / f);
473 ret->SetBinError(i, h->GetBinError(i) / f);
474 }
475 if (ret->GetEntries() <= 0) { delete ret; ret = 0; }
476 return ret;
477 }
478
479 //__________________________________________________________________
480 /**
481 * Make an extension of @a h to make it symmetric about 0
482 *
483 * @param h Histogram to symmertrice
484 *
485 * @return Symmetric extension of @a h
486 */
487 TH1* Symmetrice(const TH1* h) const
488 {
489 Int_t nBins = h->GetNbinsX();
490 TH1* s = new TH1D(Form("%s_mirror", h->GetName()),
491 Form("%s (mirrored)", h->GetTitle()),
492 nBins,
493 -h->GetXaxis()->GetXmax(),
494 -h->GetXaxis()->GetXmin());
495 s->SetMarkerColor(h->GetMarkerColor());
496 s->SetMarkerSize(h->GetMarkerSize());
497 s->SetMarkerStyle(h->GetMarkerStyle()+4);
498 s->SetFillColor(h->GetFillColor());
499 s->SetFillStyle(h->GetFillStyle());
500
501 // Find the first and last bin with data
502 Int_t first = nBins+1;
503 Int_t last = 0;
504 for (Int_t i = 1; i <= nBins; i++) {
505 if (h->GetBinContent(i) <= 0) continue;
506 first = TMath::Min(first, i);
507 last = TMath::Max(last, i);
508 }
509
510 Double_t xfirst = h->GetBinCenter(first-1);
511 Int_t f1 = h->GetXaxis()->FindBin(-xfirst);
512 Int_t l2 = s->GetXaxis()->FindBin(xfirst);
513 for (Int_t i = f1, j=l2; i <= last; i++,j--) {
514 s->SetBinContent(j, h->GetBinContent(i));
515 s->SetBinError(j, h->GetBinError(i));
516 }
517 return s;
518 }
519 //__________________________________________________________________
520 /**
521 * Rebin a histogram
522 *
523 * @param h Histogram to rebin
524 * @param rebin Rebinning factor
525 *
526 * @return
527 */
528 virtual void Rebin(TH1* h, Int_t rebin) const
529 {
530 if (rebin <= 1) return;
531
532 Int_t nBins = h->GetNbinsX();
533 if(nBins % rebin != 0) {
534 Warning("Rebin", "Rebin factor %d is not a devisor of current number "
535 "of bins %d in the histogram %s", rebin, nBins, h->GetName());
536 return;
537 }
538
539 // Make a copy
540 TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
541 tmp->Rebin(rebin);
542 tmp->SetDirectory(0);
543
544 // The new number of bins
545 Int_t nBinsNew = nBins / rebin;
546 for(Int_t i = 1;i<= nBinsNew; i++) {
547 Double_t content = 0;
548 Double_t sumw = 0;
549 Double_t wsum = 0;
550 Int_t nbins = 0;
551 for(Int_t j = 1; j<=rebin;j++) {
552 Int_t bin = (i-1)*rebin + j;
553 if(h->GetBinContent(bin) <= 0) continue;
554 Double_t c = h->GetBinContent(bin);
555 Double_t w = 1 / TMath::Power(c,2);
556 content += c;
557 sumw += w;
558 wsum += w * c;
559 nbins++;
560 }
561
562 if(content > 0 ) {
563 tmp->SetBinContent(i, wsum / sumw);
564 tmp->SetBinError(i,TMath::Sqrt(sumw));
565 }
566 }
567
568 // Finally, rebin the histogram, and set new content
569 h->Rebin(rebin);
570 for(Int_t i =1;i<=nBinsNew; i++) {
571 h->SetBinContent(i,tmp->GetBinContent(i));
572 // h->SetBinError(i,tmp->GetBinError(i));
573 }
574
575 delete tmp;
576 }
577};
578
579//____________________________________________________________________
580//
581// EOF
582//
583