New method for reading input files (taking fragmented files int account)
[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>
65a1e0cd 16#include <TSystem.h>
7e4038b5 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 */
38struct DrawRes
39{
40public:
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 {
fa4236ed 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
7e4038b5 180
181 for (int i = 0; i < nAvailable; i++) {
182 fTree->GetEntry(i);
65a1e0cd 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 }
7e4038b5 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++;
fa4236ed 205
206 // Check if there's a vertex
207 if (!fAOD->HasIpZ()) continue;
208 nWithVertex++;
209
7e4038b5 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 }
65a1e0cd 217 printf("\n");
fa4236ed 218 fVtxEff = Double_t(nWithVertex)/nTriggered;
7e4038b5 219
fa4236ed 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);
7e4038b5 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 }
fa4236ed 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
7e4038b5 324 // Make a legend
325 TString trigs(AliAODForwardMult::GetTriggerString(mask));
65a1e0cd 326 TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
327 l->SetNColumns(2);
7e4038b5 328 l->SetFillColor(0);
329 l->SetFillStyle(0);
330 l->SetBorderSize(0);
331 p1->cd();
65a1e0cd 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();
7e4038b5 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);
65a1e0cd 352 stack->SetMinimum(0);
353 FixAxis(stack, (1-yd)/1, "#frac{1}{N} #frac{dN_{ch}}{#eta}",10,false);
7e4038b5 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();
65a1e0cd 366 FixAxis(ratios, 1/yd/1.5, "Ratios", 5);
367 ratios->SetMinimum(.58);
368 ratios->SetMaximum(1.22);
7e4038b5 369 p2->Clear();
370 ratios->DrawClone("nostack e1");
371
65a1e0cd 372 TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,.6);
373 l2->SetNColumns(2);
7e4038b5 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;
65a1e0cd 420 if (gSystem->AccessPathName(fn)) {
421 Warning("GetHHD", "Output of HHD analysis (%s) not available", fn);
422 return 0;
423 }
7e4038b5 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,
65a1e0cd 462 Int_t ynDiv=10, Bool_t force=true)
7e4038b5 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);
65a1e0cd 482 ya->SetDecimals();
7e4038b5 483 // ya->SetTicks("+-");
65a1e0cd 484 ya->SetNdivisions(ynDiv);
7e4038b5 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());
65a1e0cd 509 ret->SetMarkerSize(0.9*g->GetMarkerSize());
7e4038b5 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