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