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