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