Various improvements
[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     if (fTree->GetBranch("ForwardMC"))
280       fTree->SetBranchAddress("ForwardMC", &fMCAOD);
281
282     // Set the branch pointer 
283     if (fTree->GetBranch("primary"))
284       fTree->SetBranchAddress("primary", &fPrimary);
285
286     fOut = TFile::Open(outname, "RECREATE");
287     if (!fOut) { 
288       Error("Open", "Couldn't open %s", outname);
289       return kFALSE;
290     }
291     return kTRUE;
292     
293   }
294
295 #if 0
296   //__________________________________________________________________
297   /** 
298    * Open the files @a fname containg the tree with AliAODEvent objects, 
299    * and also open the output file @a outname for writting. 
300    * 
301    * @param fname   Input file name 
302    * @param outname Output file name 
303    * 
304    * @return true on success, false otherwise 
305    */
306   Bool_t Open(const char* fname, const char* outname) 
307   {
308     Clear("");
309     
310     TFile* file = TFile::Open(fname, "READ");
311     if (!file) { 
312       Error("Open", "Cannot open file %s", fname);
313       return kFALSE;
314     }
315
316     // Get the AOD tree 
317     fTree = static_cast<TTree*>(file->Get("aodTree"));
318     if (!fTree) {
319       Error("Init", "Couldn't get the tree");
320       return kFALSE;
321     }
322
323     // Set the branch pointer 
324     fTree->SetBranchAddress("Forward", &fAOD);
325
326     // Set the branch pointer 
327     if (fTree->GetBranch("ForwardMC")) 
328       fTree->SetBranchAddress("ForwardMC", &fMCAOD);
329
330     // Set the branch pointer 
331     if (fTree->GetBranch("primary")) 
332       fTree->SetBranchAddress("primary", &fPrimary);
333
334     fOut = TFile::Open(outname, "RECREATE");
335     if (!fOut) { 
336       Error("Open", "Couldn't open %s", outname);
337       return kFALSE;
338     }
339     return kTRUE;
340   }
341 #endif
342
343   //__________________________________________________________________
344   /** 
345    * Process the events 
346    * 
347    * @param vzMin  Minimum interaction point z coordinate 
348    * @param vzMax  Maximum interaction point z coordinate 
349    * @param mask   Trigger mask 
350    *
351    * @return true on success, false otherwise 
352    */
353   Bool_t Process(Int_t mask) 
354   {
355     fNTriggered       = 0;                    // # of triggered ev.
356     fNWithVertex      = 0;                    // # of ev. w/vertex
357     fNAccepted        = 0;                    // # of ev. used
358     Int_t nAvailable  = fTree->GetEntries();  // How many entries
359  
360     for (int i = 0; i < nAvailable; i++) { 
361       fTree->GetEntry(i);
362       if (((i+1) % 100) == 0) {
363         fprintf(stdout,"Event # %9d of %9d, %9d accepted so far\r", 
364                i+1, nAvailable, fNAccepted);
365         fflush(stdout);
366       }
367
368       // Create sum histogram on first event - to match binning to input
369       if (!fSum) {
370         fSum = static_cast<TH2D*>(fAOD->GetHistogram().Clone("d2ndetadphi"));
371         Info("Process", "Created sum histogram (%d,%f,%f)x(%d,%f,%f)", 
372              fSum->GetNbinsX(), 
373              fSum->GetXaxis()->GetXmin(),
374              fSum->GetXaxis()->GetXmax(),
375              fSum->GetNbinsY(), 
376              fSum->GetYaxis()->GetXmin(),
377              fSum->GetYaxis()->GetXmax());
378       }
379       if (!fMCSum && fTree->GetBranch("ForwardMC")) { 
380         fMCSum = 
381           static_cast<TH2D*>(fMCAOD->GetHistogram().Clone("d2ndetadphiMC"));
382         Info("Process", "Created MC sum histogram (%d,%f,%f)x(%d,%f,%f)", 
383              fMCSum->GetNbinsX(), 
384              fMCSum->GetXaxis()->GetXmin(),
385              fMCSum->GetXaxis()->GetXmax(),
386              fMCSum->GetNbinsY(), 
387              fMCSum->GetYaxis()->GetXmin(),
388              fMCSum->GetYaxis()->GetXmax());
389       }
390       if (!fSumPrimary && fTree->GetBranch("primary")) { 
391         fSumPrimary = 
392           static_cast<TH2D*>(fPrimary->Clone("primarySum"));
393         Info("Process", "Created MC truth histogram (%d,%f,%f)x(%d,%f,%f)", 
394              fMCSum->GetNbinsX(), 
395              fMCSum->GetXaxis()->GetXmin(),
396              fMCSum->GetXaxis()->GetXmax(),
397              fMCSum->GetNbinsY(), 
398              fMCSum->GetYaxis()->GetXmin(),
399              fMCSum->GetYaxis()->GetXmax());
400       }
401       
402       // Add contribution from this event 
403       if (fSumPrimary) fSumPrimary->Add(fPrimary);
404      
405       fNAll++;
406       
407       // fAOD->Print();
408       // Other trigger/event requirements could be defined 
409       if (!fAOD->IsTriggerBits(mask)) continue; 
410       fNTriggered++;
411
412       // Check if there's a vertex 
413       if (!fAOD->HasIpZ()) continue; 
414       fNWithVertex++;
415
416       // Select vertex range (in centimeters) 
417       if (!fAOD->InRange(fVtxMin, fVtxMax)) continue; 
418       fNAccepted++;
419  
420       // Add contribution from this event
421       fSum->Add(&(fAOD->GetHistogram()));
422
423       // Add contribution from this event
424       if (fMCSum) fMCSum->Add(&(fMCAOD->GetHistogram()));
425
426       
427     }
428     printf("\n");
429     fVtxEff = Double_t(fNWithVertex)/fNTriggered;
430
431     Info("Process", "Total of %9d events\n"
432          "               %9d has a trigger\n" 
433          "               %9d has a vertex\n" 
434          "               %9d was used\n",
435          nAvailable, fNTriggered, fNWithVertex, fNAccepted);
436
437     return kTRUE;
438   }
439   //__________________________________________________________________
440   /** 
441    * Finish the stuff and draw 
442    * 
443    * @param rebin  How many times to rebin
444    * @param energy Collision energy 
445    * @param mask   Trigger mask 
446    * @param doHHD  Whether to do HHD comparison
447    * @param doComp Whether to do comparisons 
448    * 
449    * @return true on success, false otherwise 
450    */
451   Bool_t Finish(Int_t rebin, Int_t mask, Int_t energy, 
452                 Bool_t doHHD, Bool_t doComp)
453   {
454     fOut->cd();
455
456     // Get acceptance normalisation from underflow bins 
457     TH1D* norm   = fSum->ProjectionX("norm", 0, 1, "");
458     // Project onto eta axis - _ignoring_underflow_bins_!
459     TH1D* dndeta = fSum->ProjectionX("dndeta", 1, -1, "e");
460     dndeta->SetTitle("ALICE Forward");
461     // Normalize to the acceptance 
462     dndeta->Divide(norm);
463     // Scale by the vertex efficiency 
464     dndeta->Scale(fVtxEff, "width");
465     dndeta->SetMarkerColor(kRed+1);
466     dndeta->SetMarkerStyle(20);
467     dndeta->SetMarkerSize(1);
468     dndeta->SetFillStyle(0);
469     Rebin(dndeta, rebin);
470
471     TH1D* dndetaMC = 0;
472     if (fMCSum) { 
473       // Get acceptance normalisation from underflow bins 
474       norm = fMCSum->ProjectionX("norm", 0, 1, "");
475       // Project onto eta axis - _ignoring_underflow_bins_!
476       dndetaMC = fMCSum->ProjectionX("dndetaMC", 1, -1, "e");
477       dndetaMC->SetTitle("ALICE Forward (MC)");
478       // Normalize to the acceptance 
479       dndetaMC->Divide(norm);
480       // Scale by the vertex efficiency 
481       dndetaMC->Scale(fVtxEff, "width");
482       dndetaMC->SetMarkerColor(kRed+3);
483       dndetaMC->SetMarkerStyle(21);
484       dndetaMC->SetMarkerSize(1);
485       dndetaMC->SetFillStyle(0);
486       Rebin(dndetaMC, rebin);
487     }
488
489     TH1D* dndetaTruth = 0;
490     if (fSumPrimary) { 
491       dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth", -1, -1, "e");
492       //dndetaTruth->Scale(1./fNTriggered, "width");
493       dndetaTruth->Scale(1./fNAll, "width");
494       dndetaTruth->SetMarkerColor(kGray+3);
495       dndetaTruth->SetMarkerStyle(22);
496       dndetaTruth->SetMarkerSize(1);
497       dndetaTruth->SetFillStyle(0);
498       Rebin(dndetaTruth, rebin);
499     }
500
501     DrawIt(dndeta, dndetaMC, dndetaTruth, mask, energy, doHHD, doComp);
502
503     return kTRUE;
504   }
505   //__________________________________________________________________
506   /** 
507    */
508   void DrawIt(TH1* dndeta, TH1* dndetaMC, TH1* dndetaTruth, 
509               Int_t mask, Int_t energy, Bool_t doHHD, 
510               Bool_t doComp)
511   {
512     // --- 1st part - prepare data -----------------------------------
513     TH1* dndetaSym = Symmetrice(dndeta);
514
515     Double_t max = dndeta->GetMaximum();
516
517     // Make our histogram stack 
518     THStack* stack = new THStack("results", "Results");
519
520     TH1* dndetaTruthSym = 0;
521     if (dndetaTruth) {
522       dndetaTruth->SetFillColor(kGray);
523       dndetaTruth->SetFillStyle(3001);
524       dndetaTruthSym = Symmetrice(dndetaTruth);
525       stack->Add(dndetaTruthSym, "e5 p");
526       stack->Add(dndetaTruth, "e5 p");
527       Info("DrawIt", "Maximum of MC dndeta (truth): %f, was %f", 
528            dndetaTruth->GetMaximum(),dndeta->GetMaximum());
529       max = TMath::Max(dndetaTruth->GetMaximum(),max);
530     }
531
532     // Get the data from HHD's analysis - if available 
533     TH1* dndetaHHD    = 0;
534     TH1* dndetaHHDSym = 0;
535     // Info("DrawIt", "Will %sdraw HHD results", (doHHD ? "" : "not "));
536     if (doHHD) {
537       TString hhdName(fOut->GetName());
538       hhdName.ReplaceAll("dndeta", "hhd");    
539       dndetaHHD    = GetHHD(hhdName.Data(), mask & AliAODForwardMult::kNSD);
540       dndetaHHDSym = 0;
541       if (dndetaHHD) { 
542         // Symmetrice and add to stack 
543         dndetaHHD->SetTitle("ALICE Forward (HHD)");
544         dndetaHHDSym = Symmetrice(dndetaHHD);
545         stack->Add(dndetaHHDSym);
546         stack->Add(dndetaHHD);
547         max = TMath::Max(dndetaHHD->GetMaximum(),max);
548       }
549       else 
550         Warning("DrawIt", "No HHD data found");
551       fOut->cd();
552     }
553
554     // If we have MC analysed data, plot it 
555     if (dndetaMC) { 
556       TH1* dndetaMCSym = Symmetrice(dndetaMC);
557       stack->Add(dndetaMCSym);
558       stack->Add(dndetaMC);
559       max = TMath::Max(dndetaMC->GetMaximum(),max);
560     }
561
562     // Add the analysis results to the list 
563     stack->Add(dndetaSym);
564     stack->Add(dndeta);
565
566     // Get graph of 'other' data - e.g., UA5, CMS, ... - and check if
567     // there's any graphs.  Set the pad division based on that.
568     // Info("DrawIt", "Will %sdraw other results", (doComp ? "" : "not "));
569     TMultiGraph* other    = (doComp ? GetData(energy, mask) : 0);
570     THStack*     ratios   = MakeRatios(dndeta,      dndetaSym, 
571                                        dndetaHHD,   dndetaHHDSym, 
572                                        dndetaTruth, dndetaTruthSym, 
573                                        other);
574
575     // Check if we have ratios 
576
577     // --- 2nd part - Draw in canvas ---------------------------------
578     // Make a canvas 
579     gStyle->SetOptTitle(0);
580     gStyle->SetTitleFont(132, "xyz");
581     gStyle->SetLabelFont(132, "xyz");
582     TCanvas* c = new TCanvas("c", "C", 900, 800);
583     c->SetFillColor(0);
584     c->SetBorderMode(0);
585     c->SetBorderSize(0);
586     c->cd();
587
588     Double_t yd = (ratios ? 0.3 : 0);
589
590     // Make a sub-pad for the result itself
591     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0);
592     p1->SetTopMargin(0.05);
593     p1->SetBottomMargin(ratios ? 0.001 : 0.1);
594     p1->SetRightMargin(0.05);
595     p1->SetGridx();
596     p1->SetTicks(1,1);
597     p1->Draw();
598     p1->cd();
599
600     // Fix the apperance of the stack and redraw. 
601     if (other) {
602       TGraphAsymmErrors* o      = 0;
603       TIter              nextG(other->GetListOfGraphs());
604       while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
605         Double_t gmax = TMath::MaxElement(o->GetN(), o->GetY());
606         //Info("DrawIt", "Maximum of %s is %f, was %f", o->GetName(),gmax,max);
607         max = TMath::Max(max, gmax);
608       }
609     }
610     max *= 1.1;
611     // Info("DrawIt", "Setting maximum to %f", max);
612     stack->SetMinimum(ratios ? -0.1 : 0);
613     stack->SetMaximum(max);
614     FixAxis(stack, 1/(1-yd)/1.5, 
615             "#frac{1}{#it{N}} #frac{#it{dN}_{ch}}{#it{d#eta}}");
616     p1->Clear();
617     stack->DrawClone("nostack e1");
618
619     // Draw other data 
620     if (other) {
621       TGraphAsymmErrors* o      = 0;
622       TIter              nextG(other->GetListOfGraphs());
623       while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
624         o->Draw("same p");
625     }
626
627     // Make a legend in the main result pad 
628     TString trigs(AliAODForwardMult::GetTriggerString(mask));
629     TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
630     l->SetNColumns(2);
631     l->SetFillColor(0);
632     l->SetFillStyle(0);
633     l->SetBorderSize(0);
634     l->SetTextFont(132);
635     p1->cd();
636
637     // Put a title on top 
638     TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
639     tit->SetNDC();
640     tit->SetTextFont(132);
641     tit->SetTextSize(0.05);
642     tit->Draw();
643
644     // Put a nice label in the plot 
645     TString eS;
646     if (energy > 1000) eS = Form("%4.2fTeV", float(energy)/1000);
647     else               eS = Form("%3dGeV", energy);
648     TLatex* tt = new TLatex(.93, .93, 
649                             Form("#sqrt{s}=%s, %s", eS.Data(),
650                                  AliAODForwardMult::GetTriggerString(mask)));
651     tt->SetNDC();
652     tt->SetTextFont(132);
653     tt->SetTextAlign(33);
654     tt->Draw();
655
656     // Put number of accepted events on the plot 
657     TLatex* et = new TLatex(.93, .83, Form("%d events", fNAccepted));
658     et->SetNDC();
659     et->SetTextFont(132);
660     et->SetTextAlign(33);
661     et->Draw();
662
663     // Put number of accepted events on the plot 
664     TLatex* vt = new TLatex(.93, .88, 
665                             Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin, fVtxMax));
666     vt->SetNDC();
667     vt->SetTextFont(132);
668     vt->SetTextAlign(33);
669     vt->Draw();
670
671     // Mark the plot as preliminary
672     TLatex* pt = new TLatex(.12, .93, "Preliminary");
673     pt->SetNDC();
674     pt->SetTextFont(22);
675     pt->SetTextSize(0.07);
676     pt->SetTextColor(kRed+1);
677     pt->SetTextAlign(13);
678     pt->Draw();
679     c->cd();
680
681     // If we do not have the ratios, fix up the display 
682     // p1->SetPad(0, 0, 1, 1);
683     // p1->SetBottomMargin(0.1);
684     // l->SetY1(0.11);
685     // stack->SetMinimum(0);
686     // FixAxis(stack, (1-yd)/1,  "#frac{1}{N} #frac{dN_{ch}}{#eta}",10,false);
687     if (ratios) {
688       // If we do have the ratios, then make a new pad and draw the 
689       // ratios there 
690       c->cd();
691       TPad* p2 = new TPad("p2", "p2", 0, 0.0, 1.0, yd, 0, 0);
692       p2->SetTopMargin(0.001);
693       p2->SetRightMargin(0.05);
694       p2->SetBottomMargin(1/yd * 0.07);
695       p2->SetGridx();
696       p2->SetTicks(1,1);
697       p2->Draw();
698       p2->cd();
699
700       // Fix up axis 
701       FixAxis(ratios, 1/yd/1.5, "Ratios", 5);
702
703       // Fix up y range and redraw 
704       ratios->SetMinimum(.58);
705       ratios->SetMaximum(1.22);
706       p2->Clear();
707       ratios->DrawClone("nostack e1");
708       
709       // Make a legend 
710       TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,.6);
711       l2->SetNColumns(2);
712       l2->SetFillColor(0);
713       l2->SetFillStyle(0);
714       l2->SetBorderSize(0);
715       l2->SetTextFont(132);
716
717       // Make a nice band from 0.9 to 1.1 
718       TGraphErrors* band = new TGraphErrors(2);
719       band->SetPoint(0, dndetaSym->GetXaxis()->GetXmin(), 1);
720       band->SetPoint(1, dndeta->GetXaxis()->GetXmax(), 1);
721       band->SetPointError(0, 0, .1);
722       band->SetPointError(1, 0, .1);
723       band->SetFillColor(kYellow+2);
724       band->SetFillStyle(3002);
725       band->Draw("3 same");
726
727       // Replot the ratios on top 
728       ratios->DrawClone("nostack e1 same");
729
730       c->cd();
731     }
732     
733     // Plot to disk
734     TString imgName(fOut->GetName());
735     imgName.ReplaceAll(".root", ".png");
736     c->SaveAs(imgName.Data());
737     imgName.ReplaceAll(".png", ".C");
738     c->SaveAs(imgName.Data());
739     
740     fOut->cd();
741     stack->Write();
742     if (other)  other->Write();
743     if (ratios) ratios->Write();
744
745     // Close our file 
746     fOut->Close();
747   }
748   //__________________________________________________________________
749   /** 
750    * Get the result from previous analysis code 
751    * 
752    * @param fn  File to open 
753    * @param nsd Whether this is NSD
754    * 
755    * @return null or result of previous analysis code 
756    */
757   TH1* GetHHD(const char* fn="fmd_dNdeta_mult.root", bool nsd=false) 
758   {
759     TDirectory* savdir = gDirectory;
760     if (gSystem->AccessPathName(fn)) { 
761       Warning("GetHHD", "Output of HHD analysis (%s) not available", fn);
762       return 0;
763     }
764     TFile* file = TFile::Open(fn, "READ");
765     if (!file) { 
766       Warning("GetHHD", "couldn't open HHD file %s", fn);
767       savdir->cd();
768       return 0;
769     }
770     TString hist(Form("dNdeta_dNdeta%s", (nsd ? "NSD" : "")));
771     TH1* h = static_cast<TH1*>(file->Get(hist.Data()));
772     if (!h) { 
773       Warning("GetHHD", "Couldn't find HHD histogram %s in %s", 
774               hist.Data(), fn);
775       file->Close();
776       savdir->cd();
777       return 0;
778     }
779     TH1* r = static_cast<TH1*>(h->Clone("dndeta_hhd"));
780     r->SetTitle("1/N dN_{ch}/d#eta (HHD)");
781     r->SetFillStyle(0);
782     r->SetFillColor(0);
783     r->SetMarkerStyle(21);
784     r->SetMarkerColor(kPink+1);
785     r->SetDirectory(savdir);
786
787     file->Close();
788     savdir->cd();
789     return r;
790   }
791   //__________________________________________________________________
792   /** 
793    */ 
794   THStack* MakeRatios(const TH1* dndeta, const TH1* sym, 
795                       const TH1* hhd,    const TH1* hhdsym, 
796                       const TH1* mc,     const TH1* mcsym,
797                       TMultiGraph* other) const 
798   {
799     // If we have 'other' data, then do the ratio of the results to that
800     Bool_t hasOther = (other && other->GetListOfGraphs() && 
801                        other->GetListOfGraphs()->GetEntries() > 0);
802     Bool_t hasHhd   = (hhd && hhdsym);
803     if (!hasOther && !hasHhd && !mc && !mcsym) return 0;
804
805     THStack* ratios = new THStack("ratios", "Ratios");
806     if (hasOther) {
807       TGraphAsymmErrors* o      = 0;
808       TIter              nextG(other->GetListOfGraphs());
809       while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
810         ratios->Add(Ratio(dndeta, o));
811         ratios->Add(Ratio(sym,    o));
812         ratios->Add(Ratio(hhd,    o));
813         ratios->Add(Ratio(hhdsym, o));
814       }
815     }
816
817     // If we have data from HHD's analysis, then do the ratio of 
818     // our result to that data. 
819     if (hasHhd) { 
820       TH1F* t1 = static_cast<TH1F*>(dndeta->Clone(Form("%s_%s", 
821                                                        dndeta->GetName(), 
822                                                        hhd->GetName())));
823       TH1F* t2 = static_cast<TH1F*>(sym->Clone(Form("%s_%s", 
824                                                     sym->GetName(), 
825                                                     hhdsym->GetName())));
826       t1->SetTitle(Form("%s / %s", dndeta->GetTitle(), hhd->GetTitle()));
827       t2->SetTitle(Form("%s / %s", sym->GetTitle(), hhdsym->GetTitle()));
828       t1->Divide(hhd);
829       t2->Divide(hhdsym);
830       t1->SetMarkerColor(hhd->GetMarkerColor());
831       t2->SetMarkerColor(hhdsym->GetMarkerColor());
832       ratios->Add(t1);
833       ratios->Add(t2);
834     }
835
836     // Do comparison to MC 
837     if (mc) { 
838       TH1D* t1 = static_cast<TH1D*>(dndeta->Clone(Form("%s_%s", 
839                                                        dndeta->GetName(), 
840                                                        mc->GetName())));
841       TH1D* t2 = static_cast<TH1D*>(sym->Clone(Form("%s_%s", 
842                                                     sym->GetName(), 
843                                                     mcsym->GetName())));
844       t1->SetTitle(Form("%s / %s", dndeta->GetTitle(), mc->GetTitle()));
845       t2->SetTitle(Form("%s / %s", sym->GetTitle(), mcsym->GetTitle()));
846       t1->Divide(mc);
847       t2->Divide(mcsym);
848       t1->SetMarkerColor(mc->GetMarkerColor());
849       t2->SetMarkerColor(mcsym->GetMarkerColor());
850       ratios->Add(t1);
851       ratios->Add(t2);
852     }
853
854     // Check if we have ratios 
855     Bool_t   hasRatios = (ratios->GetHists() && 
856                           (ratios->GetHists()->GetEntries() > 0));
857 #if 0
858     Info("MakeRatios", "Got a total of %d ratios", !hasRatios ? 0 :
859          ratios->GetHists()->GetEntries());
860 #endif
861
862     if (!hasRatios) { delete ratios; ratios = 0; }
863     return ratios;
864   }
865
866   //__________________________________________________________________
867   /** 
868    * Fix the apperance of the axis in a stack 
869    * 
870    * @param stack  stack of histogram
871    * @param s      Scaling factor 
872    * @param ytitle Y axis title 
873    * @param force  Whether to draw the stack first or not 
874    * @param ynDiv  Divisions on Y axis 
875    */
876   void FixAxis(THStack* stack, Double_t s, const char* ytitle, 
877                Int_t ynDiv=10, Bool_t force=true) 
878   {
879     if (force) stack->Draw("nostack e1");
880
881     TH1* h = stack->GetHistogram();
882     if (!h) return;
883
884     h->SetXTitle("#eta");
885     h->SetYTitle(ytitle);
886     TAxis* xa = h->GetXaxis();
887     TAxis* ya = h->GetYaxis();
888     if (xa) { 
889       xa->SetTitle("#eta");
890       // xa->SetTicks("+-");
891       xa->SetTitleSize(s*xa->GetTitleSize());
892       xa->SetLabelSize(s*xa->GetLabelSize());
893       xa->SetTickLength(s*xa->GetTickLength());
894     }
895     if (ya) { 
896       ya->SetTitle(ytitle);
897       ya->SetDecimals();
898       // ya->SetTicks("+-");
899       ya->SetNdivisions(ynDiv);
900       ya->SetTitleSize(s*ya->GetTitleSize());
901       ya->SetLabelSize(s*ya->GetLabelSize());
902     }      
903   }
904   //__________________________________________________________________
905   /** 
906    * Compute the ratio of @a h to @a g.  @a g is evaluated at the bin
907    * centers of @a h
908    * 
909    * @param h  Numerator 
910    * @param g  Divisor 
911    * 
912    * @return h/g 
913    */
914   TH1* Ratio(const TH1* h, const TGraph* g) const 
915   {
916     if (!h || !g) return 0;
917
918     TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
919     ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
920     ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
921     ret->Reset();
922     ret->SetMarkerStyle(h->GetMarkerStyle());
923     ret->SetMarkerColor(g->GetMarkerColor());
924     ret->SetMarkerSize(0.9*g->GetMarkerSize());
925     Double_t xlow  = g->GetX()[0];
926     Double_t xhigh = g->GetX()[g->GetN()-1];
927     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
928
929     for (Int_t i = 1; i <= h->GetNbinsX(); i++) { 
930       Double_t c = h->GetBinContent(i);
931       if (c <= 0) continue;
932
933       Double_t x = h->GetBinCenter(i);
934       if (x < xlow || x > xhigh) continue; 
935
936       Double_t f = g->Eval(x);
937       if (f <= 0) continue; 
938
939       ret->SetBinContent(i, c / f);
940       ret->SetBinError(i, h->GetBinError(i) / f);
941     }
942     if (ret->GetEntries() <= 0) { delete ret; ret = 0; }
943     return ret;
944   }
945
946   //__________________________________________________________________
947   /** 
948    * Make an extension of @a h to make it symmetric about 0 
949    * 
950    * @param h Histogram to symmertrice 
951    * 
952    * @return Symmetric extension of @a h 
953    */
954   TH1* Symmetrice(const TH1* h) const
955   {
956     fOut->cd();
957
958     Int_t nBins = h->GetNbinsX();
959     TH1*  s     = new TH1D(Form("%s_mirror", h->GetName()),
960                            Form("%s (mirrored)", h->GetTitle()), 
961                            nBins, 
962                            -h->GetXaxis()->GetXmax(), 
963                            -h->GetXaxis()->GetXmin());
964     s->SetMarkerColor(h->GetMarkerColor());
965     s->SetMarkerSize(h->GetMarkerSize());
966     s->SetMarkerStyle(h->GetMarkerStyle()+4);
967     s->SetFillColor(h->GetFillColor());
968     s->SetFillStyle(h->GetFillStyle());
969     // s->SetDirectory(0);
970
971     // Find the first and last bin with data 
972     Int_t first = nBins+1;
973     Int_t last  = 0;
974     for (Int_t i = 1; i <= nBins; i++) { 
975       if (h->GetBinContent(i) <= 0) continue; 
976       first = TMath::Min(first, i);
977       last  = TMath::Max(last,  i);
978     }
979     
980     Double_t xfirst = h->GetBinCenter(first-1);
981     Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
982     Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
983     for (Int_t i = f1, j=l2; i <= last; i++,j--) { 
984       s->SetBinContent(j, h->GetBinContent(i));
985       s->SetBinError(j, h->GetBinError(i));
986     }
987     // Fill in overlap bin 
988     s->SetBinContent(l2+1, h->GetBinContent(first));
989     s->SetBinError(l2+1, h->GetBinError(first));
990     return s;
991   }
992   //__________________________________________________________________
993   /** 
994    * Rebin a histogram 
995    * 
996    * @param h     Histogram to rebin
997    * @param rebin Rebinning factor 
998    * 
999    * @return 
1000    */
1001   virtual void Rebin(TH1* h, Int_t rebin) const
1002   { 
1003     if (rebin <= 1) return;
1004
1005     Int_t nBins = h->GetNbinsX();
1006     if(nBins % rebin != 0) {
1007       Warning("Rebin", "Rebin factor %d is not a devisor of current number "
1008               "of bins %d in the histogram %s", rebin, nBins, h->GetName());
1009       return;
1010     }
1011     
1012     // Make a copy 
1013     TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
1014     tmp->Rebin(rebin);
1015     tmp->SetDirectory(0);
1016
1017     // The new number of bins 
1018     Int_t nBinsNew = nBins / rebin;
1019     for(Int_t i = 1;i<= nBinsNew; i++) {
1020       Double_t content = 0;
1021       Double_t sumw    = 0;
1022       Double_t wsum    = 0;
1023       Int_t    nbins   = 0;
1024       for(Int_t j = 1; j<=rebin;j++) {
1025         Int_t bin = (i-1)*rebin + j;
1026         if(h->GetBinContent(bin) <= 0) continue;
1027         Double_t c =  h->GetBinContent(bin);
1028         Double_t e =  h->GetBinError(bin);
1029         Double_t w =  1 / (e*e); // 1/c/c
1030         content    += c;
1031         sumw       += w;
1032         wsum       += w * c;
1033         nbins++;
1034       }
1035       
1036       if(content > 0 && nbins > 0) {
1037         tmp->SetBinContent(i, wsum / sumw);
1038         tmp->SetBinError(i,1./TMath::Sqrt(sumw));
1039       }
1040     }
1041
1042     // Finally, rebin the histogram, and set new content
1043     h->Rebin(rebin);
1044     h->Reset();
1045     for(Int_t i = 1; i<= nBinsNew; i++) {
1046       h->SetBinContent(i,tmp->GetBinContent(i));
1047       h->SetBinError(i,  tmp->GetBinError(i));
1048     }
1049     
1050     delete tmp;
1051   }
1052 };
1053
1054 //____________________________________________________________________
1055 //
1056 // EOF
1057 //
1058