PWG2/FORWARD/analysis2/older/DrawResKeep.C

   1 #include <TH1D.h>
   2 #include <TH2D.h>
   3 #include <TH1I.h>
   4 #include <TProfile.h>
   5 #include <TList.h>
   6 #include <TAxis.h>
   7 #include <TCanvas.h>
   8 #include <TPad.h>
   9 #include <TFile.h>
  10 #include <TTree.h>
  11 #include <TError.h>
  12 #include <TStyle.h>
  13 #include <THStack.h>
  14 #include <TLegend.h>
  15 #include <TMath.h>
  16 #include "AliAODForwardMult.h"
  17
  18 /**
  19  * @defgroup pwg2_forward_analysis_scripts PWG2 Forward analysis - scripts
  20  *
  21  * @ingroup pwg2_forward_analysis
  22  */
  23 /**
  24  * Example macro to loop over the event-by-event 2D histogram of
  25  * @f[
  26  *   \frac{d^{2}N_{ch}}{d\eta\,d\phi}
  27  * @f]
  28  * stored in an AOD.
  29  *
  30  * The class needs the files &lt;<i>base</i>&gt;<tt>_hists.root</tt>
  31  * containing the histograms generated by AliForwardMultiplicity and
  32  * the file &lt;<i>base</i>&gt;<tt>_aods.root</tt> containing the tree
  33  * with AliAODEvent objects where AliAODForwardMult objects have been
  34  * added to in the branch <tt>Forward</tt>
  35  *
  36  * @ingroup pwg2_forward_analysis_scripts
  37  */
  38 class DrawRes
  39 {
  40 public:
  41   /**
  42    * Constructor
  43    *
  44    * @param special If true, add to the list of 'specials'
  45    */
  46   DrawRes(Bool_t special=true)
  47     : fVzMin(-10),
  48       fVzMax(10),
  49       fBinVzMin(0),
  50       fBinVzMax(0),
  51       fRebin(1),
  52       fTree(0),
  53       fAOD(0),
  54       fEvents(0),
  55       fCollect(0),
  56       fByVtx1D(0),
  57       fTotal1D(0),
  58       fTotal2D(0),
  59       fFromVtx(0),
  60       fNorm(0)
  61   {}
  62   //__________________________________________________________________
  63   /**
  64    * Open the files &lt;<i>base</i>&gt;<tt>_hists.root</tt>
  65    * containing the histograms generated by AliForwardMultiplicity and
  66    * the file &lt;<i>base</i>&gt;<tt>_aods.root</tt> containing the tree
  67    * with AliAODEvent objects.
  68    *
  69    * @param base  Base name of files
  70    * @param vzMin Minimum collision vertex z position to use
  71    * @param vzMax Maximum collision vertex z position to use
  72    * @param rebin Rebinning factor
  73    *
  74    * @return true on success, false otherwise
  75    */
  76   Bool_t Open(const char* base,
  77               Double_t    vzMin=-10,
  78               Double_t    vzMax=10,
  79               Int_t       rebin=1)
  80   {
  81     // Set our cuts etc.
  82     fVzMin = vzMin;
  83     fVzMax = vzMax;
  84     if (fVzMax < fVzMin && fVzMin < 0) fVzMax = -fVzMin;
  85     fRebin = rebin;
  86
  87     // Clear previously created data objects
  88     fByVtx1D.Delete();
  89     if (fTotal1D) { delete fTotal1D; fTotal1D = 0; }
  90     if (fTotal2D) { delete fTotal2D; fTotal2D = 0; }
  91     if (fTree && fTree->GetCurrentFile()) {
  92       fTree->GetCurrentFile()->Close();
  93     }
  94     fCollect = 0;
  95     fTree    = 0;
  96
  97     // Open the AOD file
  98     TString fn   = TString::Format("%s_aods.root", base);
  99     TFile*  file = TFile::Open(fn.Data(), "READ");
 100     if (!file) {
 101       Error("Init", "Couldn't open %s", fn.Data());
 102       return kFALSE;
 103     }
 104
 105     // Get the AOD tree
 106     fTree = static_cast<TTree*>(file->Get("aodTree"));
 107     if (!fTree) {
 108       Error("Init", "Couldn't get the tree");
 109       return kFALSE;
 110     }
 111
 112     // Set the branch pointer
 113     fTree->SetBranchAddress("Forward", &fAOD);
 114
 115     // Open the histogram file
 116     fn   = TString::Format("%s_hists.root", base);
 117     file = TFile::Open(fn.Data(), "READ");
 118     if (!file) {
 119       Error("Init", "Couldn't open %s", fn.Data());
 120       return kFALSE;
 121     }
 122
 123     // Get the list stored in the file
 124     TList* forward =
 125       static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
 126     if (!forward)  {
 127       Error("Init", "Couldn't get forward list");
 128       return kFALSE;
 129     }
 130
 131     // Get the list from the collector
 132     fCollect =
 133       static_cast<TList*>(forward->FindObject("fmdHistCollector"));
 134     if (!fCollect)  {
 135       Error("Init", "Couldn't get collector list");
 136       return kFALSE;
 137     }
 138
 139     // Get the event (by vertex bin) histogram
 140     fEvents = static_cast<TH1I*>(forward->FindObject("nEventsTrVtx"));
 141     if (!fEvents) {
 142       Error("Init", "Couldn't get the event histogram");
 143       return kFALSE;
 144     }
 145
 146     // Find the min/max bins to use based on the cuts given
 147     fBinVzMin = fEvents->FindBin(fVzMin);
 148     fBinVzMax = fEvents->FindBin(fVzMax-.0000001);
 149     Info("Open", "Selected vertex bins are [%d,%d]", fBinVzMin, fBinVzMax);
 150
 151     return kTRUE;
 152   }
 153   //__________________________________________________________________
 154   /**
 155    * Check if the passed vertex bin number [1,nVtxBins] is within our
 156    * cut.
 157    *
 158    * @param bin Vertex bin [1,nVtxBins]
 159    *
 160    * @return true if within cut, false otherwise
 161    */
 162   Bool_t IsInsideVtxCut(Int_t bin) const
 163   {
 164     return bin >= fBinVzMin && bin <= fBinVzMax;
 165   }
 166
 167   //__________________________________________________________________
 168   /**
 169    * Make a 1D histogram of @f$ dN_{ch}/d\eta@f$ with the given name and
 170    * title
 171    *
 172    * @param name   Name of histogram
 173    * @param title  Title of histogram
 174    * @param a1     Eta axis
 175    *
 176    * @return Newly allocated 1D histogram object
 177    */
 178   TH1D* Make1D(const char* name, const char* title, const TAxis* a1)
 179   {
 180     TH1D* ret = new TH1D(name, title,
 181                          a1->GetNbins(), a1->GetXmin(), a1->GetXmax());
 182     ret->SetXTitle("#eta");
 183     ret->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
 184     ret->Sumw2();
 185     ret->SetMarkerColor(kRed+1);
 186     ret->SetLineColor(kRed+1);
 187     ret->SetMarkerStyle(24);
 188     ret->SetMarkerSize(1);
 189     ret->SetStats(0);
 190     ret->SetDirectory(0);
 191
 192     return ret;
 193   }
 194   //__________________________________________________________________
 195   /**
 196    * Make a 2D histogram of @f$ d^{2}N_{ch}/d\eta\,d\phi@f$
 197    *
 198    * @param name  Name of histogram
 199    * @param title Title of histogram
 200    * @param a1    Eta axis
 201    * @param a2    Phi axis
 202    *
 203    * @return
 204    */
 205   TH2D* Make2D(const char* name, const char* title,
 206                const TAxis* a1, const TAxis* a2)
 207   {
 208     TH2D* ret = new TH2D(name, title,
 209                          a1->GetNbins(), a1->GetXmin(), a1->GetXmax(),
 210                          a2->GetNbins(), a2->GetXmin(), a2->GetXmax());
 211     ret->SetXTitle("#eta");
 212     ret->SetYTitle("#varphi");
 213     ret->SetZTitle("#frac{1}{N}#frac{dN^{2}_{ch}}{d#etad#varphi}");
 214     ret->Sumw2();
 215     ret->SetStats(0);
 216     ret->SetDirectory(0);
 217
 218     return ret;
 219   }
 220   //__________________________________________________________________
 221   /**
 222    * Utility function to set-up histograms based on the input
 223    * @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram.   This member function
 224    * is called on the first event so that we have the proper binning
 225    *
 226    * @param templ Input histogram
 227    *
 228    * @return true on succcess
 229    */
 230   Bool_t FirstEvent(const TH2D& templ)
 231   {
 232     const TAxis* etaAxis = templ.GetXaxis();
 233     const TAxis* phiAxis = templ.GetYaxis();
 234
 235     // Generate sum histograms.
 236     // - fTotal1D will be the sum of projections on the X axis of the input
 237     //   histograms
 238     // - fTotal2D will be the direct sum of the input histograms.
 239     // - fFromVtx will be the sum of the per vertex histograms after
 240     //   processing all events
 241     fTotal1D = Make1D("dndeta",
 242                       "#frac{1}{N}#frac{dN_{ch}}{d#eta} direct sum", etaAxis);
 243     fTotal2D = Make2D("d2ndetadphi",
 244                       "#frac{1}{N}#frac{dN^{2}_{ch}}{d#etad#phi} direct sum",
 245                       etaAxis, phiAxis);
 246     fFromVtx = Make1D("dndeta_test", "#frac{1}{N}#frac{dN_{ch}}{d#eta} "
 247                       "from VTX", etaAxis);
 248     fTotal1D->SetMarkerStyle(25);
 249     fTotal1D->SetMarkerSize(1.1);
 250     fNorm    = new TH1D("norm", "Normalisation",
 251                         etaAxis->GetNbins(),
 252                         etaAxis->GetXmin(),
 253                         etaAxis->GetXmax());
 254     fNorm->SetFillColor(kRed);
 255     fNorm->SetFillStyle(3001);
 256     fNorm->SetXTitle("#eta");
 257     fNorm->SetYTitle("Normalisation");
 258     fNorm->Sumw2();
 259     fNorm->SetDirectory(0);
 260     fVtx     = new TH1D("vtx", "Events per vertex bin",
 261                         fEvents->GetXaxis()->GetNbins(),
 262                         fEvents->GetXaxis()->GetXmin(),
 263                         fEvents->GetXaxis()->GetXmax());
 264     fVtx->SetXTitle("v_{z} [cm]");
 265     fVtx->SetYTitle("Events");
 266     fVtx->SetDirectory(0);
 267     fVtx->SetFillColor(kRed+1);
 268     fVtx->SetFillStyle(3001);
 269
 270
 271     // Generate the per-vertex bin histograms.  These will be the sum
 272     // of each of the per-event input histograms for a given vertex range.
 273     for (Int_t i = fBinVzMin; i <= fBinVzMax; i++) {
 274       TH1* h1 = Make1D(Form("dndeta_vz%02d", i),
 275                        Form("#frac{1}{N}#frac{dN_{ch}}{d#eta}, vtxbin=%d", i),
 276                        etaAxis);
 277       fByVtx1D.AddAtAndExpand(h1, i);
 278     }
 279     return kTRUE;
 280   }
 281
 282
 283   //__________________________________________________________________
 284   /**
 285    * Process the events
 286    *
 287    *
 288    * @return true on success, false otherwise
 289    */
 290   Bool_t Process()
 291   {
 292     // Get the number of events in the tree
 293     Int_t nEntries  = fTree->GetEntries();
 294     fNAccepted = 0;
 295
 296     // Loop over the events in the tree
 297     for (Int_t event = 0; event < nEntries; event++) {
 298       fTree->GetEntry(event);
 299
 300       // Get our input histogram
 301       const TH2D& hist = fAOD->GetHistogram();
 302
 303
 304       // If fTotal1D or fTotal2D are not made yet, do so (first event)
 305       if (!fTotal2D || !fTotal1D) {
 306         if (!FirstEvent(hist)) {
 307           Error("Process", "Failed to initialize on first event");
 308           return kFALSE;
 309         }
 310       }
 311
 312       // Check the trigger
 313       if (!fAOD->IsTriggerBits(AliAODForwardMult::kInel)) {
 314         // Info("Process", "Not an INEL event");
 315         continue;
 316       }
 317
 318       // Get the vertex bin - add 1 as we are using histogram bin
 319       // numbers in this class
 320       Int_t vtxBin = fAOD->GetVtxBin()+1;
 321
 322       // Check if we're within vertex cut
 323       if (!IsInsideVtxCut(vtxBin)) {
 324         Info("Process", "In event # %d, %d not within vertex cut "
 325              "[%d,%d] (%f,%f)",
 326              event, vtxBin, fBinVzMin, fBinVzMax, fVzMin, fVzMax);
 327         continue;
 328       }
 329       fVtx->AddBinContent(vtxBin);
 330
 331       // Increment our accepted counter
 332       fNAccepted++;
 333
 334       // Get the scale of this vertex range
 335       Double_t vtxScale = fEvents->GetBinContent(vtxBin);
 336
 337       // Get 'by-vertex' histograms
 338       TH1* tmp1 = static_cast<TH1D*>(fByVtx1D.At(vtxBin));
 339       if (!tmp1) {
 340         Warning("Process", "No histogram for vertex bin %d)",vtxBin);
 341         continue;
 342       }
 343
 344       Double_t scale = 1. / vtxScale;
 345       if (tmp1->InheritsFrom(TProfile::Class()))
 346         scale = 1;
 347
 348       if (!AddContrib(hist, *tmp1, scale)) return kFALSE;
 349
 350     }
 351     for (Int_t iVz = fBinVzMin; iVz <= fBinVzMax; iVz++) {
 352       TList* l = static_cast<TList*>(fCollect
 353                                      ->FindObject(Form("vtxbin%02d",iVz-1)));
 354       if (!l) {
 355         Error("Process", "List vtxbin%02d not found in %s",
 356               iVz-1, fCollect->GetName());
 357         continue;
 358       }
 359       TH1F* ve = static_cast<TH1F*>(l->FindObject("etaAcceptance"));
 360       if (!ve){
 361         Error("Process", "No eta acceptance histogram found in  "
 362               "vtxbin%02d/etaAcceptance", iVz-1);
 363         continue;
 364       }
 365       ve->Sumw2();
 366       fNorm->Add(ve, fVtx->GetBinContent(iVz));
 367     }
 368     for (Int_t i = 1; i <= fNorm->GetNbinsX(); i++)
 369       fNorm->SetBinError(i, 0);
 370     // fNorm->Scale(1.);
 371
 372     Info("Process", "Accepted %6d out of %6d events", fNAccepted, nEntries);
 373     return kTRUE;
 374   }
 375   //__________________________________________________________________
 376   Bool_t AddContrib(const TH2D& hist, TH1& vtx1D, Double_t scale)
 377   {
 378     // Make a projection on the X axis of the input histogram
 379     TH1D*    proj  = hist.ProjectionX("_px", 0, -1, "e");
 380
 381     // Add to per-vertex histogram
 382     vtx1D.Add(proj); // , scale);
 383
 384 #if defined(USE_NORM_HIST)
 385     scale = 1;
 386 #endif
 387
 388 #if defined(USE_WEIGHTED_MEAN)
 389     Int_t nBinPhi = hist.GetXaxis()->GetNbins();
 390     for (Int_t binEta = 1; binEta <= nBinEta; binEta++) {
 391       AddBinContrib(binEta, *proj, *fTotal1D);
 392       for (Int_t binPhi = 1; binPhi < nBinPhi; binPhi++)
 393         AddBinContrib(hist.GetBin(binEta, binPhi), hist, *fTotal2D);
 394     }
 395 #else
 396     // Add to 1D summed histogram
 397     fTotal1D->Add(proj); // , scale);
 398
 399     // Add to 2D summed histogram
 400     fTotal2D->Add(&hist); // , scale);
 401 #endif
 402
 403     // Remove the projection
 404     delete proj;
 405
 406     return kTRUE;
 407   }
 408   //__________________________________________________________________
 409   Bool_t AddBinContrib(Int_t bin, const TH1& in, TH1& out)
 410   {
 411     Double_t ic = in.GetBinContent(bin);
 412     Double_t ie = in.GetBinError(bin);
 413     if (ie <= 0.00001) return kTRUE;
 414
 415     Double_t iw = 1 / (ie*ie);
 416     Double_t oc = out.GetBinContent(bin);
 417     Double_t oe = out.GetBinError(bin);
 418
 419     // Store weighted sum in histogram
 420     out.SetBinContent(bin, oc + iw * ic);
 421     // Store sum of weights in histogram
 422     out.SetBinError(bin, oe + iw);
 423
 424     return kTRUE;
 425   }
 426   //__________________________________________________________________
 427   Bool_t SetResult(TH1*)
 428   {
 429 #if defined(USE_WEIGHTED_MEAN)
 430     Double_t scale   = 1.;
 431     Int_t    nBinEta = fTotal2D->GetXaxis()->GetNbins();
 432     Int_t    nBinPhi = fTotal1D->GetXaxis()->GetNbins();
 433     for (Int_t binEta = 1; binEta <= nBinEta; binEta++) {
 434       SetBinResult(binEta, *fTotal1D);
 435       for (Int_t binPhi = 1; binPhi < nBinPhi; binPhi++)
 436         SetBinResult(fTotal2D->GetBin(binEta, binPhi), *fTotal2D);
 437     }
 438 #elif defined (USE_NORM_HIST)
 439     fTotal1D->Divide(fNorm);
 440     Int_t    nBinEta = fTotal2D->GetXaxis()->GetNbins();
 441     Int_t    nBinPhi = fTotal1D->GetXaxis()->GetNbins();
 442     for (Int_t binEta = 1; binEta <= nBinEta; binEta++) {
 443       Double_t n = fNorm->GetBinContent(binEta);
 444       for (Int_t binPhi = 1; binPhi < nBinPhi; binPhi++) {
 445         if (n <= 0) {
 446           fTotal2D->SetBinContent(binEta,binPhi,0);
 447           continue;
 448         }
 449         Double_t c = fTotal2D->GetBinContent(binEta,binPhi);
 450         fTotal2D->SetBinContent(binEta,binPhi,c/n);
 451       }
 452     }
 453 #else
 454 #endif
 455     // fTotal2D->Scale(1, "width");
 456     return kTRUE;
 457   }
 458
 459   //__________________________________________________________________
 460   Bool_t SetBinResult(Int_t bin, TH1& in)
 461   {
 462     Double_t ic = in.GetBinContent(bin);
 463     Double_t ie = in.GetBinError(bin);
 464     if (ie <= 0.00001) return kTRUE;
 465
 466     Double_t av = ic / ie;
 467     Double_t er = TMath::Sqrt(1/ie);
 468
 469     // Set bin to be
 470     // @f[
 471     //    \frac{\sum_i w_i c_i}{\sum_i w_i}
 472     // @f]
 473     // where @f$ w_i = 1/e_i^2@f$, and set the error to be
 474     // @f[
 475     //    \sqrt{\frac{1}{\sum_i w_i}}
 476     // @f]
 477     in.SetBinContent(bin, av);
 478     in.SetBinError(bin, er);
 479
 480     return kTRUE;
 481   }
 482   //__________________________________________________________________
 483   /**
 484    * Called at the end of the job.
 485    *
 486    * It plots the result of the analysis in tree canvases.
 487    * - One shows the per-vertex accumalted histograms and compares them
 488    *   to the per-vertex histograms made in the analysis.
 489    * - Another shows the final @f$ dN_{ch}/d\eta@f$ calculated in
 490    *   different ways and compare those to the  @f$ dN_{ch}/d\eta@f$ made
 491    *   during the analysis.
 492    * - The last canvas shows the @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram.
 493    *
 494    * @return true on succes, false otherwise
 495    */
 496   Bool_t Finish()
 497   {
 498     TFile* out = TFile::Open("out.root", "RECREATE");
 499
 500     // Set the style
 501     gStyle->SetOptTitle(0);
 502     gStyle->SetPadColor(0);
 503     gStyle->SetPadBorderSize(0);
 504     gStyle->SetPadBorderMode(0);
 505     gStyle->SetPadRightMargin(0.05);
 506     gStyle->SetPadTopMargin(0.05);
 507     gStyle->SetPalette(1);
 508
 509     // Get number of bins
 510     Int_t nBin = fBinVzMax-fBinVzMin+1;
 511
 512     // Make first canvas
 513     TCanvas* cVtx = new TCanvas("cVtx", "By Vertex", 1000, 700);
 514     cVtx->SetBorderSize(0);
 515     cVtx->SetBorderMode(0);
 516     cVtx->Divide((nBin+.5)/2, 2, 0.0005, 0.0005);
 517
 518     // Loop over vertex histograms
 519     Int_t nHists = 0;
 520     for (Int_t i = fBinVzMin; i <= fBinVzMax; i++) {
 521       TDirectory* vtxDir = out->mkdir(Form("vtx%02d", i));
 522       vtxDir->cd();
 523
 524       TH1D* vh1 = static_cast<TH1D*>(fByVtx1D.At(i));
 525       if (!vh1) {
 526         Error("Finish", "No histogram at %d", i);
 527         continue;
 528       }
 529
 530       fFromVtx->Add(vh1);
 531
 532       // Scale and add to output
 533       if (fVtx->GetBinContent(i) > 0)
 534         vh1->Scale(1. / fVtx->GetBinContent(i), "width");
 535
 536       // Write to output file
 537       vh1->Write();
 538
 539       // Divide the pad
 540       TVirtualPad* pad = cVtx->cd(i-fBinVzMin+1);
 541       pad->Divide(1,2,0.005,0);
 542
 543       // Draw the per-vertex histogram
 544       pad->cd(1);
 545       vh1->Draw();
 546
 547       // Get the same histogram from the analysis and draw it
 548       TProfile* p =
 549         static_cast<TProfile*>(fCollect
 550                                ->FindObject(Form("dndeta_vtx%02d",i-1)));
 551       p->SetMarkerColor(kBlue+1);
 552       p->SetLineColor(kBlue+1);
 553       p->SetMarkerSize(0.8);
 554       p->SetMarkerStyle(20);
 555       p->Draw("same");
 556       p->Write();
 557
 558       // Make the ratio of the two histograms and draw it
 559       pad->cd(2);
 560       TH1* ratio = static_cast<TH1*>(vh1->Clone(Form("ratio_vtx%02d", i-1)));
 561       // ratio->Scale(1./fEvents->GetBinContent(i));
 562       ratio->Divide(p);
 563       ratio->SetMarkerSize(.8);
 564       ratio->SetLineColor(kGray);
 565       ratio->Draw("P");
 566     }
 567     cVtx->cd();
 568     out->cd();
 569
 570     // Get the number of events selected
 571     // Int_t nEvSelected = fEvents->Integral(fBinVzMin, fBinVzMax);
 572
 573     // The second canvas
 574     TCanvas* cTotal1D =  new TCanvas("total", "Result", 800, 800);
 575     cTotal1D->SetFillColor(0);
 576     cTotal1D->SetBorderMode(0);
 577     cTotal1D->SetBorderSize(0);
 578     cTotal1D->cd();
 579
 580     // Make our main pad
 581     TPad* p1 = new TPad("p1", "p1", 0, 0.3, 1.0, 1.0, 0, 0);
 582     p1->Draw();
 583     p1->cd();
 584
 585     // Make a stack to 'auto-scale' when plotting more than 1 histogram.
 586     THStack* stack = new THStack("results", "Results for #frac{dN_{ch}{d#eta}");
 587
 588     // Scale the histogram summed over the vertex bins.  This must be
 589     // divided by the number of bins we have summed over.  If we do
 590     // rebinning, we should scale it again.
 591     fFromVtx->Divide(fNorm);
 592     // fFromVtx->Scale(1.);
 593     fFromVtx->Scale(1, "width");
 594
 595
 596     if (fRebin > 1) { fFromVtx->Rebin(fRebin); fFromVtx->Scale(1. / fRebin); }
 597     stack->Add(fFromVtx, "");
 598
 599     // Generate the projection from the direct sum of the input histograms,
 600     // and scale it to the number of vertex bins and the bin width
 601     TH1* proj = fTotal2D->ProjectionX("dndeta_proj", 0, -1, "e");
 602     proj->SetLineColor(kGreen+1);
 603     proj->SetMarkerColor(kGreen+1);
 604     proj->SetMarkerStyle(24);
 605     proj->SetMarkerSize(1.5);
 606     proj->SetTitle(Form("%s directly", proj->GetTitle()));
 607     proj->Divide(fNorm);
 608     proj->Scale(1., "width");
 609     stack->Add(proj, "");
 610
 611     // Scale our direct sum of the projects of the input histograms to
 612     // the number of vertex bins and the bin width. If we do rebinning,
 613     // we must scale it one more time.
 614     // fTotal1D->Scale(1. / fNAccepted, "width");
 615     fTotal1D->Divide(fNorm);
 616     fTotal1D->Scale(1., "width");
 617     if (fRebin > 1) { fTotal1D->Rebin(fRebin); fTotal1D->Scale(1. / fRebin); }
 618     stack->Add(fTotal1D);
 619
 620     // Get the result from the analysis and plit that too (after modifying
 621     // the style and possible rebinning)
 622     TProfile* dtotal = static_cast<TProfile*>(fCollect->FindObject("dndeta"));
 623     if (!dtotal) {
 624       Error("Finish", "Couldn't get the event histogram");
 625       return kFALSE;
 626     }
 627     dtotal->SetTitle(Form("%s directly", dtotal->GetTitle()));
 628     dtotal->SetLineColor(kBlue+1);
 629     dtotal->SetMarkerColor(kBlue+1);
 630     dtotal->SetMarkerStyle(20);
 631     dtotal->SetMarkerSize(0.8);
 632     if (fRebin > 1) { dtotal->Rebin(fRebin); }
 633     stack->Add(dtotal, "");
 634
 635     // Draw stack
 636     stack->Draw("nostack e1");
 637     stack->Write();
 638
 639     // Make a legend
 640     TLegend* l = p1->BuildLegend(0.31, 0.15, 0.5, 0.6);
 641     l->SetFillColor(0);
 642     l->SetBorderSize(0);
 643     l->SetTextSize(0.04); // l->GetTextSize()/3);
 644     cTotal1D->cd();
 645
 646     // Generate our second pad
 647     TPad* p2 = new TPad("p2", "p2", 0, 0.0, 1.0, 0.3, 0, 0);
 648     p2->Draw();
 649     p2->cd();
 650
 651     // Create a new stack
 652     stack = new THStack("ratios", "Ratios to old method");
 653     // Calculate the ratio of our summed over vertex bins to the
 654     // result from the analysis and draw it
 655     TH1* ratioVtx = static_cast<TH1*>(fFromVtx->Clone("ratioVtx"));
 656     ratioVtx->SetDirectory(0);
 657     ratioVtx->Divide(dtotal);
 658     stack->Add(ratioVtx, "");
 659
 660     // Calculate the ratio of our direct sum of input histograms
 661     // result from the analysis  and draw it
 662     TH1* ratioSum = static_cast<TH1*>(fTotal1D->Clone("ratioSum"));
 663     ratioSum->SetDirectory(0);
 664     ratioSum->Divide(dtotal);
 665     stack->Add(ratioSum, "");
 666
 667     // Calculate the ratio of our direct sum of input histograms
 668     // result from the analysis  and draw it
 669     TH1* ratioProj = static_cast<TH1*>(proj->Clone("ratioProj"));
 670     ratioProj->SetDirectory(0);
 671     ratioProj->Divide(dtotal);
 672     stack->Add(ratioProj, "");
 673
 674     stack->Draw("nostack e1");
 675     stack->Write();
 676
 677     // update
 678     cTotal1D->cd();
 679
 680     // Generate the last canvas and show the summed input histogram
 681     TCanvas* other = new TCanvas("other", "Other", 800, 600);
 682     other->SetFillColor(0);
 683     other->SetBorderMode(0);
 684     other->SetBorderSize(0);
 685     other->Divide(2,1);
 686
 687     other->cd(1);
 688     fNorm->Draw("HIST");
 689     fNorm->Write();
 690
 691     other->cd(2);
 692     fEvents->Draw();
 693     fVtx->Draw("same");
 694     fEvents->Write();
 695     // fTotal2D->Draw("lego2 e");
 696
 697     return kTRUE;
 698   }
 699   //__________________________________________________________________
 700   /**
 701    * Run the post-processing.
 702    *
 703    * This will open the files &lt;<i>base</i>&gt;<tt>_hists.root</tt>
 704    * containing the histograms generated by AliForwardMultiplicity and
 705    * the file &lt;<i>base</i>&gt;<tt>_aods.root</tt> containing the
 706    * tree with AliAODEvent objects.
 707    *
 708    * Then it will loop over the events, accepting only INEL events
 709    * that have a primary collision point along z within @a vzMin and
 710    * @a vzMax.
 711    *
 712    * After the processing, the result will be shown in 3 canvases,
 713    * possibly  rebinning the result by the factor @a rebin.
 714    *
 715    * @param base  Base name of files
 716    * @param vzMin Minimum collision vertex z position to use
 717    * @param vzMax Maximum collision vertex z position to use
 718    * @param rebin Rebinning factor
 719    *
 720    * @return true on success, false otherwise
 721    */
 722   Bool_t Run(const char* base,
 723              Double_t    vzMin=-10,
 724              Double_t    vzMax=10,
 725              Int_t       rebin=1)
 726   {
 727     if (!Open(base,vzMin,vzMax,rebin)) return kFALSE;
 728     if (!Process()) return kFALSE;
 729     if (!Finish()) return kFALSE;
 730
 731     return kTRUE;
 732   }
 733 protected:
 734   Double_t           fVzMin;     // Minimum vertex
 735   Double_t           fVzMax;     // Maximum vertex
 736   Int_t              fBinVzMin;  // Corresponding bin to min vertex
 737   Int_t              fBinVzMax;  // Corresponding bin to max vertex
 738   Int_t              fRebin;     // Rebin factor
 739   TTree*             fTree;      // Event tree
 740   AliAODForwardMult* fAOD;       // Our event-by-event data
 741   TH1I*              fEvents;    // histogram of event counts per vertex
 742   TList*             fCollect;   // List of analysis histograms
 743   TObjArray          fByVtx1D;   // List of per-vertex 1D histograms
 744   TH1D*              fTotal1D;   // Direct sum of input histograms
 745   TH2D*              fTotal2D;   // Direct sum of input histograms
 746   TH1D*              fFromVtx;   // Sum of per-vertex histograms
 747   TH1D*              fNorm;      // Histogram of # events with data per bin
 748   TH1D*              fVtx;
 749   Int_t              fNAccepted; // # of accepted events
 750 };
 751
 752 //
 753 // EOF
 754 //