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