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