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