8acd19a7f622afa6a451e7f045a9d17c601f63e7
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AnalyseAOD.C
1 /**
2  * @file 
3  * 
4  * @ingroup pwg2_forward_scripts
5  */
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 <TChain.h>
13 #include <TError.h>
14 #include <TSystemDirectory.h>
15 #include <TROOT.h>
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>
24 #include <TSystem.h>
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
35  *   Root> Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/AnalyseAOD.C")
36  *   Root> AnalyseAOD dr
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  * 
44  * @ingroup pwg2_forward_scripts
45  */
46 struct AnalyseAOD 
47 {
48 public: 
49   /** AOD tree */
50   TTree*             fTree;
51   /** AOD object */
52   AliAODForwardMult* fAOD;
53   /** AOD object */
54   AliAODForwardMult* fMCAOD;
55   /** Output file */
56   TFile*             fOut;
57   /** Summed histogram */
58   TH2D*              fSum;
59   /** Summed histogram */
60   TH2D*              fMCSum;
61   /** Primary information */
62   TH2D*              fPrimary;
63   /** Sum primary information */
64   TH2D*              fSumPrimary;
65   /** Vertex efficiency */
66   Double_t           fVtxEff;
67   /** Title to put on the plot */
68   TString            fTitle;
69   /** Do HHD comparison */
70   Bool_t             fDoHHD;
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;
77   /** Number of events accepted */
78   Int_t              fNAll;
79   /** Min vertex */
80   Double_t           fVtxMin;
81   /** Max vertex */
82   Double_t           fVtxMax;
83   
84     
85   //__________________________________________________________________
86   /** 
87    * Constructor 
88    * 
89    */
90   AnalyseAOD()
91     : fTree(0), 
92       fAOD(0),
93       fMCAOD(0),
94       fOut(0), 
95       fSum(0),
96       fMCSum(0),
97       fPrimary(0), 
98       fSumPrimary(0),
99       fVtxEff(0),
100       fTitle(""),
101       fDoHHD(kTRUE)
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     }
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;
127     
128   }
129   //__________________________________________________________________
130   /** 
131    * Run 
132    * 
133    * @param file   Input path for files with AOD trees
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
140    * @param doHHD  Whether to do HHD comparison
141    * @param doComp Whether to do comparisons 
142    * 
143    * @return True on success, false otherwise 
144    */
145   Bool_t Run(const char* file=".", 
146              Double_t vzMin=-10, Double_t vzMax=10, Int_t rebin=1,
147              Int_t mask=AliAODForwardMult::kInel, Int_t energy=900,
148              const char* title="", Bool_t doHHD=false, Bool_t doComp=false)
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());
162     fTitle  = title;
163     fVtxMin = vzMin;
164     fVtxMax = vzMax;
165     if (!Open(file, outName)) return kFALSE;
166     if (!Process(mask)) return kFALSE;
167     if (!Finish(rebin, mask, energy,doHHD,doComp)) return kFALSE;
168
169     return kTRUE;
170   }
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
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
324     // Set the branch pointer 
325     if (fTree->GetBranch("ForwardMC")) 
326       fTree->SetBranchAddress("ForwardMC", &fMCAOD);
327
328     // Set the branch pointer 
329     if (fTree->GetBranch("primary")) 
330       fTree->SetBranchAddress("primary", &fPrimary);
331
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   }
339 #endif
340
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    */
351   Bool_t Process(Int_t mask) 
352   {
353     fNTriggered       = 0;                    // # of triggered ev.
354     fNWithVertex      = 0;                    // # of ev. w/vertex
355     fNAccepted        = 0;                    // # of ev. used
356     Int_t nAvailable  = fTree->GetEntries();  // How many entries
357  
358     for (int i = 0; i < nAvailable; i++) { 
359       fTree->GetEntry(i);
360       if (((i+1) % 100) == 0) {
361         fprintf(stdout,"Event # %9d of %9d, %9d accepted so far\r", 
362                i+1, nAvailable, fNAccepted);
363         fflush(stdout);
364       }
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       }
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       }
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       }
399       
400       // Add contribution from this event 
401       if (fSumPrimary) fSumPrimary->Add(fPrimary);
402      
403       fNAll++;
404       
405       // fAOD->Print();
406       // Other trigger/event requirements could be defined 
407       if (!fAOD->IsTriggerBits(mask)) continue; 
408       fNTriggered++;
409
410       // Check if there's a vertex 
411       if (!fAOD->HasIpZ()) continue; 
412       fNWithVertex++;
413
414       // Select vertex range (in centimeters) 
415       if (!fAOD->InRange(fVtxMin, fVtxMax)) continue; 
416       fNAccepted++;
417  
418       // Add contribution from this event
419       fSum->Add(&(fAOD->GetHistogram()));
420
421       // Add contribution from this event
422       if (fMCSum) fMCSum->Add(&(fMCAOD->GetHistogram()));
423
424       
425     }
426     printf("\n");
427     fVtxEff = Double_t(fNWithVertex)/fNTriggered;
428
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",
433          nAvailable, fNTriggered, fNWithVertex, fNAccepted);
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 
444    * @param doHHD  Whether to do HHD comparison
445    * @param doComp Whether to do comparisons 
446    * 
447    * @return true on success, false otherwise 
448    */
449   Bool_t Finish(Int_t rebin, Int_t mask, Int_t energy, 
450                 Bool_t doHHD, Bool_t doComp)
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");
458     dndeta->SetTitle("ALICE Forward");
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);
467     Rebin(dndeta, rebin);
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
487     TH1D* dndetaTruth = 0;
488     if (fSumPrimary) { 
489       dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth", -1, -1, "e");
490       //dndetaTruth->Scale(1./fNTriggered, "width");
491       dndetaTruth->Scale(1./fNAll, "width");
492       dndetaTruth->SetMarkerColor(kGray+3);
493       dndetaTruth->SetMarkerStyle(22);
494       dndetaTruth->SetMarkerSize(1);
495       dndetaTruth->SetFillStyle(0);
496       Rebin(dndetaTruth, rebin);
497     }
498
499     DrawIt(dndeta, dndetaMC, dndetaTruth, mask, energy, doHHD, doComp);
500
501     return kTRUE;
502   }
503   //__________________________________________________________________
504   /** 
505    */
506   void DrawIt(TH1* dndeta, TH1* dndetaMC, TH1* dndetaTruth, 
507               Int_t mask, Int_t energy, Bool_t doHHD, 
508               Bool_t doComp)
509   {
510     // --- 1st part - prepare data -----------------------------------
511     TH1* dndetaSym = Symmetrice(dndeta);
512
513     Double_t max = dndeta->GetMaximum();
514
515     // Make our histogram stack 
516     THStack* stack = new THStack("results", "Results");
517
518     TH1* dndetaTruthSym = 0;
519     if (dndetaTruth) {
520       dndetaTruth->SetFillColor(kGray);
521       dndetaTruth->SetFillStyle(3001);
522       dndetaTruthSym = Symmetrice(dndetaTruth);
523       stack->Add(dndetaTruthSym, "e5 p");
524       stack->Add(dndetaTruth, "e5 p");
525       Info("DrawIt", "Maximum of MC dndeta (truth): %f, was %f", 
526            dndetaTruth->GetMaximum(),dndeta->GetMaximum());
527       max = TMath::Max(dndetaTruth->GetMaximum(),max);
528     }
529
530     // Get the data from HHD's analysis - if available 
531     TH1* dndetaHHD    = 0;
532     TH1* dndetaHHDSym = 0;
533     // Info("DrawIt", "Will %sdraw HHD results", (doHHD ? "" : "not "));
534     if (doHHD) {
535       TString hhdName(fOut->GetName());
536       hhdName.ReplaceAll("dndeta", "hhd");    
537       dndetaHHD    = GetHHD(hhdName.Data(), mask & AliAODForwardMult::kNSD);
538       dndetaHHDSym = 0;
539       if (dndetaHHD) { 
540         // Symmetrice and add to stack 
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);
546       }
547       else 
548         Warning("DrawIt", "No HHD data found");
549       fOut->cd();
550     }
551
552     // If we have MC analysed data, plot it 
553     if (dndetaMC) { 
554       TH1* dndetaMCSym = Symmetrice(dndetaMC);
555       stack->Add(dndetaMCSym);
556       stack->Add(dndetaMC);
557       max = TMath::Max(dndetaMC->GetMaximum(),max);
558     }
559
560     // Add the analysis results to the list 
561     stack->Add(dndetaSym);
562     stack->Add(dndeta);
563
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.
566     // Info("DrawIt", "Will %sdraw other results", (doComp ? "" : "not "));
567     TMultiGraph* other    = (doComp ? GetData(energy, mask) : 0);
568     THStack*     ratios   = MakeRatios(dndeta,      dndetaSym, 
569                                        dndetaHHD,   dndetaHHDSym, 
570                                        dndetaTruth, dndetaTruthSym, 
571                                        other);
572
573     // Check if we have ratios 
574
575     // --- 2nd part - Draw in canvas ---------------------------------
576     // Make a canvas 
577     gStyle->SetOptTitle(0);
578     gStyle->SetTitleFont(132, "xyz");
579     gStyle->SetLabelFont(132, "xyz");
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
586     Double_t yd = (ratios ? 0.3 : 0);
587
588     // Make a sub-pad for the result itself
589     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0);
590     p1->SetTopMargin(0.05);
591     p1->SetBottomMargin(ratios ? 0.001 : 0.1);
592     p1->SetRightMargin(0.05);
593     p1->SetGridx();
594     p1->SetTicks(1,1);
595     p1->Draw();
596     p1->cd();
597
598     // Fix the apperance of the stack and redraw. 
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());
604         //Info("DrawIt", "Maximum of %s is %f, was %f", o->GetName(),gmax,max);
605         max = TMath::Max(max, gmax);
606       }
607     }
608     max *= 1.1;
609     // Info("DrawIt", "Setting maximum to %f", max);
610     stack->SetMinimum(ratios ? -0.1 : 0);
611     stack->SetMaximum(max);
612     FixAxis(stack, 1/(1-yd)/1.5, 
613             "#frac{1}{#it{N}} #frac{#it{dN}_{ch}}{#it{d#eta}}");
614     p1->Clear();
615     stack->DrawClone("nostack e1");
616
617     // Draw other data 
618     if (other) {
619       TGraphAsymmErrors* o      = 0;
620       TIter              nextG(other->GetListOfGraphs());
621       while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
622         o->Draw("same p");
623     }
624
625     // Make a legend in the main result pad 
626     TString trigs(AliAODForwardMult::GetTriggerString(mask));
627     TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
628     l->SetNColumns(2);
629     l->SetFillColor(0);
630     l->SetFillStyle(0);
631     l->SetBorderSize(0);
632     l->SetTextFont(132);
633     p1->cd();
634
635     // Put a title on top 
636     TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
637     tit->SetNDC();
638     tit->SetTextFont(132);
639     tit->SetTextSize(0.05);
640     tit->Draw();
641
642     // Put a nice label in the plot 
643     TString eS;
644     if (energy > 1000) eS = Form("%4.2fTeV", float(energy)/1000);
645     else               eS = Form("%3dGeV", energy);
646     TLatex* tt = new TLatex(.93, .93, 
647                             Form("#sqrt{s}=%s, %s", eS.Data(),
648                                  AliAODForwardMult::GetTriggerString(mask)));
649     tt->SetNDC();
650     tt->SetTextFont(132);
651     tt->SetTextAlign(33);
652     tt->Draw();
653
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");
671     pt->SetNDC();
672     pt->SetTextFont(22);
673     pt->SetTextSize(0.07);
674     pt->SetTextColor(kRed+1);
675     pt->SetTextAlign(13);
676     pt->Draw();
677     c->cd();
678
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 
688       c->cd();
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();
697
698       // Fix up axis 
699       FixAxis(ratios, 1/yd/1.5, "Ratios", 5);
700
701       // Fix up y range and redraw 
702       ratios->SetMinimum(.58);
703       ratios->SetMaximum(1.22);
704       p2->Clear();
705       ratios->DrawClone("nostack e1");
706       
707       // Make a legend 
708       TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,.6);
709       l2->SetNColumns(2);
710       l2->SetFillColor(0);
711       l2->SetFillStyle(0);
712       l2->SetBorderSize(0);
713       l2->SetTextFont(132);
714
715       // Make a nice band from 0.9 to 1.1 
716       TGraphErrors* band = new TGraphErrors(2);
717       band->SetPoint(0, dndetaSym->GetXaxis()->GetXmin(), 1);
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");
724
725       // Replot the ratios on top 
726       ratios->DrawClone("nostack e1 same");
727
728       c->cd();
729     }
730     
731     // Plot to disk
732     TString imgName(fOut->GetName());
733     imgName.ReplaceAll(".root", ".png");
734     c->SaveAs(imgName.Data());
735     imgName.ReplaceAll(".png", ".C");
736     c->SaveAs(imgName.Data());
737     
738     fOut->cd();
739     stack->Write();
740     if (other)  other->Write();
741     if (ratios) ratios->Write();
742
743     // Close our file 
744     fOut->Close();
745   }
746   //__________________________________________________________________
747   /** 
748    * Get the result from previous analysis code 
749    * 
750    * @param fn  File to open 
751    * @param nsd Whether this is NSD
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;
758     if (gSystem->AccessPathName(fn)) { 
759       Warning("GetHHD", "Output of HHD analysis (%s) not available", fn);
760       return 0;
761     }
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);
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   }
789   //__________________________________________________________________
790   /** 
791    */ 
792   THStack* MakeRatios(const TH1* dndeta, const TH1* sym, 
793                       const TH1* hhd,    const TH1* hhdsym, 
794                       const TH1* mc,     const TH1* mcsym,
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);
801     if (!hasOther && !hasHhd && !mc && !mcsym) return 0;
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())));
821       TH1F* t2 = static_cast<TH1F*>(sym->Clone(Form("%s_%s", 
822                                                     sym->GetName(), 
823                                                     hhdsym->GetName())));
824       t1->SetTitle(Form("%s / %s", dndeta->GetTitle(), hhd->GetTitle()));
825       t2->SetTitle(Form("%s / %s", sym->GetTitle(), hhdsym->GetTitle()));
826       t1->Divide(hhd);
827       t2->Divide(hhdsym);
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());
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));
855 #if 0
856     Info("MakeRatios", "Got a total of %d ratios", !hasRatios ? 0 :
857          ratios->GetHists()->GetEntries());
858 #endif
859
860     if (!hasRatios) { delete ratios; ratios = 0; }
861     return ratios;
862   }
863
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 
872    * @param ynDiv  Divisions on Y axis 
873    */
874   void FixAxis(THStack* stack, Double_t s, const char* ytitle, 
875                Int_t ynDiv=10, Bool_t force=true) 
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);
895       ya->SetDecimals();
896       // ya->SetTicks("+-");
897       ya->SetNdivisions(ynDiv);
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());
922     ret->SetMarkerSize(0.9*g->GetMarkerSize());
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   {
954     fOut->cd();
955
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());
967     // s->SetDirectory(0);
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     }
985     // Fill in overlap bin 
986     s->SetBinContent(l2+1, h->GetBinContent(first));
987     s->SetBinError(l2+1, h->GetBinError(first));
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);
1026         
1027         //Double_t w = 1 / TMath::Power(c,2);
1028         Double_t w = 1 / TMath::Power(h->GetBinError(bin),2);
1029         content    += c;
1030         sumw       += w;
1031         wsum       += w * c;
1032         nbins++;
1033       }
1034       
1035       if(content > 0 ) {
1036         tmp->SetBinContent(i, wsum / sumw);
1037         tmp->SetBinError(i,1/TMath::Sqrt(sumw));
1038       }
1039     }
1040
1041     // Finally, rebin the histogram, and set new content
1042     h->Rebin(rebin);
1043     h->Reset();
1044     for(Int_t i =1;i<=nBinsNew; i++) {
1045       h->SetBinContent(i,tmp->GetBinContent(i));
1046       h->SetBinError(i,tmp->GetBinError(i));
1047     }
1048     
1049     delete tmp;
1050   }
1051 };
1052
1053 //____________________________________________________________________
1054 //
1055 // EOF
1056 //
1057