Mior fixes
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / DrawdNdeta.C
1 /**
2  * @file   DrawdNdeta.C
3  * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
4  * @date   Wed Mar 23 14:07:10 2011
5  * 
6  * @brief  Script to visualise the dN/deta for pp and PbPb
7  *
8  * This script is independent of any AliROOT code - and should stay
9  * that way.
10  * 
11  * The script is <emph>very</emph> long - sigh - the joy of drawing
12  * things nicely in ROOT
13  * 
14  * @ingroup pwg2_forward_dndeta
15  */
16 #include <TH1.h>
17 #include <TColor.h>
18 #include <THStack.h>
19 #include <TGraphErrors.h>
20 #include <TGraphAsymmErrors.h>
21 #include <TMultiGraph.h>
22 #include <TFile.h>
23 #include <TList.h>
24 #include <TString.h>
25 #include <TError.h>
26 #include <TSystem.h>
27 #include <TROOT.h>
28 #include <TMath.h>
29 #include <TCanvas.h>
30 #include <TPad.h>
31 #include <TStyle.h>
32 #include <TLegend.h>
33 #include <TLegendEntry.h>
34 #include <TLatex.h>
35 #include <TImage.h>
36 #include <TRandom.h>
37 #include <fstream>
38 #define SYSERR_COLOR kBlue-10
39 #define SYSERR_STYLE 1001
40
41 Double_t myFunc(Double_t* xp, Double_t* pp);
42
43 /**
44  * Class to draw dN/deta results 
45  * 
46  * @ingroup pwg2_forward_tasks_dndeta
47  * @ingroup pwg2_forward_dndeta
48  */
49 struct dNdetaDrawer 
50 {
51   /**
52    * POD of data for range zooming 
53    */
54   struct RangeParam 
55   {
56     TAxis*       fMasterAxis; // Master axis 
57     TAxis*       fSlave1Axis; // First slave axis 
58     TVirtualPad* fSlave1Pad;  // First slave pad 
59     TAxis*       fSlave2Axis; // Second slave axis 
60     TVirtualPad* fSlave2Pad;  // Second slave pad 
61   };
62   //__________________________________________________________________
63   /** 
64    * Constructor 
65    * 
66    */
67   dNdetaDrawer()
68     : // Options 
69     fShowRatios(false),    // Show ratios 
70     fShowLeftRight(false), // Show asymmetry 
71     fShowRings(false),     // Show rings too
72     fExport(false),        // Export data to script
73     fCutEdges(false),      // Whether to cut edges
74     fRemoveOuters(false),  // Whether to remove outers
75     fShowOthers(0),        // Show other data
76     // Settings 
77     fRebin(0),             // Rebinning factor 
78     fFwdSysErr(0.076),     // Systematic error in forward range
79     fCenSysErr(0),         // Systematic error in central range 
80     fTitle(""),            // Title on plot
81     fClusterScale(""),     // Scaling of clusters to tracklets      
82     // Read (or set) information 
83     fTrigString(0),        // Trigger string (read, or set)
84     fNormString(0),        // Normalisation string (read, or set)
85     fSNNString(0),         // Energy string (read, or set)
86     fSysString(0),         // Collision system string (read or set)
87     fVtxAxis(0),           // Vertex cuts (read or set)
88     fCentAxis(0),          // Centrality axis
89     // Resulting plots 
90     fResults(0),           // Stack of results 
91     fRatios(0),            // Stack of ratios 
92     fLeftRight(0),         // Left-right asymmetry
93     fOthers(0),            // Older data 
94     fTriggers(0),          // Number of triggers
95     fTruth(0),             // Pointer to truth 
96     fRangeParam(0)         // Parameter object for range zoom 
97   {
98     fRangeParam = new RangeParam;
99     fRangeParam->fMasterAxis = 0;
100     fRangeParam->fSlave1Axis = 0;
101     fRangeParam->fSlave1Pad  = 0;
102     fRangeParam->fSlave2Axis = 0;
103     fRangeParam->fSlave2Pad  = 0;
104   }
105   dNdetaDrawer(const dNdetaDrawer&) {}
106   dNdetaDrawer& operator=(const dNdetaDrawer&) { return *this; }
107
108   //__________________________________________________________________
109   virtual ~dNdetaDrawer()
110   {
111     if (fRatios  && fRatios->GetHists())  fRatios->GetHists()->Delete();
112     if (fResults && fResults->GetHists()) fResults->GetHists()->Delete();
113
114     if (fTrigString) { delete fTrigString; fTrigString = 0; }
115     if (fSNNString)  { delete fSNNString;  fSNNString  = 0; }
116     if (fSysString)  { delete fSysString;  fSysString  = 0; }
117     if (fVtxAxis)    { delete fVtxAxis;    fVtxAxis    = 0; }
118     if (fCentAxis)   { delete fCentAxis;   fCentAxis   = 0; }
119     if (fResults)    { delete fResults;    fResults    = 0; }
120     if (fRatios)     { delete fRatios;     fRatios     = 0; }
121     if (fOthers)     { delete fOthers;     fOthers     = 0; }
122     if (fTriggers)   { delete fTriggers;   fTriggers   = 0; } 
123     fRangeParam = 0;
124   }
125
126   //==================================================================
127   /** 
128    * @{ 
129    * @name Set parameters 
130    */
131   /** 
132    * Show other (UA5, CMS, ...) data 
133    * 
134    * @param x Whether to show or not 
135    */
136   void SetShowOthers(UShort_t x)    { fShowOthers = x; }
137   //__________________________________________________________________
138   /** 
139    * Whether to show ratios or not.  If there's nothing to compare to,
140    * the ratio panel will be implicitly disabled
141    * 
142    * @param x Whether to show or not 
143    */
144   void SetShowRatios(Bool_t x)    { fShowRatios = x; }
145   //__________________________________________________________________
146   /** 
147    * 
148    * Whether to show the left/right asymmetry 
149    *
150    * @param x To show or not 
151    */
152   void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
153   //__________________________________________________________________
154   /** 
155    * Whether to show rings 
156    * 
157    * @param x To show or not 
158    */
159   void SetShowRings(Bool_t x) { fShowRings = x; }
160   //__________________________________________________________________
161   /** 
162    * Whether to export results to a script 
163    *
164    * @param x Wheter to export results to a script
165    */
166   void SetExport(Bool_t x)     { fExport = x; }
167   //__________________________________________________________________
168   /** 
169    * Set the rebinning factor 
170    * 
171    * @param x Rebinning factor (must be a divisor in the number of bins) 
172    */
173   void SetRebin(UShort_t x)       { fRebin = x; }
174   //__________________________________________________________________
175   /** 
176    * Wheter to cut away the edges 
177    * 
178    * @param x Whether or not to cut away edges 
179    */
180   void SetCutEdges(Bool_t x)      { fCutEdges = x; }
181   //__________________________________________________________________
182   /** 
183    * Set the title of the plot
184    * 
185    * @param x Title
186    */
187   void SetTitle(TString x)        { fTitle = x; }
188   //__________________________________________________________________
189   /** 
190    * Set the systematic error in the forward region
191    * 
192    * @param e Systematic error in the forward region 
193    */
194   void SetForwardSysError(Double_t e=0) { fFwdSysErr = e; }
195   //__________________________________________________________________
196   /** 
197    * Set the systematic error in the forward region
198    * 
199    * @param e Systematic error in the forward region 
200    */
201   void SetCentralSysError(Double_t e=0) { fCenSysErr = e; }
202   /* @} */
203   //==================================================================  
204   /** 
205    * @{ 
206    * @name Override settings from input 
207    */
208   /** 
209    * Override setting from file 
210    * 
211    * @param sNN Center of mass energy per nucleon pair (GeV)
212    */
213   void SetSNN(UShort_t sNN) 
214   {
215     fSNNString = new TNamed("sNN", Form("%04dGeV", sNN));
216     fSNNString->SetUniqueID(sNN);
217   }
218   //__________________________________________________________________
219   /** 
220    * Set the collision system 
221    * - 1: pp 
222    * - 2: PbPb
223    * 
224    * @param sys collision system
225    */
226   void SetSys(UShort_t sys)
227   {
228     fSysString = new TNamed("sys", (sys == 1 ? "pp" : 
229                                     sys == 2 ? "PbPb" : "unknown"));
230     fSysString->SetUniqueID(sys);
231   }
232   //__________________________________________________________________
233   /** 
234    * Set the vertex range in centimeters 
235    * 
236    * @param vzMin Min @f$ v_z@f$
237    * @param vzMax Max @f$ v_z@f$
238    */
239   void SetVertexRange(Double_t vzMin, Double_t vzMax) 
240   {
241     fVtxAxis = new TAxis(10, vzMin, vzMax);
242     fVtxAxis->SetName("vtxAxis");
243     fVtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", vzMin, vzMax));
244   }
245   //__________________________________________________________________
246   void SetTrigger(UShort_t trig)
247   {
248     fTrigString = new TNamed("trigString", (trig & 0x1 ? "INEL" : 
249                                             trig & 0x2 ? "INEL>0" : 
250                                             trig & 0x4 ? "NSD" : 
251                                             "unknown"));
252     fTrigString->SetUniqueID(trig);
253   }
254
255
256   //==================================================================
257   /** 
258    * @{ 
259    * @name Main steering functions 
260    */
261   /** 
262    * Run the code to produce the final result. 
263    * 
264    * @param filename  File containing the data 
265    */
266   void Run(const char* filename="forward_dndeta.root") 
267   {
268     Double_t max = 0, rmax=0, amax=0;
269
270     gStyle->SetPalette(1);
271
272     // --- Open input file -------------------------------------------
273     TFile* file = TFile::Open(filename, "READ");
274     if (!file) { 
275       Error("Open", "Cannot open %s", filename);
276       return;
277     }
278     // --- Get forward list ------------------------------------------
279     TList* forward = static_cast<TList*>(file->Get("ForwardResults"));
280     if (!forward) { 
281       Error("Open", "Couldn't find list ForwardResults");
282       return;
283     }
284     // --- Get information on the run --------------------------------
285     FetchInformation(forward);
286     // --- Set the macro pathand load other data script --------------
287     gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
288                              gROOT->GetMacroPath()));
289     gROOT->LoadMacro("OtherData.C");
290
291     // --- Get the central results -----------------------------------
292     TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
293     if (!clusters) Warning("Open", "Couldn't find list CentralResults");
294
295     // --- Get the central results -----------------------------------
296     TList* mcTruth = static_cast<TList*>(file->Get("MCTruthResults"));
297     if (!mcTruth) Warning("Open", "Couldn't find list MCTruthResults");
298
299     // --- Make our containtes ---------------------------------------
300     fResults   = new THStack("results", "Results");
301     fRatios    = new THStack("ratios",  "Ratios");
302     fLeftRight = new THStack("asymmetry", "Left-right asymmetry");
303     fOthers    = new TMultiGraph();
304     
305     // --- Loop over input data --------------------------------------
306     /*TObjArray* mcA =*/ FetchResults(mcTruth,  "MCTruth", max, rmax, amax);
307     TObjArray* fwdA = FetchResults(forward,  "Forward", max, rmax, amax);
308     TObjArray* cenA = FetchResults(clusters, "Central", max, rmax, amax);
309
310     // --- Get trigger information -----------------------------------
311     TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
312     if (sums) {
313       TList* all = static_cast<TList*>(sums->FindObject("all"));
314       if (all) {
315         fTriggers = FetchResult(all, "triggers");
316         if (!fTriggers) all->ls();
317       }
318       else  {
319         Warning("Run", "List all not found in ForwardSums");
320         sums->ls();
321       }
322     }
323     else { 
324       Warning("Run", "No ForwardSums directory found in %s", file->GetName());
325       file->ls();
326     }
327     
328     // --- Check our stacks ------------------------------------------
329     if (!fResults->GetHists() || 
330         fResults->GetHists()->GetEntries() <= 0) { 
331       Error("Run", "No histograms in result stack!");
332       return;
333     }
334     if (!fOthers->GetListOfGraphs() || 
335         fOthers->GetListOfGraphs()->GetEntries() <= 0) { 
336       Warning("Run", "No other data found - disabling that");
337       fShowOthers = 0;
338     }
339     if (!fRatios->GetHists() || 
340         fRatios->GetHists()->GetEntries() <= 0) { 
341       Warning("Run", "No ratio data found - disabling that");
342       // fRatios->ls();
343       fShowRatios = false;
344     }
345     if (!fLeftRight->GetHists() || 
346         fLeftRight->GetHists()->GetEntries() <= 0) { 
347       Warning("Run", "No left/right data found - disabling that");
348       // fLeftRight->ls();
349       fShowLeftRight = false;
350     }
351     if (fFwdSysErr > 0) { 
352       if (fCenSysErr <= 0) fCenSysErr = fFwdSysErr;
353       for (Int_t i = 0; i < fwdA->GetEntriesFast(); i++) {
354         TH1* fwd = static_cast<TH1*>(fwdA->At(i));
355         TH1* cen = static_cast<TH1*>(cenA->At(i));
356         CorrectForward(fwd);
357         CorrectCentral(cen);
358         Double_t low, high;
359         TH1* tmp = Merge(cen, fwd, low, high);
360         TF1* f   = FitMerged(tmp, low, high);
361         MakeSysError(tmp, cen, fwd, f);
362         delete f;
363         Info("", "Adding systematic error histogram %s", 
364              tmp->GetName());
365         fResults->GetHists()->AddFirst(tmp, "e5");
366         TH1* tmp2 = Symmetrice(tmp);
367         tmp2->SetFillColor(tmp->GetFillColor());
368         tmp2->SetFillStyle(tmp->GetFillStyle());
369         tmp2->SetMarkerStyle(tmp->GetMarkerStyle());
370         tmp2->SetLineWidth(tmp->GetLineWidth());
371         fResults->GetHists()->AddFirst(tmp2, "e5");
372         fResults->Modified();
373       }
374     }
375     delete fwdA;
376     delete cenA;
377     
378     // --- Close the input file --------------------------------------
379     file->Close();
380
381     
382
383     // --- Plot results ----------------------------------------------
384     Plot(max, rmax, amax);
385   }
386
387   //__________________________________________________________________
388   /** 
389    * Fetch the information on the run from the results list
390    * 
391    * @param results  Results list
392    */
393   void FetchInformation(const TList* results)
394   {
395     if (!fTrigString) 
396       fTrigString = static_cast<TNamed*>(results->FindObject("trigger"));
397     if (!fNormString) 
398       fNormString = static_cast<TNamed*>(results->FindObject("scheme"));
399     if (!fSNNString) 
400       fSNNString  = static_cast<TNamed*>(results->FindObject("sNN"));
401     if (!fSysString) 
402       fSysString  = static_cast<TNamed*>(results->FindObject("sys"));
403     if (!fVtxAxis)
404       fVtxAxis    = static_cast<TAxis*>(results->FindObject("vtxAxis"));
405     if (!fCentAxis) 
406       fCentAxis   = static_cast<TAxis*>(results->FindObject("centAxis"));
407
408     TNamed* options = static_cast<TAxis*>(results->FindObject("options"));
409     if (!fTrigString) fTrigString = new TNamed("trigger", "unknown");
410     if (!fNormString) fNormString = new TNamed("scheme", "unknown");
411     if (!fSNNString)  fSNNString  = new TNamed("sNN", "unknown");
412     if (!fSysString)  fSysString  = new TNamed("sys", "unknown");
413     if (!fVtxAxis) { 
414       fVtxAxis    = new TAxis(1,0,0);
415       fVtxAxis->SetName("vtxAxis");
416       fVtxAxis->SetTitle("v_{z} range unspecified");
417     }
418
419     TString centTxt("none");
420     if (fCentAxis) { 
421       Int_t nCent = fCentAxis->GetNbins();
422       centTxt = Form("%d bins", nCent);
423       for (Int_t i = 0; i <= nCent; i++) 
424         centTxt.Append(Form("%c%d", i == 0 ? ' ' : '-', 
425                             int(fCentAxis->GetXbins()->At(i))));
426     }
427     Info("FetchInformation", 
428          "Initialized for\n"
429          "   Trigger:       %-30s  (%d)\n"
430          "   sqrt(sNN):     %-30s  (%dGeV)\n"
431          "   System:        %-30s  (%d)\n"
432          "   Vz range:      %-30s  (%f,%f)\n"
433          "   Normalization: %-30s  (%d)\n"
434          "   Centrality:    %s\n"
435          "   Options:       %s",
436          fTrigString->GetTitle(), fTrigString->GetUniqueID(), 
437          fSNNString->GetTitle(),  fSNNString->GetUniqueID(), 
438          fSysString->GetTitle(),  fSysString->GetUniqueID(), 
439          fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax(),
440          fNormString->GetTitle(), fNormString->GetUniqueID(),
441          centTxt.Data(), (options ? options->GetTitle() : "none"));
442   }
443   //__________________________________________________________________
444   TMultiGraph* FetchOthers(UShort_t centLow, UShort_t centHigh)
445   {
446     TMultiGraph* thisOther = 0;
447     if (fShowOthers == 0) return 0;
448
449     UShort_t sys   = (fSysString  ? fSysString->GetUniqueID() : 0);
450     UShort_t trg   = (fTrigString ? fTrigString->GetUniqueID() : 0);
451     UShort_t snn   = (fSNNString  ? fSNNString->GetUniqueID() : 0);
452     Long_t   ret   = gROOT->ProcessLine(Form("GetData(%d,%d,%d,%d,%d,%d);",
453                                              sys,snn,trg,
454                                              centLow,centHigh,
455                                              fShowOthers));
456     if (!ret) return 0;
457
458     thisOther = reinterpret_cast<TMultiGraph*>(ret);    
459     return thisOther;
460   }
461   //__________________________________________________________________
462   /** 
463    * Get the results from the top-level list 
464    * 
465    * @param list  List 
466    * @param name  name 
467    * @param max   On return, maximum of data 
468    * @param rmax  On return, maximum of ratios
469    * @param amax  On return, maximum of left-right comparisons
470    */
471   TObjArray* 
472   FetchResults(const TList* list, 
473                const char*  name, 
474                Double_t&    max,
475                Double_t&    rmax,
476                Double_t&    amax)
477   {
478     UShort_t   n = fCentAxis ? fCentAxis->GetNbins() : 0;
479     if (n == 0) {
480       TList* all = static_cast<TList*>(list->FindObject("all"));
481       if (!all) {
482         Error("FetchResults", "Couldn't find list 'all' in %s", 
483               list->GetName());
484         return 0;
485       }
486       TObjArray* a = new TObjArray;
487       TH1*       h = FetchResults(all, name, FetchOthers(0,0), 
488                                   -1000, 0, max, rmax, amax);
489       a->AddAt(h, 0);
490       return a;
491     }
492     
493     TObjArray* a = new TObjArray;
494     for (UShort_t i = 0; i < n; i++) { 
495       UShort_t centLow  = fCentAxis->GetBinLowEdge(i+1);
496       UShort_t centHigh = fCentAxis->GetBinUpEdge(i+1);
497       TString  lname    = Form("cent%03d_%03d", centLow, centHigh);
498       TList*   thisCent = static_cast<TList*>(list->FindObject(lname));
499       Int_t    col      = GetCentralityColor(i+1);
500
501       TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
502       if (!thisCent) {
503         Error("FetchResults", "Couldn't find list '%s' in %s", 
504               lname.Data(), list->GetName());
505         continue;
506       }
507       TH1* h = FetchResults(thisCent, name, FetchOthers(centLow, centHigh), 
508                             col, centTxt.Data(), max, rmax, amax);
509       a->AddAt(h, i);
510     }
511     return a;
512   } 
513   //__________________________________________________________________
514   Int_t GetCentralityColor(Int_t bin) const
515   {
516     UShort_t centLow  = fCentAxis->GetBinLowEdge(bin);
517     UShort_t centHigh = fCentAxis->GetBinUpEdge(bin);
518     Float_t  fc       = (centLow+double(centHigh-centLow)/2) / 100;
519     Int_t    nCol     = gStyle->GetNumberOfColors();
520     Int_t    icol     = TMath::Min(nCol-1,int(fc * nCol + .5));
521     Int_t    col      = gStyle->GetColorPalette(icol);
522     //Info("GetCentralityColor","%3d: %3d-%3d -> %3d",bin,centLow,centHigh,col);
523     return col;
524   }
525   //__________________________________________________________________
526   void SetAttributes(TH1* h, Int_t color)
527   {
528     if (!h) return;
529     if (color < 0) return;
530     // h->SetLineColor(color);
531     h->SetMarkerColor(color);
532     // h->SetFillColor(color);
533   }
534   //__________________________________________________________________
535   void SetAttributes(TGraph* g, Int_t color)
536   {
537     if (!g) return;
538     if (color < 0) return;
539     // g->SetLineColor(color);
540     g->SetMarkerColor(color);
541     // g->SetFillColor(color);
542   }
543   //__________________________________________________________________
544   void ModifyTitle(TNamed* /*h*/, const char* /*centTxt*/)
545   {
546     return;
547     // if (!centTxt || !h) return;
548     // h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
549   }
550
551   //__________________________________________________________________
552   /** 
553    * Fetch results for a particular centrality bin
554    * 
555    * @param list       List 
556    * @param name       Name 
557    * @param thisOther  Other graphs 
558    * @param color      Color 
559    * @param centTxt    Centrality text
560    * @param max        On return, data maximum
561    * @param rmax       On return, ratio maximum 
562    * @param amax       On return, left-right maximum 
563    */
564   TH1* FetchResults(const TList* list, 
565                     const char*  name, 
566                     TMultiGraph* thisOther,
567                     Int_t        color,
568                     const char*  centTxt,
569                     Double_t&    max,
570                     Double_t&    rmax,
571                     Double_t&    amax)
572   {
573     TH1* dndeta      = FetchResult(list, Form("dndeta%s", name));
574     TH1* dndetaMC    = FetchResult(list, Form("dndeta%sMC", name));
575     TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
576     TH1* dndetaSym   = 0;
577     TH1* dndetaMCSym = 0;
578     SetAttributes(dndeta,     color);
579     SetAttributes(dndetaMC,   fCentAxis ? color : color+2);
580     SetAttributes(dndetaTruth,color);
581     SetAttributes(dndetaSym,  color);
582     SetAttributes(dndetaMCSym,fCentAxis ? color : color+2);
583     if (dndetaMC && fCentAxis) 
584       dndetaMC->SetMarkerStyle(dndetaMC->GetMarkerStyle()+2);
585     if (dndetaMCSym && fCentAxis) 
586       dndetaMCSym->SetMarkerStyle(dndetaMCSym->GetMarkerStyle()+2);
587     if (dndetaTruth && fCentAxis) {
588       dndetaTruth->SetMarkerStyle(34);
589       dndetaTruth->SetMarkerColor(kYellow-1);
590     }
591     ModifyTitle(dndeta,     centTxt);
592     ModifyTitle(dndetaMC,   centTxt);
593     ModifyTitle(dndetaTruth,centTxt);
594     ModifyTitle(dndetaSym,  centTxt);
595     ModifyTitle(dndetaMCSym,centTxt);
596       
597
598     max = TMath::Max(max, AddHistogram(fResults, dndetaTruth, "e5 p"));
599     max = TMath::Max(max, AddHistogram(fResults, dndetaMC,    dndetaMCSym));
600     max = TMath::Max(max, AddHistogram(fResults, dndeta,      dndetaSym));
601
602     if (dndetaTruth) {
603       fTruth = dndetaTruth;
604     }
605     else {
606       if (fShowRings) {
607         THStack* rings = static_cast<THStack*>(list->FindObject("dndetaRings"));
608         if (rings) { 
609           TIter next(rings->GetHists());
610           TH1*  hist = 0;
611           while ((hist = static_cast<TH1*>(next()))) 
612             max = TMath::Max(max, AddHistogram(fResults, hist));
613         }
614       }
615       // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f", 
616       //      dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
617       
618       if (fShowLeftRight) {
619         fLeftRight->Add(Asymmetry(dndeta,    amax));
620         fLeftRight->Add(Asymmetry(dndetaMC,  amax));
621       }
622       
623       if (thisOther) {
624         TIter next(thisOther->GetListOfGraphs());
625         TGraph* g = 0;
626         while ((g = static_cast<TGraph*>(next()))) {
627           fRatios->Add(Ratio(dndeta,    g, rmax));
628           fRatios->Add(Ratio(dndetaSym, g, rmax));
629           SetAttributes(g, color);
630           ModifyTitle(g, centTxt);
631           if (!fOthers->GetListOfGraphs() || 
632               !fOthers->GetListOfGraphs()->FindObject(g->GetName())) {
633             max = TMath::Max(max,TMath::MaxElement(g->GetN(), g->GetY()));
634             fOthers->Add(g);
635           }
636         }
637         // fOthers->Add(thisOther);
638       }
639     }
640     if (dndetaMC) { 
641       fRatios->Add(Ratio(dndeta,    dndetaMC,    rmax));
642       fRatios->Add(Ratio(dndetaSym, dndetaMCSym, rmax));
643     }
644     if (fTruth) {
645       fRatios->Add(Ratio(dndeta,      fTruth, rmax));
646       fRatios->Add(Ratio(dndetaSym,   fTruth, rmax));
647     }
648     return dndeta;
649   }
650   //__________________________________________________________________
651   /** 
652    * Plot the results
653    * @param max        Max value 
654    * @param rmax       Maximum diviation from 1 of ratios 
655    * @param amax       Maximum diviation from 1 of asymmetries 
656    */
657   void Plot(Double_t     max, 
658             Double_t     rmax,
659             Double_t     amax)
660   {
661     gStyle->SetOptTitle(0);
662     gStyle->SetTitleFont(132, "xyz");
663     gStyle->SetLabelFont(132, "xyz");
664     
665     Int_t    h  = 800;
666     Int_t    w  = 800; // h / TMath::Sqrt(2);
667     Double_t y1 = 0;
668     Double_t y2 = 0;
669     Double_t y3 = 0;
670     if (!fShowRatios)    w  *= 1.3;
671     else                 y1 =  0.3;
672     if (!fShowLeftRight) w  *= 1.3;
673     else { 
674       Double_t y11 = y1;
675       y1 = (y11 > 0.0001 ? 0.4 : 0.2);
676       y2 = (y11 > 0.0001 ? 0.2 : 0.3);
677     }
678     TCanvas* c = new TCanvas("Results", "Results", w, h);
679     c->SetFillColor(0);
680     c->SetBorderSize(0);
681     c->SetBorderMode(0);
682
683     PlotResults(max, y1);
684     c->cd();
685
686     PlotRatios(rmax, y2, y1);
687     c->cd( );
688
689     PlotLeftRight(amax, y3, y2);
690     c->cd();
691
692     
693     Int_t   vMin = fVtxAxis->GetXmin();
694     Int_t   vMax = fVtxAxis->GetXmax();    
695     TString trg(fTrigString->GetTitle());
696     Int_t   nev  = 0;
697     if (fTriggers) nev = fTriggers->GetBinContent(1);
698     trg          = trg.Strip(TString::kBoth);
699     trg.ReplaceAll(" ", "_");
700     trg.ReplaceAll(">", "Gt");
701     TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev",
702                       fSysString->GetTitle(), 
703                       fSNNString->GetTitle(), 
704                       trg.Data(),
705                       vMin < 0 ? 'm' : 'p',  TMath::Abs(vMin),
706                       vMax < 0 ? 'm' : 'p',  TMath::Abs(vMax),
707                       nev));
708     c->SaveAs(Form("%s.png",  base.Data()));
709     c->SaveAs(Form("%s.root", base.Data()));
710     c->SaveAs(Form("%s.C",    base.Data()));
711     base.ReplaceAll("dndeta", "export");
712     Export(base);
713   }
714   //__________________________________________________________________
715   /** 
716    * Build main legend 
717    * 
718    * @param stack   Stack to include 
719    * @param mg      (optional) Multi graph to include 
720    * @param x1      Lower X coordinate in the range [0,1]
721    * @param y1      Lower Y coordinate in the range [0,1]
722    * @param x2      Upper X coordinate in the range [0,1]
723    * @param y2      Upper Y coordinate in the range [0,1]
724    */
725   void BuildLegend(THStack* stack, TMultiGraph* mg, 
726                    Double_t x1, Double_t y1, Double_t x2, Double_t y2)
727   {
728     TLegend* l = new TLegend(x1,y1,x2,y2);
729     l->SetNColumns(fCentAxis ? 1 : 2);
730     l->SetFillColor(0);
731     l->SetFillStyle(0);
732     l->SetBorderSize(0);
733     l->SetTextFont(132);
734
735     // Loop over items in stack and get unique items, while ignoring
736     // mirrored data and systematic error bands 
737     TIter    next(stack->GetHists());
738     TH1*     hist = 0;
739     TObjArray unique;
740     unique.SetOwner();
741     Bool_t   sysErrSeen = false;
742     while ((hist = static_cast<TH1*>(next()))) { 
743       TString t(hist->GetTitle());
744       TString n(hist->GetName());
745       n.ToLower();
746       if (t.Contains("mirrored")) continue;
747       if (n.Contains("syserror")) { sysErrSeen = true; continue; }
748       if (unique.FindObject(t.Data())) continue;
749       TObjString* s = new TObjString(hist->GetTitle());
750       s->SetUniqueID(((hist->GetMarkerStyle() & 0xFFFF) << 16) |
751                      ((hist->GetMarkerColor() & 0xFFFF) <<  0));
752       unique.Add(s);
753       // l->AddEntry(hist, hist->GetTitle(), "pl");
754     }
755     if (mg) {
756       // If we have other stuff, scan for unique names 
757       TIter nexto(mg->GetListOfGraphs());
758       TGraph* g = 0;
759       while ((g = static_cast<TGraph*>(nexto()))) { 
760         TString n(g->GetTitle());
761         if (n.Contains("mirrored")) continue;
762         if (unique.FindObject(n.Data())) continue;
763         TObjString* s = new TObjString(n);
764         s->SetUniqueID(((g->GetMarkerStyle() & 0xFFFF) << 16) |
765                        ((g->GetMarkerColor() & 0xFFFF) <<  0));
766         unique.Add(s);
767         // l->AddEntry(hist, hist->GetTitle(), "pl");
768       }
769     }
770
771     // Add legend entries for unique items only
772     TIter nextu(&unique);
773     TObject* s = 0;
774     Int_t i = 0;
775     while ((s = nextu())) { 
776       TLegendEntry* dd = l->AddEntry(Form("data%2d", i++), 
777                                      s->GetName(), "lp");
778       Int_t style = (s->GetUniqueID() >> 16) & 0xFFFF;
779       Int_t color = (s->GetUniqueID() >>  0) & 0xFFFF;
780       dd->SetLineColor(kBlack);
781       if (fCentAxis) dd->SetMarkerColor(kBlack);
782       else           dd->SetMarkerColor(color);
783       dd->SetMarkerStyle(style);
784     }
785     if (sysErrSeen) {
786       // Add entry for systematic errors 
787       TLegendEntry* d0 = l->AddEntry("d0", Form("%4.1f%% Systematic error", 
788                                                 100*fFwdSysErr), "f");
789       d0->SetLineColor(SYSERR_COLOR);
790       d0->SetMarkerColor(SYSERR_COLOR);
791       d0->SetFillColor(SYSERR_COLOR);
792       d0->SetFillStyle(SYSERR_STYLE);
793       d0->SetMarkerStyle(0);
794       d0->SetLineWidth(0);
795       i++;
796     }
797     if (!fCentAxis && i % 2 == 1)  {
798       // To make sure the 'data' and 'mirrored' entries are on a line
799       // by themselves 
800       TLegendEntry* dd = l->AddEntry("dd", "   ", "");
801       dd->SetTextSize(0);
802       dd->SetFillStyle(0);
803       dd->SetFillColor(0);
804       dd->SetLineWidth(0);
805       dd->SetLineColor(0);
806       dd->SetMarkerSize(0);
807     }
808     // Add entry for 'data'
809     TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
810     d1->SetLineColor(kBlack);
811     d1->SetMarkerColor(kBlack);
812     d1->SetMarkerStyle(20);
813     // Add entry for 'mirrored data'
814     TLegendEntry* d2 = l->AddEntry("d2", "Mirrored data", "lp");
815     d2->SetLineColor(kBlack);
816     d2->SetMarkerColor(kBlack);
817     d2->SetMarkerStyle(24);
818     
819     l->Draw();
820   }
821   //__________________________________________________________________
822   /** 
823    * Build centrality legend 
824    * 
825    * @param axis    Stack to include 
826    * @param x1      Lower X coordinate in the range [0,1]
827    * @param y1      Lower Y coordinate in the range [0,1]
828    * @param x2      Upper X coordinate in the range [0,1]
829    * @param y2      Upper Y coordinate in the range [0,1]
830    */
831   void BuildCentLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
832   {
833     if (!fCentAxis) return;
834
835     TLegend* l = new TLegend(x1,y1,x2,y2);
836     l->SetNColumns(1);
837     l->SetFillColor(0);
838     l->SetFillStyle(0);
839     l->SetBorderSize(0);
840     l->SetTextFont(132);
841
842     Int_t n = fCentAxis->GetNbins();
843     for (Int_t i = 1; i <= n; i++) { 
844       Double_t low = fCentAxis->GetBinLowEdge(i);
845       Double_t upp = fCentAxis->GetBinUpEdge(i);
846       TLegendEntry* e = l->AddEntry(Form("dummy%02d", i),
847                                     Form("%3d%% - %3d%%", 
848                                          int(low), int(upp)), "pl");
849       e->SetMarkerColor(GetCentralityColor(i));
850     }
851     l->Draw();
852   }
853   //__________________________________________________________________
854   /** 
855    * Plot the results
856    *    
857    * @param max       Maximum 
858    * @param yd        Bottom position of pad 
859    */
860   void PlotResults(Double_t max, Double_t yd) 
861   {
862     // Make a sub-pad for the result itself
863     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
864     p1->SetTopMargin(0.05);
865     p1->SetBorderSize(0);
866     p1->SetBorderMode(0);
867     p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
868     p1->SetRightMargin(0.03);
869     if (fShowLeftRight || fShowRatios) p1->SetGridx();
870     p1->SetTicks(1,1);
871     p1->SetNumber(1);
872     p1->Draw();
873     p1->cd();
874
875     // Info("PlotResults", "Plotting results with max=%f", max);
876     fResults->SetMaximum(1.15*max);
877     fResults->SetMinimum(yd > 0.00001 ? -0.1 : 0);
878
879     FixAxis(fResults, (1-yd)*(yd > .001 ? 1 : .9 / 1.2), 
880             "#font[12]{#frac{1}{N} "
881             "#frac{dN_{#font[132]{ch}}}{d#font[152]{#eta}}}");
882
883     p1->Clear();
884     fResults->DrawClone("nostack e1");
885
886     fRangeParam->fSlave1Axis = fResults->GetXaxis();
887     fRangeParam->fSlave1Pad  = p1;
888
889     // Draw other data
890     if (fShowOthers != 0) {
891       TGraphAsymmErrors* o      = 0;
892       TIter              nextG(fOthers->GetListOfGraphs());
893       while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
894         o->DrawClone("same p");
895     }
896
897     // Make a legend in the main result pad
898     BuildCentLegend(.12, 1-p1->GetTopMargin()-.01-.5,  
899                     .35, 1-p1->GetTopMargin()-.01-.1);
900     Double_t x1 = (fCentAxis ? .7 : .15); 
901     Double_t x2 = (fCentAxis ? 1-p1->GetRightMargin()-.01: .90);
902     Double_t y1 = (fCentAxis ? .5: p1->GetBottomMargin()+.01); 
903     Double_t y2 = (fCentAxis ? 1-p1->GetTopMargin()-.01-.15 : .35);
904                    
905     BuildLegend(fResults, fOthers, x1, y1, x2, y2);
906
907     // Put a title on top
908     fTitle.ReplaceAll("@", " ");
909     TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
910     tit->SetNDC();
911     tit->SetTextFont(132);
912     tit->SetTextSize(0.05);
913     tit->Draw();
914
915     Int_t aliceBlue = TColor::GetColor(41,73,156);
916     // Put a nice label in the plot
917     TString     eS;
918     UShort_t    snn = fSNNString->GetUniqueID();
919     const char* sys = fSysString->GetTitle();
920     if (snn == 2750) snn = 2760;
921     if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
922     else            eS = Form("%3dGeV", snn);
923     TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s%s}=%s, %s", 
924                                            sys, 
925                                            (fCentAxis ? "_{NN}" : ""),
926                                            eS.Data(), 
927                                            fCentAxis ? "by centrality" : 
928                                            fTrigString->GetTitle()));
929     tt->SetTextColor(aliceBlue);
930     tt->SetNDC();
931     tt->SetTextFont(132);
932     tt->SetTextAlign(33);
933     tt->Draw();
934
935     // Put number of accepted events on the plot
936     Int_t nev = 0;
937     if (fTriggers) nev = fTriggers->GetBinContent(1);
938     TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
939     et->SetTextColor(aliceBlue);
940     et->SetNDC();
941     et->SetTextFont(132);
942     et->SetTextAlign(33);
943     et->Draw();
944
945     // Put number of accepted events on the plot
946     if (fVtxAxis) { 
947       TLatex* vt = new TLatex(.93, .88, fVtxAxis->GetTitle());
948       vt->SetNDC();
949       vt->SetTextFont(132);
950       vt->SetTextAlign(33);
951       vt->SetTextColor(aliceBlue);
952       vt->Draw();
953     }
954     // results->Draw("nostack e1 same");
955
956     fRangeParam->fSlave1Axis = FindXAxis(p1, fResults->GetName());
957     fRangeParam->fSlave1Pad  = p1;
958
959
960     // Mark the plot as preliminary
961     TLatex* pt = new TLatex(.12, .93, "Work in progress");
962     pt->SetNDC();
963     pt->SetTextFont(22);
964     // pt->SetTextSize();
965     pt->SetTextColor(TColor::GetColor(234,26,46));
966     pt->SetTextAlign(13);
967     pt->Draw();
968
969     const char*  logos[] = { "ALICE.png", "FMD.png", 0 };
970     const char** logo    = logos;
971     while (*logo) {
972       if (gSystem->AccessPathName(*logo)) {
973         logo++;
974         continue;
975       }
976       TPad* pad = new TPad("logo", "logo", .12, .7, .25, .9, 0, 0, 0);
977       pad->SetFillStyle(0);
978       pad->Draw();
979       pad->cd();
980       TImage* i = TImage::Create();
981       i->ReadImage(*logo);
982       i->Draw();
983       break;
984     }
985     p1->cd();
986   }
987   //__________________________________________________________________
988   /** 
989    * Plot the ratios 
990    * 
991    * @param max     Maximum diviation from 1 
992    * @param y1      Lower y coordinate of pad
993    * @param y2      Upper y coordinate of pad
994    */
995   void PlotRatios(Double_t max, Double_t y1, Double_t y2) 
996   {
997     if (fShowRatios == 0) return;
998
999     bool isBottom = (y1 < 0.0001);
1000     Double_t yd = y2 - y1;
1001     // Make a sub-pad for the result itself
1002     TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
1003     p2->SetTopMargin(0.001);
1004     p2->SetRightMargin(0.03);
1005     p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
1006     p2->SetGridx();
1007     p2->SetTicks(1,1);
1008     p2->SetNumber(2);
1009     p2->Draw();
1010     p2->cd();
1011
1012     // Fix up axis
1013     FixAxis(fRatios, yd, "Ratios", 7);
1014
1015     fRatios->SetMaximum(1+TMath::Max(.22,1.05*max));
1016     fRatios->SetMinimum(1-TMath::Max(.32,1.05*max));
1017     p2->Clear();
1018     fRatios->DrawClone("nostack e1");
1019
1020     
1021     // Make a legend
1022     BuildLegend(fRatios, 0, .15,p2->GetBottomMargin()+.01,.9,
1023                 isBottom ? .6 : .4);
1024 #if 0
1025     TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,
1026                                   isBottom ? .6 : .4);
1027     l2->SetNColumns(2);
1028     l2->SetFillColor(0);
1029     l2->SetFillStyle(0);
1030     l2->SetBorderSize(0);
1031     l2->SetTextFont(132);
1032 #endif
1033     // Make a nice band from 0.9 to 1.1
1034     TGraphErrors* band = new TGraphErrors(2);
1035     band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
1036     band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
1037     band->SetPointError(0, 0, .1);
1038     band->SetPointError(1, 0, .1);
1039     band->SetFillColor(kYellow+2);
1040     band->SetFillStyle(3002);
1041     band->SetLineStyle(2);
1042     band->SetLineWidth(1);
1043     band->Draw("3 same");
1044     band->DrawClone("X L same");
1045     
1046     // Replot the ratios on top
1047     fRatios->DrawClone("nostack e1 same");
1048
1049     if (isBottom) {
1050       fRangeParam->fMasterAxis = FindXAxis(p2, fRatios->GetName());
1051       p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)", 
1052                                 fRangeParam));
1053     }
1054     else { 
1055       fRangeParam->fSlave2Axis = FindXAxis(p2, fRatios->GetName());
1056       fRangeParam->fSlave2Pad  = p2;
1057     }
1058   }
1059   //__________________________________________________________________
1060   /** 
1061    * Plot the asymmetries
1062    * 
1063    * @param max     Maximum diviation from 1 
1064    * @param y1      Lower y coordinate of pad
1065    * @param y2      Upper y coordinate of pad
1066    */
1067   void PlotLeftRight(Double_t max, Double_t y1, Double_t y2) 
1068   {
1069     if (!fShowLeftRight) return;
1070
1071     bool isBottom = (y1 < 0.0001);
1072     Double_t yd = y2 - y1;
1073     // Make a sub-pad for the result itself
1074     TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
1075     p3->SetTopMargin(0.001);
1076     p3->SetRightMargin(0.03);
1077     p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
1078     p3->SetGridx();
1079     p3->SetTicks(1,1);
1080     p3->SetNumber(2);
1081     p3->Draw();
1082     p3->cd();
1083
1084     // Fix up axis
1085     FixAxis(fLeftRight, yd, "Right/Left", 4);
1086
1087     fLeftRight->SetMaximum(1+TMath::Max(.12,1.05*max));
1088     fLeftRight->SetMinimum(1-TMath::Max(.15,1.05*max));
1089     p3->Clear();
1090     fLeftRight->DrawClone("nostack e1");
1091
1092     
1093     // Make a legend
1094     Double_t xx1 = (fCentAxis ? .7                           : .15); 
1095     Double_t xx2 = (fCentAxis ? 1-p3->GetRightMargin()-.01   : .90);
1096     Double_t yy1 = p3->GetBottomMargin()+.01;
1097     Double_t yy2 = (fCentAxis ? 1-p3->GetTopMargin()-.01-.15 : .5);
1098     BuildLegend(fLeftRight, 0, xx1, yy1, xx2, yy2);
1099     // TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
1100     // l2->SetNColumns(2);
1101     // l2->SetFillColor(0);
1102     // l2->SetFillStyle(0);
1103     // l2->SetBorderSize(0);
1104     // l2->SetTextFont(132);
1105
1106     // Make a nice band from 0.9 to 1.1
1107     TGraphErrors* band = new TGraphErrors(2);
1108     band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
1109     band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
1110     band->SetPointError(0, 0, .05);
1111     band->SetPointError(1, 0, .05);
1112     band->SetFillColor(kYellow+2);
1113     band->SetFillStyle(3002);
1114     band->SetLineStyle(2);
1115     band->SetLineWidth(1);
1116     band->Draw("3 same");
1117     band->DrawClone("X L same");
1118
1119     fLeftRight->DrawClone("nostack e1 same");
1120     if (isBottom) {
1121       fRangeParam->fMasterAxis = FindXAxis(p3, fLeftRight->GetName());
1122       p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)", 
1123                                 fRangeParam));
1124     }
1125     else { 
1126       fRangeParam->fSlave2Axis = FindXAxis(p3, fLeftRight->GetName());
1127       fRangeParam->fSlave2Pad  = p3;
1128     }
1129   }
1130   /** @} */
1131   //==================================================================
1132   /** 
1133    * @{ 
1134    * @name Data utility functions 
1135    */
1136   /** 
1137    * Get a result from the passed list
1138    * 
1139    * @param list List to search 
1140    * @param name Object name to search for 
1141    * 
1142    * @return 
1143    */
1144   TH1* FetchResult(const TList* list, const char* name) const 
1145   {
1146     if (!list) return 0;
1147     
1148     TH1* ret = static_cast<TH1*>(list->FindObject(name));
1149 #if 0
1150     if (!ret) {
1151       // all->ls();
1152       Warning("GetResult", "Histogram %s not found", name);
1153     }
1154 #endif
1155
1156     return ret;
1157   }
1158   //__________________________________________________________________
1159   /** 
1160    * Add a histogram to the stack after possibly rebinning it  
1161    * 
1162    * @param stack   Stack to add to 
1163    * @param hist    histogram
1164    * @param option  Draw options 
1165    * 
1166    * @return Maximum of histogram 
1167    */
1168   Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option="") const 
1169   {
1170     // Check if we have input 
1171     if (!hist) return 0;
1172
1173     // Rebin if needed 
1174     Rebin(hist);
1175
1176     stack->Add(hist, option);
1177     return hist->GetMaximum();
1178   }
1179   //__________________________________________________________________
1180   /** 
1181    * Add a histogram to the stack after possibly rebinning it  
1182    * 
1183    * @param stack   Stack to add to 
1184    * @param hist    histogram
1185    * @param option  Draw options 
1186    * @param sym     On return, the data symmetriced (added to stack)
1187    * 
1188    * @return Maximum of histogram 
1189    */
1190   Double_t AddHistogram(THStack* stack, TH1* hist, TH1*& sym,
1191                         Option_t* option="") const 
1192   {
1193     // Check if we have input 
1194     if (!hist) return 0;
1195
1196     // Rebin if needed 
1197     Rebin(hist);
1198     stack->Add(hist, option);
1199
1200     // Now symmetrice the histogram 
1201     sym = Symmetrice(hist);
1202     stack->Add(sym, option);
1203
1204     return hist->GetMaximum();
1205   }
1206
1207   //__________________________________________________________________
1208   /** 
1209    * Rebin a histogram 
1210    * 
1211    * @param h     Histogram to rebin
1212    * 
1213    * @return 
1214    */
1215   virtual void Rebin(TH1* h) const
1216   { 
1217     if (fRebin <= 1) return;
1218
1219     Int_t nBins = h->GetNbinsX();
1220     if(nBins % fRebin != 0) {
1221       Warning("Rebin", "Rebin factor %d is not a devisor of current number "
1222               "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
1223       return;
1224     }
1225     
1226     // Make a copy 
1227     TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
1228     tmp->Rebin(fRebin);
1229     tmp->SetDirectory(0);
1230     tmp->Reset();
1231     // The new number of bins 
1232     Int_t nBinsNew = nBins / fRebin;
1233     for(Int_t i = 1;i<= nBinsNew; i++) {
1234       Double_t content = 0;
1235       Double_t sumw    = 0;
1236       Double_t wsum    = 0;
1237       Int_t    nbins   = 0;
1238       for(Int_t j = 1; j<=fRebin;j++) {
1239         Int_t    bin = (i-1)*fRebin + j;
1240         Double_t c   =  h->GetBinContent(bin);
1241
1242         if (c <= 0) continue;
1243
1244         if (fCutEdges) {
1245           if (h->GetBinContent(bin+1)<=0 || 
1246               h->GetBinContent(bin-1)) {
1247             Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)", 
1248                     bin, c, h->GetName(), 
1249                     bin+1, h->GetBinContent(bin+1), 
1250                     bin-1, h->GetBinContent(bin-1));
1251             continue;
1252           }     
1253         }
1254         Double_t e =  h->GetBinError(bin);
1255         Double_t w =  1 / (e*e); // 1/c/c
1256         content    += c;
1257         sumw       += w;
1258         wsum       += w * c;
1259         nbins++;
1260       }
1261       
1262       if(content > 0 && nbins > 1 ) {
1263         tmp->SetBinContent(i, wsum / sumw);
1264         tmp->SetBinError(i,1./TMath::Sqrt(sumw));
1265       }
1266     }
1267
1268     // Finally, rebin the histogram, and set new content
1269     h->Rebin(fRebin);
1270     h->Reset();
1271     for(Int_t i = 1; i<= nBinsNew; i++) {
1272       h->SetBinContent(i,tmp->GetBinContent(i));
1273       h->SetBinError(i,  tmp->GetBinError(i));
1274     }
1275     
1276     delete tmp;
1277   }
1278   //__________________________________________________________________
1279   /** 
1280    * Make an extension of @a h to make it symmetric about 0 
1281    * 
1282    * @param h Histogram to symmertrice 
1283    * 
1284    * @return Symmetric extension of @a h 
1285    */
1286   TH1* Symmetrice(const TH1* h) const
1287   {
1288     Int_t nBins = h->GetNbinsX();
1289     TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
1290     s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
1291     s->Reset();
1292     s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
1293     s->SetMarkerColor(h->GetMarkerColor());
1294     s->SetMarkerSize(h->GetMarkerSize());
1295     s->SetMarkerStyle(h->GetMarkerStyle()+4);
1296     s->SetFillColor(h->GetFillColor());
1297     s->SetFillStyle(h->GetFillStyle());
1298     s->SetDirectory(0);
1299
1300     // Find the first and last bin with data 
1301     Int_t first = nBins+1;
1302     Int_t last  = 0;
1303     for (Int_t i = 1; i <= nBins; i++) { 
1304       if (h->GetBinContent(i) <= 0) continue; 
1305       first = TMath::Min(first, i);
1306       last  = TMath::Max(last,  i);
1307     }
1308     
1309     Double_t xfirst = h->GetBinCenter(first-1);
1310     Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
1311     Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
1312     for (Int_t i = f1, j=l2; i <= last; i++,j--) { 
1313       s->SetBinContent(j, h->GetBinContent(i));
1314       s->SetBinError(j, h->GetBinError(i));
1315     }
1316     // Fill in overlap bin 
1317     s->SetBinContent(l2+1, h->GetBinContent(first));
1318     s->SetBinError(l2+1, h->GetBinError(first));
1319     return s;
1320   }
1321   //__________________________________________________________________
1322   /** 
1323    * Calculate the left-right asymmetry of input histogram 
1324    * 
1325    * @param h   Input histogram
1326    * @param max On return, the maximum distance from 1 of the histogram
1327    * 
1328    * @return Asymmetry 
1329    */
1330   TH1* Asymmetry(TH1* h, Double_t& max)
1331   {
1332     if (!h) return 0;
1333
1334     TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
1335     // Int_t    oBins = h->GetNbinsX();
1336     // Double_t high  = h->GetXaxis()->GetXmax();
1337     // Double_t low   = h->GetXaxis()->GetXmin();
1338     // Double_t dBin  = (high - low) / oBins;
1339     // Int_t    tBins = Int_t(2*high/dBin+.5);
1340     // ret->SetBins(tBins, -high, high);
1341     ret->SetDirectory(0);
1342     ret->Reset();
1343     ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
1344     ret->SetYTitle("Right/Left");
1345     Int_t nBins = h->GetNbinsX();
1346     for (Int_t i = 1; i <= nBins; i++) { 
1347       Double_t x = h->GetBinCenter(i);
1348       if (x > 0) break;
1349       
1350       Double_t c1 = h->GetBinContent(i);
1351       Double_t e1 = h->GetBinError(i);
1352       if (c1 <= 0) continue; 
1353       
1354       Int_t    j  = h->FindBin(-x);
1355       if (j <= 0 || j > nBins) continue;
1356
1357       Double_t c2 = h->GetBinContent(j);
1358       Double_t e2 = h->GetBinError(j);
1359
1360       Double_t c12 = c1*c1;
1361       Double_t e   = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
1362       
1363       Int_t    k   = ret->FindBin(x);
1364       ret->SetBinContent(k, c2/c1);
1365       ret->SetBinError(k, e);
1366     }
1367     max = TMath::Max(max, RatioMax(ret));
1368
1369     return ret;
1370   }
1371   //__________________________________________________________________
1372   /** 
1373    * Transform a graph into a histogram 
1374    * 
1375    * @param g 
1376    * 
1377    * @return 
1378    */
1379   TH1* Graph2Hist(const TGraphAsymmErrors* g) const
1380   {
1381     Int_t    nBins = g->GetN();
1382     TArrayF  bins(nBins+1);
1383     Double_t dx = 0;
1384     for (Int_t i = 0; i < nBins; i++) { 
1385       Double_t x   = g->GetX()[i];
1386       Double_t exl = g->GetEXlow()[i];
1387       Double_t exh = g->GetEXhigh()[i];
1388       bins.fArray[i]   = x-exl;
1389       bins.fArray[i+1] = x+exh;
1390       Double_t dxi = exh+exl;
1391       if (i == 0) dx  = dxi;
1392       else if (dxi != dx) dx = 0;
1393     }
1394     TString name(g->GetName());
1395     TString title(g->GetTitle());
1396     TH1D* h = 0;
1397     if (dx != 0) {
1398       h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
1399     }
1400     else {
1401       h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
1402     }
1403     h->SetMarkerStyle(g->GetMarkerStyle());
1404     h->SetMarkerColor(g->GetMarkerColor());
1405     h->SetMarkerSize(g->GetMarkerSize());
1406     
1407     return h;
1408   }
1409   /* @} */
1410   //==================================================================
1411   /** 
1412    * @{ 
1413    * @name Ratio utility functions 
1414    */
1415   /** 
1416    * Get the maximum diviation from 1 in the passed ratio
1417    * 
1418    * @param h Ratio histogram
1419    * 
1420    * @return Max diviation from 1 
1421    */
1422   Double_t RatioMax(TH1* h) const
1423   {
1424     Int_t    nBins = h->GetNbinsX();
1425     Double_t ret   = 0;
1426     for (Int_t i = 1; i <= nBins; i++) { 
1427       Double_t c = h->GetBinContent(i);
1428       if (c == 0) continue;
1429       Double_t e = h->GetBinError(i);
1430       Double_t d = TMath::Abs(1-c-e);
1431       ret        = TMath::Max(d, ret);
1432     }
1433     return ret;
1434   }
1435   //__________________________________________________________________
1436   /** 
1437    * Make the ratio of h1 to h2 
1438    * 
1439    * @param o1  First object (numerator) 
1440    * @param o2  Second object (denominator)
1441    * @param max Maximum diviation from 1 
1442    * 
1443    * @return o1 / o2
1444    */
1445   TH1* Ratio(const TObject* o1, const TObject* o2, Double_t& max) const
1446   {
1447     if (!o1 || !o2) return 0;
1448
1449     TH1*        r  = 0;
1450     const TAttMarker* m1 = 0;
1451     const TAttMarker* m2 = 0;
1452     const TH1* h1 = dynamic_cast<const TH1*>(o1); 
1453     if (h1) { 
1454       m1 = h1;
1455       const TH1* h2 = dynamic_cast<const TH1*>(o2); 
1456       if (h2) { 
1457         m2 = h2;
1458         r = RatioHH(h1,h2);
1459       }
1460       else {
1461         const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
1462         if (g2) { 
1463           m2 = g2;
1464           r = RatioHG(h1,g2);      
1465         }
1466       }
1467     }
1468     else {
1469       const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
1470       if (g1) { 
1471         m1 = g1;
1472         const TGraphAsymmErrors* g2 = 
1473           dynamic_cast<const TGraphAsymmErrors*>(o2);
1474         if (g2) {
1475           m2 = g2;
1476           r = RatioGG(g1, g2);
1477         }
1478       }
1479     }
1480     if (!r) {
1481       Warning("Ratio", "Don't know how to divide a %s (%s) with a %s (%s)", 
1482               o1->ClassName(), o1->GetName(), o2->ClassName(), o2->GetName());
1483       return 0;
1484     }
1485     // Check that the histogram isn't empty
1486     if (r->GetEntries() <= 0) { 
1487       delete r; 
1488       r = 0; 
1489     }
1490     if (r) {
1491       r->SetMarkerStyle(m2->GetMarkerStyle());
1492       r->SetMarkerColor(m1->GetMarkerColor());
1493       r->SetMarkerSize(0.9*m1->GetMarkerSize());
1494       r->SetName(Form("%s_over_%s", o1->GetName(), o2->GetName()));
1495       r->SetTitle(Form("%s / %s", o1->GetTitle(), o2->GetTitle()));
1496       r->SetDirectory(0);
1497       max = TMath::Max(RatioMax(r), max);
1498     }
1499
1500     return r;
1501   }
1502   //__________________________________________________________________
1503   /** 
1504    * Compute the ratio of @a h to @a g.  @a g is evaluated at the bin
1505    * centers of @a h 
1506    * 
1507    * @param h  Numerator 
1508    * @param g  Divisor 
1509    * 
1510    * @return h/g 
1511    */
1512   TH1* RatioHG(const TH1* h, const TGraph* g) const 
1513   {
1514     if (!h || !g) return 0;
1515
1516     TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
1517     ret->Reset();
1518     Double_t xlow  = g->GetX()[0];
1519     Double_t xhigh = g->GetX()[g->GetN()-1];
1520     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
1521
1522     for (Int_t i = 1; i <= h->GetNbinsX(); i++) { 
1523       Double_t c = h->GetBinContent(i);
1524       if (c <= 0) continue;
1525
1526       Double_t x = h->GetBinCenter(i);
1527       if (x < xlow || x > xhigh) continue; 
1528
1529       Double_t f = g->Eval(x);
1530       if (f <= 0) continue; 
1531
1532       ret->SetBinContent(i, c / f);
1533       ret->SetBinError(i, h->GetBinError(i) / f);
1534     }
1535     return ret;
1536   }
1537   //__________________________________________________________________
1538   /** 
1539    * Make the ratio of h1 to h2 
1540    * 
1541    * @param h1 First histogram (numerator) 
1542    * @param h2 Second histogram (denominator)
1543    * 
1544    * @return h1 / h2
1545    */
1546   TH1* RatioHH(const TH1* h1, const TH1* h2) const
1547   {
1548     if (!h1 || !h2) return 0;
1549     TH1* t1 = static_cast<TH1*>(h1->Clone("tmp"));
1550     t1->Divide(h2);
1551     return t1;
1552   }
1553   //__________________________________________________________________
1554   /** 
1555    * Calculate the ratio of two graphs - g1 / g2
1556    * 
1557    * @param g1 Numerator 
1558    * @param g2 Denominator
1559    * 
1560    * @return g1 / g2 in a histogram 
1561    */
1562   TH1* RatioGG(const TGraphAsymmErrors* g1, 
1563                const TGraphAsymmErrors* g2) const
1564   {
1565     Int_t    nBins = g1->GetN();
1566     TArrayF  bins(nBins+1);
1567     Double_t dx   = 0;
1568     for (Int_t i = 0; i < nBins; i++) {
1569       Double_t x   = g1->GetX()[i];
1570       Double_t exl = g1->GetEXlow()[i];
1571       Double_t exh = g1->GetEXhigh()[i];
1572       bins.fArray[i]   = x-exl;
1573       bins.fArray[i+1] = x+exh;
1574       Double_t dxi = exh+exl;
1575       if (i == 0) dx  = dxi;
1576       else if (dxi != dx) dx = 0;
1577     }
1578     TH1* h = 0;
1579     if (dx != 0) {
1580       h = new TH1F("tmp", "tmp", nBins, bins[0], bins[nBins]);
1581     }
1582     else {
1583       h = new TH1F("tmp", "tmp", nBins, bins.fArray);
1584     }
1585
1586     Double_t low  = g2->GetX()[0];
1587     Double_t high = g2->GetX()[g2->GetN()-1];
1588     if (low > high) { Double_t t = low; low = high; high = t; }
1589     for (Int_t i = 0; i < nBins; i++) { 
1590       Double_t x  = g1->GetX()[i];
1591       if (x < low-dx || x > high+dx) continue;
1592       Double_t c1 = g1->GetY()[i];
1593       Double_t e1 = g1->GetErrorY(i);
1594       Double_t c2 = g2->Eval(x);
1595       
1596       h->SetBinContent(i+1, c1 / c2);
1597       h->SetBinError(i+1, e1 / c2);
1598     }
1599     return h;
1600   }  
1601   /* @} */
1602   //==================================================================
1603   /** 
1604    * @{ 
1605    * @name Graphics utility functions 
1606    */
1607   /** 
1608    * Find an X axis in a pad 
1609    * 
1610    * @param p     Pad
1611    * @param name  Histogram to find axis for 
1612    * 
1613    * @return Found axis or null
1614    */
1615   TAxis* FindXAxis(TVirtualPad* p, const char* name)
1616   {
1617     TObject* o = p->GetListOfPrimitives()->FindObject(name);
1618     if (!o) { 
1619       Warning("FindXAxis", "%s not found in pad", name);
1620       return 0;
1621     }
1622     THStack* stack = dynamic_cast<THStack*>(o);
1623     if (!stack) { 
1624       Warning("FindXAxis", "%s is not a THStack", name);
1625       return 0;
1626     }
1627     if (!stack->GetHistogram()) { 
1628       Warning("FindXAxis", "%s has no histogram", name);
1629       return 0;
1630     }
1631     TAxis* ret = stack->GetHistogram()->GetXaxis();
1632     return ret;
1633   }
1634
1635   //__________________________________________________________________
1636   /**
1637    * Fix the apperance of the axis in a stack
1638    *
1639    * @param stack  stack of histogram
1640    * @param s      Scaling factor
1641    * @param ytitle Y axis title
1642    * @param force  Whether to draw the stack first or not
1643    * @param ynDiv  Divisions on Y axis
1644    */
1645   void FixAxis(THStack* stack, Double_t yd, const char* ytitle,
1646                Int_t ynDiv=210, Bool_t force=true)
1647   {
1648     if (!stack) { 
1649       Warning("FixAxis", "No stack passed for %s!", ytitle);
1650       return;
1651     }
1652     if (force) stack->Draw("nostack e1");
1653
1654     TH1* h = stack->GetHistogram();
1655     if (!h) { 
1656       Warning("FixAxis", "Stack %s has no histogram", stack->GetName());
1657       return;
1658     }
1659     Double_t s = 1/yd/1.2;
1660     // Info("FixAxis", "for %s, s=1/%f=%f", stack->GetName(), yd, s);
1661
1662     h->SetXTitle("#font[152]{#eta}");
1663     h->SetYTitle(ytitle);
1664     TAxis* xa = h->GetXaxis();
1665     TAxis* ya = h->GetYaxis();
1666
1667     // Int_t   npixels = 20
1668     // Float_t dy = gPad->AbsPixeltoY(0) - gPad->AbsPixeltoY(npixels);
1669     // Float_t ts = dy/(gPad->GetY2() - gPad->GetY1());
1670
1671     if (xa) {
1672       // xa->SetTitle(h->GetXTitle());
1673       // xa->SetTicks("+-");
1674       xa->SetTitleSize(s*xa->GetTitleSize());
1675       xa->SetLabelSize(s*xa->GetLabelSize());
1676       xa->SetTickLength(s*xa->GetTickLength());
1677       // xa->SetTitleOffset(xa->GetTitleOffset()/s);
1678
1679       if (stack != fResults) {
1680         TAxis* rxa = fResults->GetXaxis();
1681         xa->Set(rxa->GetNbins(), rxa->GetXmin(), rxa->GetXmax());
1682       }
1683     }
1684     if (ya) {
1685       // ya->SetTitle(h->GetYTitle());
1686       ya->SetDecimals();
1687       // ya->SetTicks("+-");
1688       ya->SetNdivisions(ynDiv);
1689       ya->SetTitleSize(s*ya->GetTitleSize());
1690       ya->SetTitleOffset(ya->GetTitleOffset()/s);
1691       ya->SetLabelSize(s*ya->GetLabelSize());
1692     }
1693   }
1694   //__________________________________________________________________
1695   /** 
1696    * Merge two histograms into one 
1697    * 
1698    * @param cen    Central part
1699    * @param fwd    Forward part
1700    * @param xlow   On return, lower eta bound
1701    * @param xhigh  On return, upper eta bound
1702    * 
1703    * @return Newly allocated histogram or null
1704    */
1705   TH1* 
1706   Merge(const TH1* cen, const TH1* fwd, Double_t& xlow, Double_t& xhigh)
1707   {
1708     TH1* tmp = static_cast<TH1*>(fwd->Clone("tmp"));
1709     TString name(fwd->GetName());
1710     name.ReplaceAll("Forward", "Merged");
1711     tmp->SetName(name);
1712
1713     // tmp->SetMarkerStyle(28);
1714     // tmp->SetMarkerColor(kBlack);
1715     tmp->SetDirectory(0);
1716     xlow  = 100;
1717     xhigh = -100;
1718     for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
1719       Double_t cc = cen->GetBinContent(i);
1720       Double_t cf = fwd->GetBinContent(i);
1721       Double_t ec = cen->GetBinError(i);
1722       Double_t ef = fwd->GetBinError(i);
1723       Double_t nc = cf;
1724       Double_t ne = ef;
1725       if (cc < 0.001 && cf < 0.01) continue;
1726       xlow  = TMath::Min(tmp->GetXaxis()->GetBinLowEdge(i),xlow);
1727       xhigh = TMath::Max(tmp->GetXaxis()->GetBinUpEdge(i),xhigh);
1728       if (cc > 0.001) {
1729         nc = cc;
1730         ne = ec;
1731         if (cf > 0.001) {
1732           nc  = (cf + cc) / 2;
1733           ne  = TMath::Sqrt(ec*ec + ef*ef);
1734         }
1735       }
1736       tmp->SetBinContent(i, nc);
1737       tmp->SetBinError(i, ne);
1738     }
1739     return tmp;
1740   }
1741   //____________________________________________________________________
1742   /** 
1743    * Fit  @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$ to histogram data 
1744    * 
1745    * @param tmp    Histogram
1746    * @param xlow   Lower x bound
1747    * @param xhigh  Upper x bound 
1748    *
1749    * @return Fitted function 
1750    */
1751   TF1* 
1752   FitMerged(TH1* tmp, Double_t xlow, Double_t xhigh)
1753   {
1754     TF1* tmpf  = new TF1("tmpf", "gaus", xlow, xhigh);
1755     tmp->Fit(tmpf, "NQ", "");
1756     tmp->SetDirectory(0);
1757
1758     TF1* fit = new TF1("f", myFunc, xlow, xhigh, 4);
1759     fit->SetParNames("a_{1}", "a_{2}", "#sigma_{1}", "#sigma_{2}");
1760     fit->SetParameters(tmpf->GetParameter(0), 
1761                        .2, 
1762                        tmpf->GetParameter(2), 
1763                        tmpf->GetParameter(2)/4);
1764     fit->SetParLimits(3, 0, 100);
1765     fit->SetParLimits(4, 0, 100);
1766     tmp->Fit(fit,"0W","");
1767
1768     delete tmpf;
1769     return fit;
1770   }
1771   //____________________________________________________________________
1772   /** 
1773    * Make band of systematic errors 
1774    * 
1775    * @param tmp Histogram
1776    * @param fit Fit 
1777    */
1778   void
1779   MakeSysError(TH1* tmp, TH1* cen, TH1* fwd, TF1* fit)
1780   {
1781     for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
1782       Double_t tc = tmp->GetBinContent(i);
1783       if (tc < 0.01) continue;
1784       Double_t fc = fwd->GetBinContent(i);
1785       Double_t cc = cen->GetBinContent(i);
1786       Double_t sysErr = fFwdSysErr;
1787       if (cc > .01 && fc > 0.01) 
1788         sysErr = (fFwdSysErr+fCenSysErr) / 2;
1789       else if (cc > .01) 
1790         sysErr = fCenSysErr;
1791       Double_t x = tmp->GetXaxis()->GetBinCenter(i);
1792       Double_t y = fit->Eval(x);
1793       tmp->SetBinContent(i, y);
1794       tmp->SetBinError(i,sysErr*y);
1795     }
1796     TString name(tmp->GetName());
1797     name.ReplaceAll("Merged", "SysError");
1798     tmp->SetName(name);
1799     tmp->SetMarkerColor(SYSERR_COLOR);
1800     tmp->SetLineColor(SYSERR_COLOR);
1801     tmp->SetFillColor(SYSERR_COLOR);
1802     tmp->SetFillStyle(SYSERR_STYLE);
1803     tmp->SetMarkerStyle(0);
1804     tmp->SetLineWidth(0);
1805   }
1806   void CorrectForward(TH1* h) const
1807   {
1808     if (!fRemoveOuters) return;
1809     
1810     for (Int_t i = 1; i <= h->GetNbinsX(); i++) { 
1811       Double_t eta = h->GetBinCenter(i);
1812       if (TMath::Abs(eta) < 2.3) { 
1813         h->SetBinContent(i, 0);
1814         h->SetBinError(i, 0);
1815       }
1816     }
1817   }
1818   void CorrectCentral(TH1* h) const 
1819   {
1820     if (fClusterScale.IsNull()) return;
1821     TString t(h->GetTitle());
1822     Info("CorrectCentral", "Replacing Central with Tracklets in %s", t.Data());
1823     t.ReplaceAll("Central", "Tracklets");
1824     h->SetTitle(t);
1825
1826     TF1* cf = new TF1("clusterScale", fClusterScale);
1827     for (Int_t i = 1; i <= h->GetNbinsX(); i++) { 
1828       Double_t eta = h->GetBinCenter(i);
1829       Double_t f   = cf->Eval(eta);
1830       Double_t c   = h->GetBinContent(i);
1831       if (f < .1) f = 1;
1832       h->SetBinContent(i, c / f);
1833     }
1834     delete cf;
1835   }
1836   //____________________________________________________________________
1837   void Export(const char* basename)
1838   {
1839     TString bname(basename);
1840     bname.ReplaceAll(" ", "_");
1841     bname.ReplaceAll("-", "_");
1842     TString fname(Form("%s.C", bname.Data()));
1843
1844     std::ofstream out(fname.Data());
1845     if (!out) { 
1846       Error("Export", "Failed to open output file %s", fname.Data());
1847       return;
1848     }
1849     out << "// Create by dNdetaDrawer\n"
1850         << "void " << bname << "(THStack* stack, TLegend* l, Int_t m)\n"
1851         << "{"
1852         << "   Int_t ma[] = { 24, 25, 26, 32,\n"
1853         << "                  20, 21, 22, 33,\n"
1854         << "                  34, 30, 29, 0, \n"
1855         << "                  23, 27 };\n"
1856         << "   Int_t mm = ((m < 20 || m > 34) ? 0 : ma[m-20]);\n\n";
1857     TList* hists = fResults->GetHists();
1858     TIter  next(hists);
1859     TH1*   hist = 0;
1860     while ((hist = static_cast<TH1*>(next()))) { 
1861       TString hname = hist->GetName();
1862       hname.Append(Form("_%04x", (gRandom->Integer(0xffff) & 0xffff)));
1863       hist->SetName(hname);
1864       hist->GetListOfFunctions()->Clear();
1865       hist->SavePrimitive(out, "nodraw");
1866       bool mirror = hname.Contains("mirror");
1867       bool syserr = hname.Contains("SysError");
1868       if (!syserr) 
1869         out << "   " << hname << "->SetMarkerStyle(" 
1870             << (mirror ? "mm" : "m") << ");\n";
1871       else 
1872         out << "   " << hname << "->SetMarkerStyle(1);\n";
1873       out << "   stack->Add(" << hname 
1874           << (syserr ? ",\"e5\"" : "") << ");\n\n";
1875     }
1876     UShort_t    snn = fSNNString->GetUniqueID();
1877     // const char* sys = fSysString->GetTitle();
1878     TString eS;
1879     if      (snn == 2750)     snn = 2760;
1880     if      (snn < 1000)      eS = Form("%3dGeV", snn);
1881     else if (snn % 1000 == 0) eS = Form("%dTeV", snn/1000);
1882     else                      eS = Form("%4.2fTeV", float(snn)/1000);
1883     out << "  if (l) {\n"
1884         << "    TLegendEntry* e = l->AddEntry(\"\",\"" << eS << "\",\"pl\");\n"
1885         << "    e->SetMarkerStyle(m);\n"
1886         << "    e->SetMarkerColor(kBlack);\n"
1887         << "  }\n"
1888         << "}\n" << std::endl;
1889   }
1890               
1891   /* @} */
1892
1893
1894
1895   //__________________________________________________________________
1896   /** 
1897    * @{ 
1898    * @name Options 
1899    */
1900   Bool_t       fShowRatios;   // Show ratios 
1901   Bool_t       fShowLeftRight;// Show asymmetry 
1902   Bool_t       fShowRings;    // Show rings too
1903   Bool_t       fExport;       // Export results to file
1904   Bool_t       fCutEdges;     // Whether to cut edges
1905   Bool_t       fRemoveOuters; // Whether to remove outers
1906   UShort_t     fShowOthers;   // Show other data
1907   /* @} */
1908   /** 
1909    * @{ 
1910    * @name Settings 
1911    */
1912   UShort_t     fRebin;        // Rebinning factor 
1913   Double_t     fFwdSysErr;    // Systematic error in forward range
1914   Double_t     fCenSysErr;    // Systematic error in central range 
1915   TString      fTitle;        // Title on plot
1916   TString      fClusterScale; // Scaling of clusters to tracklets      
1917   /* @} */
1918   /** 
1919    * @{ 
1920    * @name Read (or set) information 
1921    */
1922   TNamed*      fTrigString;   // Trigger string (read, or set)
1923   TNamed*      fNormString;   // Normalisation string (read, or set)
1924   TNamed*      fSNNString;    // Energy string (read, or set)
1925   TNamed*      fSysString;    // Collision system string (read or set)
1926   TAxis*       fVtxAxis;      // Vertex cuts (read or set)
1927   TAxis*       fCentAxis;     // Centrality axis
1928   /* @} */
1929   /** 
1930    * @{ 
1931    * @name Resulting plots 
1932    */
1933   THStack*     fResults;      // Stack of results 
1934   THStack*     fRatios;       // Stack of ratios 
1935   THStack*     fLeftRight;    // Left-right asymmetry
1936   TMultiGraph* fOthers;       // Older data 
1937   TH1*         fTriggers;     // Number of triggers
1938   TH1*         fTruth;        // Pointer to truth 
1939   /* @} */
1940   RangeParam*  fRangeParam;   // Parameter object for range zoom 
1941 };
1942
1943 //____________________________________________________________________
1944 /** 
1945  * Function to calculate 
1946  * @f[
1947  *  g(x;A_1,A_2,\sigma_1,\sigma_2) = 
1948  *       A_1\left(\frac{1}{2\pi\sigma_1}e^{(x/\sigma_1)^2} - 
1949  *           A_2\frac{1}{2\pi\sigma_2}e^{(x/\sigma_2)^2}\right)
1950  * @f]
1951  * 
1952  * @param xp Pointer to x array
1953  * @param pp Pointer to parameter array 
1954  * 
1955  * @return @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$
1956  */
1957 Double_t myFunc(Double_t* xp, Double_t* pp)
1958 {
1959   Double_t x  = xp[0];
1960   Double_t a1 = pp[0];
1961   Double_t a2 = pp[1];
1962   Double_t s1 = pp[2];
1963   Double_t s2 = pp[3];
1964   return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
1965 }
1966
1967 //=== Stuff for auto zooming =========================================
1968 void UpdateRange(dNdetaDrawer::RangeParam* p)
1969 {
1970   if (!p) { 
1971     Warning("UpdateRange", "No parameters %p", p);
1972     return;
1973   }
1974   if (!p->fMasterAxis) { 
1975     Warning("UpdateRange", "No master axis %p", p->fMasterAxis);
1976     return;
1977   }
1978   Int_t    first = p->fMasterAxis->GetFirst();
1979   Int_t    last  = p->fMasterAxis->GetLast();
1980   Double_t x1    = p->fMasterAxis->GetBinCenter(first);
1981   Double_t x2    = p->fMasterAxis->GetBinCenter(last);
1982
1983   if (p->fSlave1Axis) { 
1984     Int_t i1 = p->fSlave1Axis->FindBin(x1);
1985     Int_t i2 = p->fSlave1Axis->FindBin(x2);
1986     p->fSlave1Axis->SetRange(i1, i2);
1987     p->fSlave1Pad->Modified();
1988     p->fSlave1Pad->Update();
1989   }
1990   if (p->fSlave2Axis) { 
1991     Int_t i1 = p->fSlave2Axis->FindBin(x1);
1992     Int_t i2 = p->fSlave2Axis->FindBin(x2);
1993     p->fSlave2Axis->SetRange(i1, i2);
1994     p->fSlave2Pad->Modified();
1995     p->fSlave2Pad->Update();
1996   }
1997   TCanvas*  c = gPad->GetCanvas();
1998   c->cd();
1999 }
2000   
2001 //____________________________________________________________________
2002 void RangeExec(dNdetaDrawer::RangeParam* p)
2003 {
2004   // Event types: 
2005   //  51:   Mouse move 
2006   //  53:   
2007   //  1:    Button down 
2008   //  21:   Mouse drag
2009   //  11:   Mouse release 
2010   // dNdetaDrawer::RangeParam* p = 
2011   //   reinterpret_cast<dNdetaDrawer::RangeParam*>(addr);
2012   Int_t event     = gPad->GetEvent();
2013   TObject *select = gPad->GetSelected();
2014   if (event == 53) { 
2015     UpdateRange(p);
2016     return;
2017   }
2018   if (event != 11 || !select || select != p->fMasterAxis) return;
2019   UpdateRange(p);
2020 }
2021
2022 //=== Steering functions
2023 //==============================================  
2024 void
2025 Usage()
2026 {
2027   Info("DrawdNdeta", "Usage: DrawdNdeta(FILE,TITLE,REBIN,OTHERS,FLAGS)\n\n"
2028        "  const char* FILE   File name to open (\"forward_root\")\n"
2029        "  const char* TITLE  Title to put on plot (\"\")\n"
2030        "  UShort_t    REBIN  Rebinning factor (1)\n"
2031        "  UShort_t    OTHERS Other data to draw - more below (0x7)\n"
2032        "  UShort_t    FLAGS  Visualisation flags - more below (0x7)\n\n"
2033        " OTHERS is a bit mask of\n\n"
2034        "  0x1   Show UA5 data (INEL,NSD, ppbar, 900GeV)\n"
2035        "  0x2   Show CMS data (NSD, pp)\n"
2036        "  0x4   Show published ALICE data (INEL,INEL>0,NSD, pp)\n"
2037        "  0x8   Show event genertor data\n\n"
2038        " FLAGS is a bit mask of\n\n"
2039        "  0x1   Show ratios of data to other data and possibly MC\n"
2040        "  0x2   Show left-right asymmetry\n"
2041        "  0x4   Show systematic error band\n"
2042        "  0x8   Show individual ring results (INEL only)\n"
2043        "  0x10  Cut edges when rebinning\n"
2044        "  0x20  Remove FMDxO points\n"
2045        "  0x40  Do not make our own canvas\n"
2046        );
2047 }
2048
2049 //____________________________________________________________________
2050 /** 
2051  * Draw @f$ dN/d\eta@f$ 
2052  * 
2053  * @param filename  File name 
2054  * @param flags     Flags 
2055  * @param title     Title 
2056  * @param rebin     Rebinning factor 
2057  * @param cutEdges  Whether to cut edges when rebinning
2058  * @param sNN       (optional) Collision energy [GeV]
2059  * @param sys       (optional) Collision system (1: pp, 2: PbPb)
2060  * @param trg       (optional) Trigger (1: INEL, 2: INEL>0, 4: NSD)   
2061  * @param vzMin     Least @f$ v_z@f$
2062  * @param vzMax     Largest @f$ v_z@f$
2063  *
2064  * @ingroup pwg2_forward_dndeta
2065  */
2066 void
2067 DrawdNdeta(const char* filename="forward_dndeta.root", 
2068            const char* title="",
2069            UShort_t    rebin=5, 
2070            UShort_t    others=0x7,
2071            UShort_t    flags=0x7,
2072            UShort_t    sNN=0, 
2073            UShort_t    sys=0,
2074            UShort_t    trg=0,
2075            Float_t     vzMin=999, 
2076            Float_t     vzMax=-999)
2077 {
2078   TString fname(filename);
2079   fname.ToLower();
2080   if (fname.CompareTo("help") == 0 || 
2081       fname.CompareTo("--help") == 0) { 
2082     Usage();
2083     return;
2084   }
2085   dNdetaDrawer* pd = new dNdetaDrawer;
2086   dNdetaDrawer& d = *pd;
2087   d.SetRebin(rebin);
2088   d.SetTitle(title);
2089   d.SetShowOthers(others);
2090   d.SetShowRatios(flags & 0x1);
2091   d.SetShowLeftRight(flags & 0x2);
2092   d.SetForwardSysError(flags & 0x4 ? 0.076 : 0);
2093   d.SetShowRings(flags & 0x8);
2094   d.SetCutEdges(flags & 0x10);
2095   d.fRemoveOuters = (flags & 0x20);
2096   d.SetExport(flags & 0x40);
2097   // d.fClusterScale = "1.06 -0.003*x +0.0119*x*x";
2098   // Do the below if your input data does not contain these settings 
2099   if (sNN > 0) d.SetSNN(sNN);     // Collision energy per nucleon pair (GeV)
2100   if (sys > 0) d.SetSys(sys);     // Collision system (1:pp, 2:PbPB)
2101   if (trg > 0) d.SetTrigger(trg); // Collision trigger (1:INEL, 2:INEL>0, 4:NSD)
2102   if (vzMin < 999 && vzMax > -999) 
2103     d.SetVertexRange(vzMin,vzMax); // Collision vertex range (cm)
2104   d.Run(filename);
2105 }
2106 //____________________________________________________________________
2107 //
2108 // EOF
2109 //
2110