]>
Commit | Line | Data |
---|---|---|
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 | */ | |
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; | |
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 |