Various fixes to deal with centrality
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / DrawdNdeta.C
1 /**
2  * Script to visualise the dN/deta 
3  *
4  * This script is independent of any AliROOT code - and should stay that way. 
5  */
6 #include <TH1.h>
7 #include <THStack.h>
8 #include <TGraphErrors.h>
9 #include <TGraphAsymmErrors.h>
10 #include <TMultiGraph.h>
11 #include <TFile.h>
12 #include <TList.h>
13 #include <TString.h>
14 #include <TError.h>
15 #include <TSystem.h>
16 #include <TROOT.h>
17 #include <TMath.h>
18 #include <TCanvas.h>
19 #include <TPad.h>
20 #include <TStyle.h>
21 #include <TLegend.h>
22 #include <TLegendEntry.h>
23 #include <TLatex.h>
24 #include <TImage.h>
25
26 /**
27  * Class to draw dN/deta results 
28  * 
29  */
30 struct dNdetaDrawer 
31 {
32   /**
33    * POD of data for range zooming 
34    */
35   struct RangeParam 
36   {
37     TAxis*       fMasterAxis; // Master axis 
38     TAxis*       fSlave1Axis; // First slave axis 
39     TVirtualPad* fSlave1Pad;  // First slave pad 
40     TAxis*       fSlave2Axis; // Second slave axis 
41     TVirtualPad* fSlave2Pad;  // Second slave pad 
42   };
43   //__________________________________________________________________
44   /** 
45    * Constructor 
46    * 
47    */
48   dNdetaDrawer()
49     : fShowOthers(false),       // Bool_t 
50       fShowAlice(false),        // Bool_t
51       fShowRatios(false),       // Bool_t 
52       fShowLeftRight(false),    // Bool_t
53       fRebin(5),                // UShort_t
54       fCutEdges(false),         // Bool_t
55       fTitle(""),               // TString
56       fTrigString(0),           // TNamed*
57       fSNNString(0),            // TNamed*
58       fSysString(0),            // TNamed*
59       fVtxAxis(0),              // TAxis*
60       // fForward(0),           // TH1*
61       // fForwardMC(0),         // TH1*
62       // fForwardHHD(0),        // TH1*
63       // fTruth(0),             // TH1*
64       // fCentral(0),           // TH1*
65       // fForwardSym(0),        // TH1*
66       // fForwardMCSym(0),      // TH1*
67       // fForwardHHDSym(0),     // TH1*
68       fTriggers(0),             // TH1*
69       fRangeParam(0)
70
71   {
72     fRangeParam = new RangeParam;
73     fRangeParam->fMasterAxis = 0;
74     fRangeParam->fSlave1Axis = 0;
75     fRangeParam->fSlave1Pad  = 0;
76     fRangeParam->fSlave2Axis = 0;
77     fRangeParam->fSlave2Pad  = 0;
78   }
79   //==================================================================
80   /** 
81    * @{ 
82    * @name Set parameters 
83    */
84   /** 
85    * Show other (UA5, CMS, ...) data 
86    * 
87    * @param x Whether to show or not 
88    */
89   void SetShowOthers(Bool_t x)    { fShowOthers = x; }
90   //__________________________________________________________________
91   /** 
92    * Show ALICE published data 
93    * 
94    * @param x Wheter to show or not 
95    */
96   void SetShowAlice(Bool_t x)     { fShowAlice = x; }
97   //__________________________________________________________________
98   /** 
99    * Whether to show ratios or not.  If there's nothing to compare to,
100    * the ratio panel will be implicitly disabled
101    * 
102    * @param x Whether to show or not 
103    */
104   void SetShowRatios(Bool_t x)    { fShowRatios = x; }
105   //__________________________________________________________________
106   /** 
107    * 
108    * Whether to show the left/right asymmetry 
109    *
110    * @param x To show or not 
111    */
112   void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
113   //__________________________________________________________________
114   /** 
115    * Set the rebinning factor 
116    * 
117    * @param x Rebinning factor (must be a divisor in the number of bins) 
118    */
119   void SetRebin(UShort_t x)       { fRebin = x; }
120   //__________________________________________________________________
121   /** 
122    * Wheter to cut away the edges 
123    * 
124    * @param x Whether or not to cut away edges 
125    */
126   void SetCutEdges(Bool_t x)      { fCutEdges = x; }
127   //__________________________________________________________________
128   /** 
129    * Set the title of the plot
130    * 
131    * @param x Title
132    */
133   void SetTitle(TString x)        { fTitle = x; }
134   /* @} */
135   //==================================================================  
136   /** 
137    * @{ 
138    * @name Override settings from input 
139    */
140   /** 
141    * Override setting from file 
142    * 
143    * @param sNN Center of mass energy per nucleon pair (GeV)
144    */
145   void SetSNN(UShort_t sNN) 
146   {
147     fSNNString = new TNamed("sNN", Form("%04dGeV", sNN));
148     fSNNString->SetUniqueID(sNN);
149   }
150   //__________________________________________________________________
151   /** 
152    * Set the collision system 
153    * - 1: pp 
154    * - 2: PbPb
155    * 
156    * @param sys collision system
157    */
158   void SetSys(UShort_t sys)
159   {
160     fSysString = new TNamed("sys", (sys == 1 ? "pp" : 
161                                     sys == 2 ? "PbPb" : "unknown"));
162     fSysString->SetUniqueID(sys);
163   }
164   //__________________________________________________________________
165   /** 
166    * Set the vertex range in centimeters 
167    * 
168    * @param vzMin Min @f$ v_z@f$
169    * @param vzMax Max @f$ v_z@f$
170    */
171   void SetVertexRange(Double_t vzMin, Double_t vzMax) 
172   {
173     fVtxAxis = new TAxis(10, vzMin, vzMax);
174     fVtxAxis->SetName("vtxAxis");
175     fVtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", vzMin, vzMax));
176   }
177   //__________________________________________________________________
178   void SetTrigger(UShort_t trig)
179   {
180     fTrigString = new TNamed("trigString", (trig & 0x1 ? "INEL" : 
181                                             trig & 0x2 ? "INEL>0" : 
182                                             trig & 0x4 ? "NSD" : 
183                                             "unknown"));
184     fTrigString->SetUniqueID(trig);
185   }
186
187
188   //==================================================================
189   /** 
190    * @{ 
191    * @name Main steering functions 
192    */
193   /** 
194    * Run the code to produce the final result. 
195    * 
196    * @param filename  File containing the data 
197    */
198   void Run(const char* filename="forward_dndeta.root") 
199   {
200     Double_t max = 0, rmax=0, amax=0;
201
202     gStyle->SetPalette(1);
203
204     // --- Open input file -------------------------------------------
205     TFile* file = TFile::Open(filename, "READ");
206     if (!file) { 
207       Error("Open", "Cannot open %s", filename);
208       return;
209     }
210     // --- Get forward list ------------------------------------------
211     TList* forward = static_cast<TList*>(file->Get("ForwardResults"));
212     if (!forward) { 
213       Error("Open", "Couldn't find list ForwardResults");
214       return;
215     }
216     // --- Get information on the run --------------------------------
217     FetchInformation(forward);
218     // --- Set the macro pathand load other data script --------------
219     gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
220                              gROOT->GetMacroPath()));
221     gROOT->LoadMacro("OtherData.C");
222
223     // --- Get the central results -----------------------------------
224     TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
225     if (!clusters) Warning("Open", "Couldn't find list CentralResults");
226
227     // --- Make our containtes ---------------------------------------
228     fResults   = new THStack("results", "Results");
229     fRatios    = new THStack("ratios",  "Ratios");
230     fLeftRight = new THStack("asymmetry", "Left-right asymmetry");
231     fOthers    = new TMultiGraph();
232     
233     // --- Loop over input data --------------------------------------
234     FetchResults(forward,  "Forward", max, rmax, amax);
235     FetchResults(clusters, "Central", max, rmax, amax);
236
237     // --- Get trigger information -----------------------------------
238     TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
239     if (sums) {
240       TList* all = static_cast<TList*>(sums->FindObject("all"));
241       if (all) {
242         fTriggers = FetchResult(all, "triggers");
243         if (!fTriggers) all->ls();
244       }
245       else  {
246         Warning("Run", "List all not found in ForwardSums");
247         sums->ls();
248       }
249     }
250     else { 
251       Warning("Run", "No ForwardSums directory found in %s", file->GetName());
252       file->ls();
253     }
254     
255     // --- Check our stacks ------------------------------------------
256     if (!fResults->GetHists() || 
257         fResults->GetHists()->GetEntries() <= 0) { 
258       Error("Run", "No histograms in result stack!");
259       return;
260     }
261     if (!fOthers->GetListOfGraphs() || 
262         fOthers->GetListOfGraphs()->GetEntries() <= 0) { 
263       Warning("Run", "No other data found - disabling that");
264       fShowOthers = false;
265     }
266     if (!fRatios->GetHists() || 
267         fRatios->GetHists()->GetEntries() <= 0) { 
268       Warning("Run", "No ratio data found - disabling that");
269       // fRatios->ls();
270       fShowRatios = false;
271     }
272     if (!fLeftRight->GetHists() || 
273         fLeftRight->GetHists()->GetEntries() <= 0) { 
274       Warning("Run", "No left/right data found - disabling that");
275       // fLeftRight->ls();
276       fShowLeftRight = false;
277     }
278     
279     // --- Close the input file --------------------------------------
280     file->Close();
281
282     
283
284     // --- Plot results ----------------------------------------------
285     Plot(max, rmax, amax);
286   }
287
288   //__________________________________________________________________
289   /** 
290    * Fetch the information on the run from the results list
291    * 
292    * @param results  Results list
293    */
294   void FetchInformation(const TList* results)
295   {
296     if (!fTrigString) 
297       fTrigString = static_cast<TNamed*>(results->FindObject("trigString"));
298     if (!fSNNString) 
299       fSNNString  = static_cast<TNamed*>(results->FindObject("sNN"));
300     if (!fSysString) 
301       fSysString  = static_cast<TNamed*>(results->FindObject("sys"));
302     if (!fVtxAxis)
303       fVtxAxis    = static_cast<TAxis*>(results->FindObject("vtxAxis"));
304     if (!fCentAxis) 
305       fCentAxis   = static_cast<TAxis*>(results->FindObject("centAxis"));
306     if (!fTrigString) fTrigString = new TNamed("trigString", "unknown");
307     if (!fSNNString)  fSNNString  = new TNamed("sNN", "unknown");
308     if (!fSysString)  fSysString  = new TNamed("sys", "unknown");
309     if (!fVtxAxis) { 
310       fVtxAxis    = new TAxis(1,0,0);
311       fVtxAxis->SetName("vtxAxis");
312       fVtxAxis->SetTitle("v_{z} range unspecified");
313     }
314
315     TString centTxt;
316     if (fCentAxis) { 
317       Int_t nCent = fCentAxis->GetNbins();
318       centTxt = Form("\n   Centrality: %d bins", nCent);
319       for (Int_t i = 0; i <= nCent; i++) 
320         centTxt.Append(Form("%c%d", i == 0 ? ' ' : ',', 
321                             fCentAxis->GetXbins()->At(i)));
322     }
323     Info("FetchInformation", 
324          "Initialized for\n"
325          "   Trigger:    %s  (%d)\n"
326          "   sqrt(sNN):  %s  (%dGeV)\n"
327          "   System:     %s  (%d)\n"
328          "   Vz range:   %s  (%f,%f)%s",
329          fTrigString->GetTitle(), fTrigString->GetUniqueID(), 
330          fSNNString->GetTitle(),  fSNNString->GetUniqueID(), 
331          fSysString->GetTitle(),  fSysString->GetUniqueID(), 
332          fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax(),
333          centTxt.Data());
334   }
335   //__________________________________________________________________
336   TMultiGraph* FetchOthers(UShort_t centLow, UShort_t centHigh)
337   {
338     TMultiGraph* thisOther = 0;
339     if (!fShowOthers && !fShowAlice) return 0;
340
341     Bool_t   onlya = (fShowOthers ? false : true);
342     UShort_t sys   = (fSysString  ? fSysString->GetUniqueID() : 0);
343     UShort_t trg   = (fTrigString ? fTrigString->GetUniqueID() : 0);
344     UShort_t snn   = (fSNNString  ? fSNNString->GetUniqueID() : 0);
345     Long_t   ret   = gROOT->ProcessLine(Form("GetData(%d,%d,%d,%d,%d,%d);",
346                                              sys,snn,trg,
347                                              centLow,centHigh,onlya));
348     if (!ret) { 
349       Warning("FetchResults", "No other data found for sys=%d, sNN=%d, "
350               "trigger=%d %d%%-%d%% central %s", 
351               sys, snn, trg, centLow, centHigh, 
352               onlya ? " (ALICE results)" : "all");
353       return 0;
354     }
355     thisOther = reinterpret_cast<TMultiGraph*>(ret);
356     
357     return thisOther;
358   }
359   //__________________________________________________________________
360   /** 
361    * Get the results from the top-level list 
362    * 
363    * @param list  List 
364    * @param name  name 
365    * @param max   On return, maximum of data 
366    * @param rmax  On return, maximum of ratios
367    * @param amax  On return, maximum of left-right comparisons
368    */
369   void FetchResults(const TList* list, 
370                     const char*  name, 
371                     Double_t&    max,
372                     Double_t&    rmax,
373                     Double_t&    amax)
374   {
375     UShort_t n = fCentAxis ? fCentAxis->GetNbins() : 0;
376     if (n == 0) {
377       TList* all = static_cast<TList*>(list->FindObject("all"));
378       if (!all)
379         Error("FetchResults", "Couldn't find list 'all' in %s", 
380               list->GetName());
381       else 
382         FetchResults(all, name, FetchOthers(0,0), -1, 0, max, rmax, amax);
383       return;
384     }
385     
386     Int_t   nCol = gStyle->GetNumberOfColors();
387     for (UShort_t i = 0; i < n; i++) { 
388       UShort_t centLow  = fCentAxis->GetXbins()->At(i);
389       UShort_t centHigh = fCentAxis->GetXbins()->At(i+1);
390       TString  lname    = Form("cent%03d_%03d", centLow, centHigh);
391       TList*   thisCent = static_cast<TList*>(list->FindObject(lname));
392
393       Float_t fc   = (centLow+double(centHigh-centLow)/2) / 100;
394       Int_t   icol = TMath::Min(nCol-1,int(fc * nCol + .5));
395       Int_t   col  = gStyle->GetColorPalette(icol);
396       Info("FetchResults","Centrality %d-%d color index %d (=%d*%f) -> %d", 
397            centLow, centHigh, icol, nCol, fc, col);
398
399       TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
400       if (!thisCent)
401         Error("FetchResults", "Couldn't find list '%s' in %s", 
402               lname.Data(), list->GetName());
403       else 
404         FetchResults(thisCent, name, FetchOthers(centLow, centHigh), 
405                      col, centTxt.Data(), max, rmax, amax);
406     }
407   } 
408   //__________________________________________________________________
409   void SetAttributes(TH1* h, Int_t color)
410   {
411     if (!h) return;
412     if (color < 0) return;
413     h->SetLineColor(color);
414     h->SetMarkerColor(color);
415     h->SetFillColor(color);
416   }
417   //__________________________________________________________________
418   void SetAttributes(TGraph* g, Int_t color)
419   {
420     if (!g) return;
421     if (color < 0) return;
422     g->SetLineColor(color);
423     g->SetMarkerColor(color);
424     g->SetFillColor(color);
425   }
426   //__________________________________________________________________
427   void ModifyTitle(TNamed* h, const char* centTxt)
428   {
429     if (!centTxt || !h) return;
430     h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
431   }
432
433   //__________________________________________________________________
434   /** 
435    * Fetch results for a particular centrality bin
436    * 
437    * @param list       List 
438    * @param name       Name 
439    * @param centLow    Low cut on centrality 
440    * @param centHigh   high cut on centrality 
441    * @param max        On return, data maximum
442    * @param rmax       On return, ratio maximum 
443    * @param amax       On return, left-right maximum 
444    */
445   void FetchResults(const TList* list, 
446                     const char*  name, 
447                     TMultiGraph* thisOther,
448                     Int_t        color,
449                     const char*  centTxt,
450                     Double_t&    max,
451                     Double_t&    rmax,
452                     Double_t&    amax)
453   {
454     TH1* dndeta      = FetchResult(list, Form("dndeta%s", name));
455     TH1* dndetaMC    = FetchResult(list, Form("dndeta%sMC", name));
456     TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
457     TH1* dndetaSym   = 0;
458     TH1* dndetaMCSym = 0;
459     TH1* tracks      = FetchResult(list, "tracks");
460     SetAttributes(dndeta,     color);
461     SetAttributes(dndetaMC,   color+2);
462     SetAttributes(dndetaTruth,color);
463     SetAttributes(dndetaSym,  color);
464     SetAttributes(dndetaMCSym,color+2);
465     SetAttributes(tracks,     color+3);
466     ModifyTitle(dndeta,     centTxt);
467     ModifyTitle(dndetaMC,   centTxt);
468     ModifyTitle(dndetaTruth,centTxt);
469     ModifyTitle(dndetaSym,  centTxt);
470     ModifyTitle(dndetaMCSym,centTxt);
471     ModifyTitle(tracks,     centTxt);
472       
473
474     max = TMath::Max(max, AddHistogram(fResults, dndetaTruth, "e5 p"));
475     max = TMath::Max(max, AddHistogram(fResults, dndetaMC,    dndetaMCSym));
476     max = TMath::Max(max, AddHistogram(fResults, dndeta,      dndetaSym));
477     max = TMath::Max(max, AddHistogram(fResults, tracks));
478     
479     // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f", 
480     //      dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
481
482     if (fShowLeftRight) {
483       fLeftRight->Add(Asymmetry(dndeta,    amax));
484       fLeftRight->Add(Asymmetry(dndetaMC,  amax));
485       fLeftRight->Add(Asymmetry(tracks,    amax));
486     }
487
488     if (thisOther) {
489       TIter next(thisOther->GetListOfGraphs());
490       TGraph* g = 0;
491       while ((g = static_cast<TGraph*>(next()))) {
492         fRatios->Add(Ratio(dndeta,    g, rmax));
493         fRatios->Add(Ratio(dndetaSym, g, rmax));
494         SetAttributes(g, color);
495         ModifyTitle(g, centTxt);
496         if (!fOthers->GetListOfGraphs() || 
497             !fOthers->GetListOfGraphs()->FindObject(g->GetName()))
498           fOthers->Add(g);
499       }
500       // fOthers->Add(thisOther);
501     }
502     if (tracks) { 
503       if (!fRatios->GetHists()->FindObject(tracks->GetName()))
504         fRatios->Add(Ratio(dndeta, tracks, rmax));
505     }
506
507     if (dndetaTruth) { 
508       fRatios->Add(Ratio(dndeta,      dndetaTruth, rmax));
509       fRatios->Add(Ratio(dndetaSym,   dndetaTruth, rmax));
510       fRatios->Add(Ratio(dndetaMC,    dndetaTruth, rmax));
511       fRatios->Add(Ratio(dndetaMCSym, dndetaTruth, rmax));
512     }
513   }
514   //__________________________________________________________________
515   /** 
516    * Plot the results
517    * @param max        Max value 
518    * @param rmax       Maximum diviation from 1 of ratios 
519    * @param amax       Maximum diviation from 1 of asymmetries 
520    */
521   void Plot(Double_t     max, 
522             Double_t     rmax,
523             Double_t     amax)
524   {
525     gStyle->SetOptTitle(0);
526     gStyle->SetTitleFont(132, "xyz");
527     gStyle->SetLabelFont(132, "xyz");
528     
529     Int_t    h  = 800;
530     Int_t    w  = 800; // h / TMath::Sqrt(2);
531     Double_t y1 = 0;
532     Double_t y2 = 0;
533     Double_t y3 = 0;
534     if (!fShowRatios)    w  *= 1.4;
535     else                 y1 =  0.3;
536     if (!fShowLeftRight) w  *= 1.4;
537     else { 
538       Double_t y11 = y1;
539       y1 = (y11 > 0.0001 ? 0.4 : 0.2);
540       y2 = (y11 > 0.0001 ? y11 : 0.2);
541     }
542     TCanvas* c = new TCanvas("Results", "Results", w, h);
543     c->SetFillColor(0);
544     c->SetBorderSize(0);
545     c->SetBorderMode(0);
546
547 #if 0
548     Info("Plot", "y1=%f, y2=%f, y3=%f extra: %s %s", 
549          y1, y2, y2, (fShowRatios ? "ratios" : ""), 
550          (fShowLeftRight ? "right/left" : ""));
551 #endif
552     PlotResults(max, y1);
553     c->cd();
554
555     PlotRatios(rmax, y2, y1);
556     c->cd( );
557
558     PlotLeftRight(amax, y3, y2);
559     c->cd();
560
561     
562     Int_t   vMin = fVtxAxis->GetXmin();
563     Int_t   vMax = fVtxAxis->GetXmax();    
564     TString trg(fTrigString->GetTitle());
565     Int_t   nev  = 0;
566     if (fTriggers) nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
567     trg          = trg.Strip(TString::kBoth);
568     TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev",
569                       fSysString->GetTitle(), 
570                       fSNNString->GetTitle(), 
571                       trg.Data(),
572                       vMin < 0 ? 'm' : 'p',  TMath::Abs(vMin),
573                       vMax < 0 ? 'm' : 'p',  TMath::Abs(vMax),
574                       nev));
575     c->SaveAs(Form("%s.png",  base.Data()));
576     c->SaveAs(Form("%s.root", base.Data()));
577     c->SaveAs(Form("%s.C",    base.Data()));
578   }
579   //__________________________________________________________________
580   /** 
581    * Build main legend 
582    * 
583    * @param x1 
584    * @param y1 
585    * @param x2 
586    * @param y2 
587    */
588   void BuildLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
589   {
590     TLegend* l = new TLegend(x1,y1,x2,y2);
591     l->SetNColumns(2);
592     l->SetFillColor(0);
593     l->SetFillStyle(0);
594     l->SetBorderSize(0);
595     l->SetTextFont(132);
596
597     TIter    next(fResults->GetHists());
598     TObject* hist = 0;
599     while ((hist = next())) { 
600       TString n(hist->GetTitle());
601       if (n.Contains("mirrored")) continue;
602       l->AddEntry(hist, hist->GetTitle(), "pl");
603     }
604     TIter nexto(fOthers->GetListOfGraphs());
605     while ((hist = nexto())) { 
606       TString n(hist->GetTitle());
607       if (n.Contains("mirrored")) continue;
608       l->AddEntry(hist, hist->GetTitle(), "pl");
609     }
610     TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
611     d1->SetLineColor(kBlack);
612     d1->SetMarkerColor(kBlack);
613     d1->SetMarkerStyle(20);
614     TLegendEntry* d2 = l->AddEntry("d2", "Mirrored data", "lp");
615     d2->SetLineColor(kBlack);
616     d2->SetMarkerColor(kBlack);
617     d2->SetMarkerStyle(24);
618     
619     l->Draw();
620   }
621   //__________________________________________________________________
622   /** 
623    * Plot the results
624    *    
625    * @param max       Maximum 
626    * @param yd        Bottom position of pad 
627    */
628   void PlotResults(Double_t max, Double_t yd) 
629   {
630     // Make a sub-pad for the result itself
631     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
632     p1->SetTopMargin(0.05);
633     p1->SetBorderSize(0);
634     p1->SetBorderMode(0);
635     p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
636     p1->SetRightMargin(0.05);
637     p1->SetGridx();
638     p1->SetTicks(1,1);
639     p1->SetNumber(1);
640     p1->Draw();
641     p1->cd();
642     
643     // Info("PlotResults", "Plotting results with max=%f", max);
644     fResults->SetMaximum(1.15*max);
645     fResults->SetMinimum(yd > 0.00001 ? -0.1 : 0);
646
647     FixAxis(fResults, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
648
649     p1->Clear();
650     fResults->DrawClone("nostack e1");
651
652     fRangeParam->fSlave1Axis = fResults->GetXaxis();
653     fRangeParam->fSlave1Pad  = p1;
654
655     // Draw other data
656     if (fShowOthers || fShowAlice) {
657       TGraphAsymmErrors* o      = 0;
658       TIter              nextG(fOthers->GetListOfGraphs());
659       while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
660         o->DrawClone("same p");
661     }
662
663     // Make a legend in the main result pad
664     BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
665 #if 0
666     TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
667     l->SetNColumns(2);
668     l->SetFillColor(0);
669     l->SetFillStyle(0);
670     l->SetBorderSize(0);
671     l->SetTextFont(132);
672 #endif
673
674     // Put a title on top
675     TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
676     tit->SetNDC();
677     tit->SetTextFont(132);
678     tit->SetTextSize(0.05);
679     tit->Draw();
680
681     // Put a nice label in the plot
682     TString     eS;
683     UShort_t    snn = fSNNString->GetUniqueID();
684     const char* sys = fSysString->GetTitle();
685     if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
686     else            eS = Form("%3dGeV", snn);
687     TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s", 
688                                            sys, 
689                                            eS.Data(), 
690                                            fTrigString->GetTitle()));
691     tt->SetNDC();
692     tt->SetTextFont(132);
693     tt->SetTextAlign(33);
694     tt->Draw();
695
696     // Put number of accepted events on the plot
697     Int_t nev = 0;
698     if (fTriggers) nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
699     TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
700     et->SetNDC();
701     et->SetTextFont(132);
702     et->SetTextAlign(33);
703     et->Draw();
704
705     // Put number of accepted events on the plot
706     if (fVtxAxis) { 
707       TLatex* vt = new TLatex(.93, .88, fVtxAxis->GetTitle());
708       vt->SetNDC();
709       vt->SetTextFont(132);
710       vt->SetTextAlign(33);
711       vt->Draw();
712     }
713     // results->Draw("nostack e1 same");
714
715     fRangeParam->fSlave1Axis = FindXAxis(p1, fResults->GetName());
716     fRangeParam->fSlave1Pad  = p1;
717
718
719     // Mark the plot as preliminary
720     TLatex* pt = new TLatex(.12, .93, "Preliminary");
721     pt->SetNDC();
722     pt->SetTextFont(22);
723     pt->SetTextSize(0.07);
724     pt->SetTextColor(kRed+1);
725     pt->SetTextAlign(13);
726     pt->Draw();
727
728     if (!gSystem->AccessPathName("ALICE.png")) { 
729       TPad* logo = new TPad("logo", "logo", .12, .65, .25, .85, 0, 0, 0);
730       logo->SetFillStyle(0);
731       logo->Draw();
732       logo->cd();
733       TImage* i = TImage::Create();
734       i->ReadImage("ALICE.png");
735       i->Draw();
736     }
737     p1->cd();
738   }
739   //__________________________________________________________________
740   /** 
741    * Plot the ratios 
742    * 
743    * @param max     Maximum diviation from 1 
744    * @param y1      Lower y coordinate of pad
745    * @param y2      Upper y coordinate of pad
746    */
747   void PlotRatios(Double_t max, Double_t y1, Double_t y2) 
748   {
749     if (!fShowRatios) return;
750
751     bool isBottom = (y1 < 0.0001);
752     Double_t yd = y2 - y1;
753     // Make a sub-pad for the result itself
754     TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
755     p2->SetTopMargin(0.001);
756     p2->SetRightMargin(0.05);
757     p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
758     p2->SetGridx();
759     p2->SetTicks(1,1);
760     p2->SetNumber(2);
761     p2->Draw();
762     p2->cd();
763
764     // Fix up axis
765     FixAxis(fRatios, 1/yd/1.7, "Ratios", 7);
766
767     fRatios->SetMaximum(1+TMath::Max(.22,1.05*max));
768     fRatios->SetMinimum(1-TMath::Max(.32,1.05*max));
769     p2->Clear();
770     fRatios->DrawClone("nostack e1");
771
772     
773     // Make a legend
774     TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,
775                                   isBottom ? .6 : .4);
776     l2->SetNColumns(2);
777     l2->SetFillColor(0);
778     l2->SetFillStyle(0);
779     l2->SetBorderSize(0);
780     l2->SetTextFont(132);
781
782     // Make a nice band from 0.9 to 1.1
783     TGraphErrors* band = new TGraphErrors(2);
784     band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
785     band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
786     band->SetPointError(0, 0, .1);
787     band->SetPointError(1, 0, .1);
788     band->SetFillColor(kYellow+2);
789     band->SetFillStyle(3002);
790     band->SetLineStyle(2);
791     band->SetLineWidth(1);
792     band->Draw("3 same");
793     band->DrawClone("X L same");
794     
795     // Replot the ratios on top
796     fRatios->DrawClone("nostack e1 same");
797
798     if (isBottom) {
799       fRangeParam->fMasterAxis = FindXAxis(p2, fRatios->GetName());
800       p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)", 
801                                 fRangeParam));
802     }
803     else { 
804       fRangeParam->fSlave2Axis = FindXAxis(p2, fRatios->GetName());
805       fRangeParam->fSlave2Pad  = p2;
806     }
807   }
808   //__________________________________________________________________
809   /** 
810    * Plot the asymmetries
811    * 
812    * @param max     Maximum diviation from 1 
813    * @param y1      Lower y coordinate of pad
814    * @param y2      Upper y coordinate of pad
815    */
816   void PlotLeftRight(Double_t max, Double_t y1, Double_t y2) 
817   {
818     if (!fShowLeftRight) return;
819
820     bool isBottom = (y1 < 0.0001);
821     Double_t yd = y2 - y1;
822     // Make a sub-pad for the result itself
823     TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
824     p3->SetTopMargin(0.001);
825     p3->SetRightMargin(0.05);
826     p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
827     p3->SetGridx();
828     p3->SetTicks(1,1);
829     p3->SetNumber(2);
830     p3->Draw();
831     p3->cd();
832
833     TH1* dummy = 0;
834     if (fLeftRight->GetHists()->GetEntries() == 1) { 
835       // Add dummy histogram
836       dummy = new TH1F("dummy","", 10, -6, 6);
837       dummy->SetLineColor(0);
838       dummy->SetFillColor(0);
839       dummy->SetMarkerColor(0);
840       fLeftRight->Add(dummy);
841     }
842
843     // Fix up axis
844     FixAxis(fLeftRight, 1/yd/1.7, "Right/Left", 4);
845
846     fLeftRight->SetMaximum(1+TMath::Max(.12,1.05*max));
847     fLeftRight->SetMinimum(1-TMath::Max(.15,1.05*max));
848     p3->Clear();
849     fLeftRight->DrawClone("nostack e1");
850
851     
852     // Make a legend
853     TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
854     l2->SetNColumns(2);
855     l2->SetFillColor(0);
856     l2->SetFillStyle(0);
857     l2->SetBorderSize(0);
858     l2->SetTextFont(132);
859 #ifndef __CINT__
860     if (dummy) {
861       TList* prims = l2->GetListOfPrimitives();
862       TIter next(prims);
863       TLegendEntry* o = 0;
864       while ((o = static_cast<TLegendEntry*>(next()))) { 
865         TString lbl(o->GetLabel());
866         if (lbl != "dummy") continue; 
867         prims->Remove(o);
868         break;
869       }
870     }
871 #endif
872     // Make a nice band from 0.9 to 1.1
873     TGraphErrors* band = new TGraphErrors(2);
874     band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
875     band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
876     band->SetPointError(0, 0, .05);
877     band->SetPointError(1, 0, .05);
878     band->SetFillColor(kYellow+2);
879     band->SetFillStyle(3002);
880     band->SetLineStyle(2);
881     band->SetLineWidth(1);
882     band->Draw("3 same");
883     band->DrawClone("X L same");
884
885     fLeftRight->DrawClone("nostack e1 same");
886     if (isBottom) {
887       fRangeParam->fMasterAxis = FindXAxis(p3, fLeftRight->GetName());
888       p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)", 
889                                 fRangeParam));
890     }
891     else { 
892       fRangeParam->fSlave2Axis = FindXAxis(p3, fLeftRight->GetName());
893       fRangeParam->fSlave2Pad  = p3;
894     }
895   }
896   /** @} */
897   //==================================================================
898   /** 
899    * @{ 
900    * @name Data utility functions 
901    */
902   /** 
903    * Get a result from the passed list
904    * 
905    * @param list List to search 
906    * @param name Object name to search for 
907    * 
908    * @return 
909    */
910   TH1* FetchResult(const TList* list, const char* name) const 
911   {
912     if (!list) return 0;
913     
914     TH1* ret = static_cast<TH1*>(list->FindObject(name));
915     if (!ret) {
916       // all->ls();
917       Warning("GetResult", "Histogram %s not found", name);
918     }
919
920     return ret;
921   }
922   //__________________________________________________________________
923   /** 
924    * Add a histogram to the stack after possibly rebinning it  
925    * 
926    * @param stack   Stack to add to 
927    * @param hist    histogram
928    * @param option  Draw options 
929    * 
930    * @return Maximum of histogram 
931    */
932   Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option="") const 
933   {
934     // Check if we have input 
935     if (!hist) return 0;
936
937     // Rebin if needed 
938     Rebin(hist);
939
940     // Info("AddHistogram", "Adding %s to %s", 
941     //      hist->GetName(), stack->GetName());
942     stack->Add(hist, option);
943     // stack->ls();
944     return hist->GetMaximum();
945   }
946   //__________________________________________________________________
947   /** 
948    * Add a histogram to the stack after possibly rebinning it  
949    * 
950    * @param stack   Stack to add to 
951    * @param hist    histogram
952    * @param option  Draw options 
953    * @param sym     On return, the data symmetriced (added to stack)
954    * 
955    * @return Maximum of histogram 
956    */
957   Double_t AddHistogram(THStack* stack, TH1* hist, TH1*& sym,
958                         Option_t* option="") const 
959   {
960     // Check if we have input 
961     if (!hist) return 0;
962
963     // Rebin if needed 
964     Rebin(hist);
965     stack->Add(hist, option);
966
967     // Now symmetrice the histogram 
968     sym = Symmetrice(hist);
969     stack->Add(sym, option);
970
971     // Info("AddHistogram", "Adding %s and %s to %s", 
972     //      hist->GetName(), sym->GetName(), stack->GetName());
973     return hist->GetMaximum();
974   }
975
976   //__________________________________________________________________
977   /** 
978    * Rebin a histogram 
979    * 
980    * @param h     Histogram to rebin
981    * @param rebin Rebinning factor 
982    * 
983    * @return 
984    */
985   virtual void Rebin(TH1* h) const
986   { 
987     if (fRebin <= 1) return;
988
989     Int_t nBins = h->GetNbinsX();
990     if(nBins % fRebin != 0) {
991       Warning("Rebin", "Rebin factor %d is not a devisor of current number "
992               "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
993       return;
994     }
995     
996     // Make a copy 
997     TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
998     tmp->Rebin(fRebin);
999     tmp->SetDirectory(0);
1000     tmp->Reset();
1001     // The new number of bins 
1002     Int_t nBinsNew = nBins / fRebin;
1003     for(Int_t i = 1;i<= nBinsNew; i++) {
1004       Double_t content = 0;
1005       Double_t sumw    = 0;
1006       Double_t wsum    = 0;
1007       Int_t    nbins   = 0;
1008       for(Int_t j = 1; j<=fRebin;j++) {
1009         Int_t    bin = (i-1)*fRebin + j;
1010         Double_t c   =  h->GetBinContent(bin);
1011
1012         if (c <= 0) continue;
1013
1014         if (fCutEdges) {
1015           if (h->GetBinContent(bin+1)<=0 || 
1016               h->GetBinContent(bin-1)) {
1017             Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)", 
1018                     bin, c, h->GetName(), 
1019                     bin+1, h->GetBinContent(bin+1), 
1020                     bin-1, h->GetBinContent(bin-1));
1021             continue;
1022           }     
1023         }
1024         Double_t e =  h->GetBinError(bin);
1025         Double_t w =  1 / (e*e); // 1/c/c
1026         content    += c;
1027         sumw       += w;
1028         wsum       += w * c;
1029         nbins++;
1030       }
1031       
1032       if(content > 0 && nbins > 1 ) {
1033         tmp->SetBinContent(i, wsum / sumw);
1034         tmp->SetBinError(i,1./TMath::Sqrt(sumw));
1035       }
1036     }
1037
1038     // Finally, rebin the histogram, and set new content
1039     h->Rebin(fRebin);
1040     h->Reset();
1041     for(Int_t i = 1; i<= nBinsNew; i++) {
1042       h->SetBinContent(i,tmp->GetBinContent(i));
1043       h->SetBinError(i,  tmp->GetBinError(i));
1044     }
1045     
1046     delete tmp;
1047   }
1048   //__________________________________________________________________
1049   /** 
1050    * Make an extension of @a h to make it symmetric about 0 
1051    * 
1052    * @param h Histogram to symmertrice 
1053    * 
1054    * @return Symmetric extension of @a h 
1055    */
1056   TH1* Symmetrice(const TH1* h) const
1057   {
1058     Int_t nBins = h->GetNbinsX();
1059     TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
1060     s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
1061     s->Reset();
1062     s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
1063     s->SetMarkerColor(h->GetMarkerColor());
1064     s->SetMarkerSize(h->GetMarkerSize());
1065     s->SetMarkerStyle(h->GetMarkerStyle()+4);
1066     s->SetFillColor(h->GetFillColor());
1067     s->SetFillStyle(h->GetFillStyle());
1068     s->SetDirectory(0);
1069
1070     // Find the first and last bin with data 
1071     Int_t first = nBins+1;
1072     Int_t last  = 0;
1073     for (Int_t i = 1; i <= nBins; i++) { 
1074       if (h->GetBinContent(i) <= 0) continue; 
1075       first = TMath::Min(first, i);
1076       last  = TMath::Max(last,  i);
1077     }
1078     
1079     Double_t xfirst = h->GetBinCenter(first-1);
1080     Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
1081     Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
1082     for (Int_t i = f1, j=l2; i <= last; i++,j--) { 
1083       s->SetBinContent(j, h->GetBinContent(i));
1084       s->SetBinError(j, h->GetBinError(i));
1085     }
1086     // Fill in overlap bin 
1087     s->SetBinContent(l2+1, h->GetBinContent(first));
1088     s->SetBinError(l2+1, h->GetBinError(first));
1089     return s;
1090   }
1091   //__________________________________________________________________
1092   /** 
1093    * Calculate the left-right asymmetry of input histogram 
1094    * 
1095    * @param h   Input histogram
1096    * @param max On return, the maximum distance from 1 of the histogram
1097    * 
1098    * @return Asymmetry 
1099    */
1100   TH1* Asymmetry(TH1* h, Double_t& max)
1101   {
1102     if (!h) return 0;
1103
1104     TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
1105     // Int_t    oBins = h->GetNbinsX();
1106     // Double_t high  = h->GetXaxis()->GetXmax();
1107     // Double_t low   = h->GetXaxis()->GetXmin();
1108     // Double_t dBin  = (high - low) / oBins;
1109     // Int_t    tBins = Int_t(2*high/dBin+.5);
1110     // ret->SetBins(tBins, -high, high);
1111     ret->SetDirectory(0);
1112     ret->Reset();
1113     ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
1114     ret->SetYTitle("Right/Left");
1115     Int_t nBins = h->GetNbinsX();
1116     for (Int_t i = 1; i <= nBins; i++) { 
1117       Double_t x = h->GetBinCenter(i);
1118       if (x > 0) break;
1119       
1120       Double_t c1 = h->GetBinContent(i);
1121       Double_t e1 = h->GetBinError(i);
1122       if (c1 <= 0) continue; 
1123       
1124       Int_t    j  = h->FindBin(-x);
1125       if (j <= 0 || j > nBins) continue;
1126
1127       Double_t c2 = h->GetBinContent(j);
1128       Double_t e2 = h->GetBinError(j);
1129
1130       Double_t c12 = c1*c1;
1131       Double_t e   = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
1132       
1133       Int_t    k   = ret->FindBin(x);
1134       ret->SetBinContent(k, c2/c1);
1135       ret->SetBinError(k, e);
1136     }
1137     max = TMath::Max(max, RatioMax(ret));
1138
1139     return ret;
1140   }
1141   //__________________________________________________________________
1142   /** 
1143    * Transform a graph into a histogram 
1144    * 
1145    * @param g 
1146    * 
1147    * @return 
1148    */
1149   TH1* Graph2Hist(const TGraphAsymmErrors* g) const
1150   {
1151     Int_t    nBins = g->GetN();
1152     TArrayF  bins(nBins+1);
1153     Double_t dx = 0;
1154     for (Int_t i = 0; i < nBins; i++) { 
1155       Double_t x   = g->GetX()[i];
1156       Double_t exl = g->GetEXlow()[i];
1157       Double_t exh = g->GetEXhigh()[i];
1158       bins.fArray[i]   = x-exl;
1159       bins.fArray[i+1] = x+exh;
1160       Double_t dxi = exh+exl;
1161       if (i == 0) dx  = dxi;
1162       else if (dxi != dx) dx = 0;
1163     }
1164     TString name(g->GetName());
1165     TString title(g->GetTitle());
1166     TH1D* h = 0;
1167     if (dx != 0) {
1168       h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
1169     }
1170     else {
1171       h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
1172     }
1173     h->SetMarkerStyle(g->GetMarkerStyle());
1174     h->SetMarkerColor(g->GetMarkerColor());
1175     h->SetMarkerSize(g->GetMarkerSize());
1176     
1177     return h;
1178   }
1179   /* @} */
1180   //==================================================================
1181   /** 
1182    * @{ 
1183    * @name Ratio utility functions 
1184    */
1185   /** 
1186    * Get the maximum diviation from 1 in the passed ratio
1187    * 
1188    * @param h Ratio histogram
1189    * 
1190    * @return Max diviation from 1 
1191    */
1192   Double_t RatioMax(TH1* h) const
1193   {
1194     Int_t    nBins = h->GetNbinsX();
1195     Double_t ret   = 0;
1196     for (Int_t i = 1; i <= nBins; i++) { 
1197       Double_t c = h->GetBinContent(i);
1198       if (c == 0) continue;
1199       Double_t e = h->GetBinError(i);
1200       Double_t d = TMath::Abs(1-c-e);
1201       ret        = TMath::Max(d, ret);
1202     }
1203     return ret;
1204   }
1205   //__________________________________________________________________
1206   /** 
1207    * Make the ratio of h1 to h2 
1208    * 
1209    * @param h1 First histogram (numerator) 
1210    * @param h2 Second histogram (denominator)
1211    * 
1212    * @return h1 / h2
1213    */
1214   TH1* Ratio(const TObject* o1, const TObject* o2, Double_t& max) const
1215   {
1216     const TH1* h1 = dynamic_cast<const TH1*>(o1); 
1217     if (h1) { 
1218       const TH1* h2 = dynamic_cast<const TH1*>(o2); 
1219       if (h2) return Ratio(h1,h2,max);
1220       
1221       const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
1222       if (g2) return Ratio(h1,g2,max);      
1223     }
1224     
1225     const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
1226     if (g1) { 
1227       const TGraphAsymmErrors* g2 = dynamic_cast<const TGraphAsymmErrors*>(o2);
1228       if (g2) return Ratio(g1, g2, max);
1229     }
1230
1231     Warning("Ratio", "Don't know how to divide a %s (%s) with a %s (%s)", 
1232             o1->ClassName(), o1->GetName(), o2->ClassName(), o2->GetName());
1233     return 0;
1234   }
1235   //__________________________________________________________________
1236   /** 
1237    * Compute the ratio of @a h to @a g.  @a g is evaluated at the bin
1238    * centers of @a h 
1239    * 
1240    * @param h  Numerator 
1241    * @param g  Divisor 
1242    * 
1243    * @return h/g 
1244    */
1245   TH1* Ratio(const TH1* h, const TGraph* g, Double_t& max) const 
1246   {
1247     if (!h || !g) return 0;
1248
1249     TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
1250     ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
1251     ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
1252     ret->Reset();
1253     ret->SetMarkerStyle(g->GetMarkerStyle());
1254     ret->SetMarkerColor(h->GetMarkerColor());
1255     ret->SetMarkerSize(0.9*g->GetMarkerSize());
1256     ret->SetDirectory(0);
1257     Double_t xlow  = g->GetX()[0];
1258     Double_t xhigh = g->GetX()[g->GetN()-1];
1259     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
1260
1261     for (Int_t i = 1; i <= h->GetNbinsX(); i++) { 
1262       Double_t c = h->GetBinContent(i);
1263       if (c <= 0) continue;
1264
1265       Double_t x = h->GetBinCenter(i);
1266       if (x < xlow || x > xhigh) continue; 
1267
1268       Double_t f = g->Eval(x);
1269       if (f <= 0) continue; 
1270
1271       ret->SetBinContent(i, c / f);
1272       ret->SetBinError(i, h->GetBinError(i) / f);
1273     }
1274     if (ret->GetEntries() <= 0) { 
1275       delete ret; 
1276       ret = 0; 
1277     }
1278     else 
1279       max = TMath::Max(RatioMax(ret), max);
1280     return ret;
1281   }
1282   //__________________________________________________________________
1283   /** 
1284    * Make the ratio of h1 to h2 
1285    * 
1286    * @param h1 First histogram (numerator) 
1287    * @param h2 Second histogram (denominator)
1288    * 
1289    * @return h1 / h2
1290    */
1291   TH1* Ratio(const TH1* h1, const TH1* h2, Double_t& max) const
1292   {
1293     if (!h1 || !h2) return 0;
1294     TH1* t1 = static_cast<TH1*>(h1->Clone(Form("%s_%s", 
1295                                                h1->GetName(), 
1296                                                h2->GetName())));
1297     t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle()));
1298     t1->Divide(h2);
1299     t1->SetMarkerColor(h1->GetMarkerColor());
1300     t1->SetMarkerStyle(h2->GetMarkerStyle());
1301     t1->SetDirectory(0);
1302     max = TMath::Max(RatioMax(t1), max);
1303     return t1;
1304   }
1305   //__________________________________________________________________
1306   /** 
1307    * Calculate the ratio of two graphs - g1 / g2
1308    * 
1309    * @param g1 Numerator 
1310    * @param g2 Denominator
1311    * 
1312    * @return g1 / g2 in a histogram 
1313    */
1314   TH1* Ratio(const TGraphAsymmErrors* g1, 
1315              const TGraphAsymmErrors* g2, Double_t& max) const
1316   {
1317     Int_t    nBins = g1->GetN();
1318     TArrayF  bins(nBins+1);
1319     Double_t dx   = 0;
1320     for (Int_t i = 0; i < nBins; i++) {
1321       Double_t x   = g1->GetX()[i];
1322       Double_t exl = g1->GetEXlow()[i];
1323       Double_t exh = g1->GetEXhigh()[i];
1324       bins.fArray[i]   = x-exl;
1325       bins.fArray[i+1] = x+exh;
1326       Double_t dxi = exh+exl;
1327       if (i == 0) dx  = dxi;
1328       else if (dxi != dx) dx = 0;
1329     }
1330     TString name(Form("%s_%s", g1->GetName(), g2->GetName()));
1331     TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle()));
1332     TH1* h = 0;
1333     if (dx != 0) {
1334       h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
1335     }
1336     else {
1337       h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray);
1338     }
1339     h->SetMarkerStyle(g2->GetMarkerStyle());
1340     h->SetMarkerColor(g1->GetMarkerColor());
1341     h->SetMarkerSize(0.9*g2->GetMarkerSize());
1342     h->SetDirectory(0);
1343
1344     Double_t low  = g2->GetX()[0];
1345     Double_t high = g2->GetX()[g2->GetN()-1];
1346     if (low > high) { Double_t t = low; low = high; high = t; }
1347     for (Int_t i = 0; i < nBins; i++) { 
1348       Double_t x  = g1->GetX()[i];
1349       if (x < low-dx || x > high+dx) continue;
1350       Double_t c1 = g1->GetY()[i];
1351       Double_t e1 = g1->GetErrorY(i);
1352       Double_t c2 = g2->Eval(x);
1353       
1354       h->SetBinContent(i+1, c1 / c2);
1355       h->SetBinError(i+1, e1 / c2);
1356     }
1357     max = TMath::Max(RatioMax(h), max);
1358     return h;
1359   }  
1360   /* @} */
1361   //==================================================================
1362   /** 
1363    * @{ 
1364    * @name Graphics utility functions 
1365    */
1366   /** 
1367    * Find an X axis in a pad 
1368    * 
1369    * @param p     Pad
1370    * @param name  Histogram to find axis for 
1371    * 
1372    * @return Found axis or null
1373    */
1374   TAxis* FindXAxis(TVirtualPad* p, const char* name)
1375   {
1376     TObject* o = p->GetListOfPrimitives()->FindObject(name);
1377     if (!o) { 
1378       Warning("FindXAxis", "%s not found in pad", name);
1379       return 0;
1380     }
1381     THStack* stack = dynamic_cast<THStack*>(o);
1382     if (!stack) { 
1383       Warning("FindXAxis", "%s is not a THStack", name);
1384       return 0;
1385     }
1386     if (!stack->GetHistogram()) { 
1387       Warning("FindXAxis", "%s has no histogram", name);
1388       return 0;
1389     }
1390     TAxis* ret = stack->GetHistogram()->GetXaxis();
1391     return ret;
1392   }
1393
1394   //__________________________________________________________________
1395   /**
1396    * Fix the apperance of the axis in a stack
1397    *
1398    * @param stack  stack of histogram
1399    * @param s      Scaling factor
1400    * @param ytitle Y axis title
1401    * @param force  Whether to draw the stack first or not
1402    * @param ynDiv  Divisions on Y axis
1403    */
1404   void FixAxis(THStack* stack, Double_t s, const char* ytitle,
1405                Int_t ynDiv=210, Bool_t force=true)
1406   {
1407     if (!stack) { 
1408       Warning("FixAxis", "No stack passed for %s!", ytitle);
1409       return;
1410     }
1411     if (force) stack->Draw("nostack e1");
1412
1413     TH1* h = stack->GetHistogram();
1414     if (!h) { 
1415       Warning("FixAxis", "Stack %s has no histogram", stack->GetName());
1416       return;
1417     }
1418     
1419     h->SetXTitle("#eta");
1420     h->SetYTitle(ytitle);
1421     TAxis* xa = h->GetXaxis();
1422     TAxis* ya = h->GetYaxis();
1423     if (xa) {
1424       xa->SetTitle("#eta");
1425       // xa->SetTicks("+-");
1426       xa->SetTitleSize(s*xa->GetTitleSize());
1427       xa->SetLabelSize(s*xa->GetLabelSize());
1428       xa->SetTickLength(s*xa->GetTickLength());
1429
1430       if (stack != fResults) {
1431         TAxis* rxa = fResults->GetXaxis();
1432         xa->Set(rxa->GetNbins(), rxa->GetXmin(), rxa->GetXmax());
1433       }
1434     }
1435     if (ya) {
1436       ya->SetTitle(ytitle);
1437       ya->SetDecimals();
1438       // ya->SetTicks("+-");
1439       ya->SetNdivisions(ynDiv);
1440       ya->SetTitleSize(s*ya->GetTitleSize());
1441       ya->SetTitleOffset(ya->GetTitleOffset()/s);
1442       ya->SetLabelSize(s*ya->GetLabelSize());
1443     }
1444   }
1445   /* @} */
1446
1447
1448
1449   //__________________________________________________________________
1450   Bool_t       fShowOthers;   // Show other data
1451   Bool_t       fShowAlice;    // Show ALICE published data
1452   Bool_t       fShowRatios;   // Show ratios 
1453   Bool_t       fShowLeftRight;// Show asymmetry 
1454   UShort_t     fRebin;        // Rebinning factor 
1455   Bool_t       fCutEdges;     // Whether to cut edges
1456   TString      fTitle;        // Title on plot
1457   TNamed*      fTrigString;   // Trigger string (read, or set)
1458   TNamed*      fSNNString;    // Energy string (read, or set)
1459   TNamed*      fSysString;    // Collision system string (read or set)
1460   TAxis*       fVtxAxis;      // Vertex cuts (read or set)
1461   TAxis*       fCentAxis;     // Centrality axis
1462   THStack*     fResults;      // Stack of results 
1463   THStack*     fRatios;       // Stack of ratios 
1464   THStack*     fLeftRight;    // Left-right asymmetry
1465   TMultiGraph* fOthers;       // Older data 
1466   TH1*         fTriggers;     // Number of triggers
1467   RangeParam*  fRangeParam;   // Parameter object for range zoom 
1468   
1469 };
1470
1471 //=== Stuff for auto zooming =========================================
1472 void UpdateRange(dNdetaDrawer::RangeParam* p)
1473 {
1474   if (!p) { 
1475     Warning("UpdateRange", "No parameters %p", p);
1476     return;
1477   }
1478   if (!p->fMasterAxis) { 
1479     Warning("UpdateRange", "No master axis %p", p->fMasterAxis);
1480     return;
1481   }
1482   Int_t    first = p->fMasterAxis->GetFirst();
1483   Int_t    last  = p->fMasterAxis->GetLast();
1484   Double_t x1    = p->fMasterAxis->GetBinCenter(first);
1485   Double_t x2    = p->fMasterAxis->GetBinCenter(last);
1486   //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2);
1487
1488   if (p->fSlave1Axis) { 
1489     Int_t i1 = p->fSlave1Axis->FindBin(x1);
1490     Int_t i2 = p->fSlave1Axis->FindBin(x2);
1491     p->fSlave1Axis->SetRange(i1, i2);
1492     p->fSlave1Pad->Modified();
1493     p->fSlave1Pad->Update();
1494   }
1495   if (p->fSlave2Axis) { 
1496     Int_t i1 = p->fSlave2Axis->FindBin(x1);
1497     Int_t i2 = p->fSlave2Axis->FindBin(x2);
1498     p->fSlave2Axis->SetRange(i1, i2);
1499     p->fSlave2Pad->Modified();
1500     p->fSlave2Pad->Update();
1501   }
1502   TCanvas*  c = gPad->GetCanvas();
1503   c->cd();
1504 }
1505   
1506 //____________________________________________________________________
1507 void RangeExec(dNdetaDrawer::RangeParam* p)
1508 {
1509   // Event types: 
1510   //  51:   Mouse move 
1511   //  53:   
1512   //  1:    Button down 
1513   //  21:   Mouse drag
1514   //  11:   Mouse release 
1515   // dNdetaDrawer::RangeParam* p = 
1516   //   reinterpret_cast<dNdetaDrawer::RangeParam*>(addr);
1517   Int_t event     = gPad->GetEvent();
1518   TObject *select = gPad->GetSelected();
1519   if (event == 53) { 
1520     UpdateRange(p);
1521     return;
1522   }
1523   if (event != 11 || !select || select != p->fMasterAxis) return;
1524   UpdateRange(p);
1525 }
1526
1527 //=== Steering function ==============================================  
1528 void
1529 DrawdNdeta(const char* filename="forward_dndeta.root", 
1530            Int_t       flags=0xf,
1531            const char* title="",
1532            UShort_t    rebin=5, 
1533            Bool_t      cutEdges=false,
1534            UShort_t    sNN=0, 
1535            UShort_t    sys=0,
1536            UShort_t    trg=1,
1537            Float_t     vzMin=999, 
1538            Float_t     vzMax=-999)
1539 {
1540   dNdetaDrawer d;
1541   d.SetRebin(rebin);
1542   d.SetCutEdges(cutEdges);
1543   d.SetTitle(title);
1544   d.SetShowOthers(flags & 0x1);
1545   d.SetShowAlice(flags & 0x2);
1546   d.SetShowRatios(flags & 0x4);
1547   d.SetShowLeftRight(flags & 0x8);
1548   // Do the below if your input data does not contain these settings 
1549   if (sNN > 0) d.SetSNN(sNN);     // Collision energy per nucleon pair (GeV)
1550   if (sys > 0) d.SetSys(sys);     // Collision system (1:pp, 2:PbPB)
1551   if (trg > 0) d.SetTrigger(trg); // Collision trigger (1:INEL, 2:INEL>0, 4:NSD)
1552   if (vzMin < 999 && vzMax > -999) 
1553     d.SetVertexRange(vzMin,vzMax); // Collision vertex range (cm)
1554   d.Run(filename);
1555 }
1556 //____________________________________________________________________
1557 //
1558 // EOF
1559 //
1560