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