Various fixes to deal with centrality
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 20 Mar 2011 16:17:03 +0000 (16:17 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 20 Mar 2011 16:17:03 +0000 (16:17 +0000)
DrawdNdeta.C can now draw centrality specific results and compare to other data (MC truth, other results, etc.)
Reduced the total number of events to normalise to to a simpler expresion.

16 files changed:
PWG2/FORWARD/analysis2/AddTaskCentraldNdeta.C
PWG2/FORWARD/analysis2/AddTaskForwarddNdeta.C
PWG2/FORWARD/analysis2/AliBasedNdetaTask.cxx
PWG2/FORWARD/analysis2/AliBasedNdetaTask.h
PWG2/FORWARD/analysis2/AliFMDDensityCalculator.cxx
PWG2/FORWARD/analysis2/AliFMDDensityCalculator.h
PWG2/FORWARD/analysis2/AliFMDEventInspector.cxx
PWG2/FORWARD/analysis2/AliFMDEventInspector.h
PWG2/FORWARD/analysis2/AliFMDMCEventInspector.cxx
PWG2/FORWARD/analysis2/AliFMDMCEventInspector.h
PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx
PWG2/FORWARD/analysis2/DrawdNdeta.C
PWG2/FORWARD/analysis2/MakeAOD.C
PWG2/FORWARD/analysis2/MakedNdeta.C
PWG2/FORWARD/analysis2/OtherData.C
PWG2/FORWARD/analysis2/Run.sh

index 116ebfa..712e3d7 100644 (file)
@@ -9,11 +9,11 @@
  * 
  */
 AliAnalysisTask*
-AddTaskCentraldNdeta(const char* trig="INEL", 
-                    Double_t vzMin=-10, 
-                    Double_t vzMax=10, 
-                    Float_t centLow=0, 
-                    Float_t centHigh=100)
+AddTaskCentraldNdeta(const char* trig     = "INEL", 
+                    Double_t    vzMin    = -10, 
+                    Double_t    vzMax    = +10, 
+                    Bool_t      useCent  = false,
+                    Bool_t      cutEdges = false)
 {
   // analysis manager
   AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
@@ -23,16 +23,12 @@ AddTaskCentraldNdeta(const char* trig="INEL",
   AliCentraldNdetaTask* task = new AliCentraldNdetaTask("Central");
   task->SetVertexRange(vzMin, vzMax);
   task->SetTriggerMask(trig);
-  // task->SetUseShapeCorrection(false);
-  task->AddCentralityBin( 0,  0); // All bin - integrate over centrality
-  task->AddCentralityBin( 0,  5);
-  task->AddCentralityBin( 5, 10);
-  task->AddCentralityBin(10, 20);
-  task->AddCentralityBin(20, 30);
-  task->AddCentralityBin(30, 40);
-  task->AddCentralityBin(40, 50);
-  task->AddCentralityBin(50, 60);
-  task->AddCentralityBin(60,100);
+  task->SetCutEdges(cutEdges);
+  task->SetUseShapeCorrection(false);
+  if (useCent) {
+    Short_t bins[] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
+    task->SetCentralityAxis(11, bins);
+  }
   mgr->AddTask(task);
 
   // create containers for input/output
index 3a41490..f4f290f 100644 (file)
@@ -9,12 +9,11 @@
  * 
  */
 AliAnalysisTask*
-AddTaskForwarddNdeta(const char* trig="INEL", 
-                    Double_t vzMin=-10, 
-                    Double_t vzMax=10, 
-                    Float_t centlow = 0, 
-                    Float_t centhigh = 100, 
-                    Bool_t cutEdges = false)
+AddTaskForwarddNdeta(const char* trig     = "INEL", 
+                    Double_t    vzMin    = -10, 
+                    Double_t    vzMax    = +10, 
+                    Bool_t      useCent  = false,
+                    Bool_t      cutEdges = false)
 {
   // analysis manager
   AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
@@ -25,16 +24,11 @@ AddTaskForwarddNdeta(const char* trig="INEL",
   task->SetVertexRange(vzMin, vzMax);
   task->SetTriggerMask(trig);
   task->SetCutEdges(cutEdges);
-  // task->SetUseShapeCorrection(false);
-  task->AddCentralityBin( 0,  0); // All bin - integrate over centrality
-  task->AddCentralityBin( 0,  5);
-  task->AddCentralityBin( 5, 10);
-  task->AddCentralityBin(10, 20);
-  task->AddCentralityBin(20, 30);
-  task->AddCentralityBin(30, 40);
-  task->AddCentralityBin(40, 50);
-  task->AddCentralityBin(50, 60);
-  task->AddCentralityBin(60,100);
+  task->SetUseShapeCorrection(false);
+  if (useCent) {
+    Short_t bins[] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
+    task->SetCentralityAxis(11, bins);
+  }
   mgr->AddTask(task);
 
   // create containers for input/output
index 3913591..d5e0c85 100644 (file)
@@ -32,7 +32,8 @@ AliBasedNdetaTask::AliBasedNdetaTask()
     fUseShapeCorr(true),
     fSNNString(0),
     fSysString(0),
-    fCent(0)
+    fCent(0),
+    fCentAxis(0)
 {
   // 
   // Constructor
@@ -57,12 +58,13 @@ AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
     fUseShapeCorr(true),
     fSNNString(0),
     fSysString(0),
-    fCent(0)
+    fCent(0),
+    fCentAxis(0)
 {
   // 
   // Constructor
   // 
-  fListOfCentralities = new TList;
+  fListOfCentralities = new TObjArray(1);
 
   // Output slot #1 writes into a TH1 container
   DefineOutput(1, TList::Class()); 
@@ -87,7 +89,8 @@ AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
     fUseShapeCorr(o.fUseShapeCorr),
     fSNNString(o.fSNNString),
     fSysString(o.fSysString),
-    fCent(o.fCent)
+    fCent(o.fCent),
+    fCentAxis(o.fCentAxis)
 {}
 
 //____________________________________________________________________
@@ -110,7 +113,22 @@ AliBasedNdetaTask::~AliBasedNdetaTask()
 
 //________________________________________________________________________
 void 
-AliBasedNdetaTask::AddCentralityBin(Short_t low, Short_t high)
+AliBasedNdetaTask::SetCentralityAxis(UShort_t n, Short_t* bins)
+{
+  if (!fCentAxis) { 
+    fCentAxis = new TAxis();
+    fCentAxis->SetName("centAxis");
+    fCentAxis->SetTitle("Centrality [%]");
+  }
+  TArrayD dbins(n+1);
+  for (UShort_t i = 0; i <= n; i++) 
+    dbins[i] = (bins[i] == 100 ? 100.1 : bins[i]);
+  fCentAxis->Set(n, dbins.GetArray());
+}
+    
+//________________________________________________________________________
+void 
+AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high)
 {
   // 
   // Add a centrality bin 
@@ -119,7 +137,9 @@ AliBasedNdetaTask::AddCentralityBin(Short_t low, Short_t high)
   //    low  Low cut
   //    high High cut
   //
-  fListOfCentralities->Add(MakeCentralityBin(GetName(), low, high));
+  CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
+  AliInfo(Form("Adding centrality bin %p [%3d,%3d] @ %d", bin, low, high, at));
+  fListOfCentralities->AddAtAndExpand(bin, at);
 }
 
 //________________________________________________________________________
@@ -198,7 +218,15 @@ AliBasedNdetaTask::UserCreateOutputObjects()
   fSums->SetOwner();
 
   // Automatically add 'all' centrality bin if nothing has been defined. 
-  if (fListOfCentralities->GetEntries() <= 0) AddCentralityBin(0,0); 
+  AddCentralityBin(0, 0, 0);
+  if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) { 
+    const TArrayD* bins = fCentAxis->GetXbins();
+    Int_t          nbin = fCentAxis->GetNbins(); 
+    for (Int_t i = 0; i < nbin; i++) 
+      AddCentralityBin(i+1,  Short_t((*bins)[i]), Short_t((*bins)[i+1]));
+  }
+  fSums->Add(fCentAxis);
+
 
   // Centrality histogram 
   fCent = new TH1D("cent", "Centrality", 100, 0, 100);
@@ -244,19 +272,28 @@ AliBasedNdetaTask::UserExec(Option_t *)
     return;
   }
   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
-
+  
   // Fill centrality histogram 
-  fCent->Fill(forward->GetCentrality());
+  Float_t cent = forward->GetCentrality();
+  fCent->Fill(cent);
 
   // Get the histogram(s) 
   TH2D* data   = GetHistogram(aod, false);
   TH2D* dataMC = GetHistogram(aod, true);
 
   // Loop over centrality bins 
-  TIter next(fListOfCentralities);
-  CentralityBin* bin = 0;
-  while ((bin = static_cast<CentralityBin*>(next()))) 
-    bin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
+  CentralityBin* allBin = static_cast<CentralityBin*>(fListOfCentralities->At(0));
+  allBin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
+  
+  // Find this centrality bin 
+  if (fCentAxis && fCentAxis->GetNbins() > 0) {
+    Int_t          icent   = fCentAxis->FindBin(cent);
+    CentralityBin* thisBin = 0;
+    if (icent >= 1 && icent <= fCentAxis->GetNbins()) 
+      thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
+    if (thisBin)
+      thisBin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
+  }
 
   // Here, we get the update 
   if (!fSNNString) { 
@@ -422,6 +459,7 @@ AliBasedNdetaTask::Terminate(Option_t *)
 
   fSNNString = static_cast<TNamed*>(fSums->FindObject("sNN"));
   fSysString = static_cast<TNamed*>(fSums->FindObject("sys"));
+  fCentAxis  = static_cast<TAxis*>(fSums->FindObject("centAxis"));
 
   if(fSysString && fSNNString && 
      fSysString->GetUniqueID() == AliForwardUtil::kPP)
@@ -442,6 +480,9 @@ AliBasedNdetaTask::Terminate(Option_t *)
   // Output collision system string 
   if (fSysString) fOutput->Add(fSysString->Clone());
 
+  // Output centrality axis 
+  if (fCentAxis) fOutput->Add(fCentAxis);
+
   // Output trigger string 
   TNamed* trigString = 
     new TNamed("trigString", AliAODForwardMult::GetTriggerString(fTriggerMask));
@@ -812,10 +853,10 @@ AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
   if (!forward) return false;
   
   // Check the centrality class unless this is the 'all' bin 
-  if (!IsAllBin()) { 
-    Double_t centrality = forward->GetCentrality();
-    if (centrality < fLow || centrality >= fHigh) return false;
-  }
+  // if (!IsAllBin()) { 
+  //   Double_t centrality = forward->GetCentrality();
+  //  if (centrality < fLow || centrality >= fHigh) return false;
+  // }
   
   fTriggers->AddBinContent(kAll);
   if (forward->IsTriggerBits(AliAODForwardMult::kB)) 
@@ -930,7 +971,7 @@ AliBasedNdetaTask::CentralityBin::End(TList*      sums,
     return;
   }
   if (!fSum) { 
-    AliError("Couldn't find histogram 'base' in list");
+    AliError(Form("Couldn't find histogram '%s' in list", GetName()));
     return;
   }
   
@@ -944,23 +985,50 @@ AliBasedNdetaTask::CentralityBin::End(TList*      sums,
   Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
   Int_t nAccepted   = Int_t(fTriggers->GetBinContent(kAccepted));
   Int_t nGood       = nB - nA - nC + 2 * nE;
+
   
-  Double_t alpha    = Double_t(nAccepted) / nWithVertex;
-  Double_t vNorm     = nAccepted + alpha*(nTriggered - nWithVertex);
-  Double_t vtxEff   = Double_t(nAccepted) / vNorm;
-  vtxEff            = vtxEff * trigEff;
-  
-  if (nTriggered <= 0) { 
+  if (nTriggered <= 0.1) { 
     AliError("Number of triggered events <= 0");
     return;
   }
-  if (nGood <= 0) { 
+  if (nWithVertex <= 0.1) { 
+    AliError("Number of events with vertex <= 0");
+    return;
+  }
+  if (nGood <= 0.1) { 
     AliWarning(Form("Number of good events=%d=%d-%d-%d+2*%d<=0",
                    nGood, nB, nA, nC, nE));
     nGood = nMB;
   }
+
+
+  // Scaling 
+  //  N_A + N_A/N_V (N_T - N_V) = N_A + N_A/N_V*N_T - N_A 
+  //                            = N_A/N_V * N_T 
+  // 
+  // where 
+  //    N_A = nAccepted 
+  //    N_V = nWithVertex 
+  //    N_T = nTriggered 
+  // so that the normalisation is simply 
+  //    N_A / E_V 
+  // where 
+  //    E_V=N_V/N_T is the vertex efficiency from data 
+  // 
+  // Double_t alpha    = Double_t(nAccepted) / nWithVertex;
+  // Double_t vNorm    = nAccepted + alpha*(nTriggered - nWithVertex);
+  AliInfo(Form("Calculating event normalisation as "
+              "1 / trigEff * nAccepted * nTriggered / nWithVertex "
+              "= 1/%f * %d * %d / %d = %f",
+              trigEff, nAccepted, nTriggered, nWithVertex,
+              1. / trigEff * nAccepted * double(nTriggered) / nWithVertex));
+  Double_t vNorm    = double(nWithVertex) / double(nTriggered);
+  Double_t vtxEff   = vNorm;
+  Double_t ntotal   = vtxEff * trigEff;
+  
+
   AliInfo(Form("Total of %9d events for %s\n"
-              "                   of these %9d are minimum bias\n"
+              "                   of these %9d are triggered minimum bias\n"
               "                   of these %9d has a %s trigger\n" 
               "                   of these %9d has a vertex\n" 
               "                   of these %9d were in [%+4.1f,%+4.1f]cm\n"
@@ -969,12 +1037,14 @@ AliBasedNdetaTask::CentralityBin::End(TList*      sums,
               "                     A|C = %9d (%9d+%-9d)\n"
               "                     E   = %9d\n"
               "                   Implies %9d good triggers\n"
-              "                   Vertex efficiency: %f (%f)",
+              "                   Vertex efficiency:      %f\n"
+              "                   Trigger efficiency:     %f\n"
+              "                   Total number of events: %f\n",
               nAll, GetTitle(), nMB, nTriggered, 
               AliAODForwardMult::GetTriggerString(triggerMask),
               nWithVertex, nAccepted,
               vzMin, vzMax, 
-              nB, nA+nC, nA, nC, nE, nGood, vtxEff, vNorm));
+              nB, nA+nC, nA, nC, nE, nGood, vtxEff, trigEff, ntotal));
   
   const char* name = GetName();
   
@@ -995,7 +1065,7 @@ AliBasedNdetaTask::CentralityBin::End(TList*      sums,
   dndeta->Divide(norm);
   
   // Scale by the vertex efficiency 
-  dndeta->Scale(vtxEff, "width");
+  dndeta->Scale(ntotal, "width");
   
   SetHistogramAttributes(dndeta, kRed+1, 20, Form("ALICE %s", name));
   SetHistogramAttributes(norm, kRed+1,20,Form("ALICE %s normalisation", name)); 
@@ -1033,6 +1103,76 @@ AliBasedNdetaTask::CentralityBin::End(TList*      sums,
     if (symmetrice)   
       fOutput->Add(Symmetrice(Rebin(dndetaMC, rebin, cutEdges)));
   }
+  if (!IsAllBin()) return;
+
+  // Temporary stuff 
+  TFile* forward = TFile::Open("forward.root", "READ");
+  if (!forward)  { 
+    AliWarning(Form("No forward.root file found"));
+    return;
+  }
+
+  TH1D* shapeCorrProj = 0;
+  if (shapeCorr) {
+    shapeCorrProj = static_cast<const TH2D*>(shapeCorr)->ProjectionX();
+    shapeCorrProj->Scale(1. / shapeCorr->GetNbinsY());
+    shapeCorrProj->SetDirectory(0);
+    fOutput->Add(shapeCorrProj);
+  }
+
+  TList* official = static_cast<TList*>(forward->Get("official"));
+  if (official) { 
+    TH1F*  histEta  = static_cast<TH1F*>(official->FindObject("fHistEta"));
+    if (histEta)  {
+      TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(), 
+                           histEta->GetXaxis()->GetXmin(), 
+                           histEta->GetXaxis()->GetXmax());
+      for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
+       oEta->SetBinContent(i, histEta->GetBinContent(i));
+       oEta->SetBinError(i, histEta->GetBinError(i));
+      }
+      if (shapeCorrProj) oEta->Divide(shapeCorrProj);
+      oEta->Scale(ntotal/nAccepted, "width");
+      oEta->SetDirectory(0);
+      oEta->SetMarkerStyle(21);
+      oEta->SetMarkerColor(kCyan+3);
+      fOutput->Add(oEta);
+      fOutput->Add(Rebin(oEta, rebin, false));
+    }
+    else 
+      AliWarning(Form("Couldn't find histogram fHistEta in list %s", 
+                     official->GetName()));
+  }
+  else 
+    AliWarning(Form("Couldn't find list 'official' in %s",forward->GetName()));
+
+  TList* tracks = static_cast<TList*>(forward->Get("tracks"));
+  if (tracks) { 
+    TH1F*  histEta  = static_cast<TH1F*>(tracks->FindObject("fHistEta"));
+    if (histEta)  {
+      TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(), 
+                           histEta->GetXaxis()->GetXmin(), 
+                           histEta->GetXaxis()->GetXmax());
+      for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
+       oEta->SetBinContent(i, histEta->GetBinContent(i));
+       oEta->SetBinError(i, histEta->GetBinError(i));
+      }
+      if (shapeCorrProj) oEta->Divide(shapeCorrProj);
+      oEta->Scale(ntotal/nAccepted, "width");
+      oEta->SetDirectory(0);
+      oEta->SetMarkerStyle(22);
+      oEta->SetMarkerColor(kMagenta+2);
+      fOutput->Add(oEta);
+      fOutput->Add(Rebin(oEta, rebin, false));
+    }
+    else 
+      AliWarning(Form("Couldn't find histogram fHistEta in list %s", 
+                     tracks->GetName()));
+  }
+  else 
+    AliWarning(Form("Couldn't find list 'tracks' in %s",forward->GetName()));
+
+  forward->Close();
 }
 
 //
index fd7bea7..7a6c2cd 100644 (file)
@@ -4,11 +4,13 @@
 #ifndef ALIBASEDNDETATASK_H
 #define ALIBASEDNDETATASK_H
 #include <AliAnalysisTaskSE.h>
+class TAxis;
 class TList;
 class TH2D;
 class TH1D;
 class AliAODEvent;
 class AliAODForwardMult;
+class TObjArray;
 
 /**
  * Task to determine the 
@@ -39,13 +41,6 @@ public:
    * @name Task configuration 
    */
   /** 
-   * Add a centrality bin 
-   * 
-   * @param low  Low cut
-   * @param high High cut
-   */
-  void AddCentralityBin(Short_t low, Short_t high);
-  /** 
    * Set the vertex range to use 
    * 
    * @param min Minimum (in centermeter)
@@ -71,6 +66,18 @@ public:
    */
   void SetTriggerMask(const char* mask);
   /** 
+   * Set the centrality bins to use. 
+   * 
+   * @code 
+   *   UShort_t bins[] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
+   *   task->SetCentralityBins(11, bins);
+   * @endcode 
+   * 
+   * @param n     Number of bins (elements in @a bins minus 1)
+   * @param bins  Bin limits 
+   */
+  void SetCentralityAxis(UShort_t n, Short_t* bins);
+  /** 
    * Trigger efficiency for selected trigger(s)
    * 
    * @param e Trigger efficiency 
@@ -199,6 +206,13 @@ protected:
    */
   virtual TH2D* GetHistogram(const AliAODEvent* aod, Bool_t mc=false) = 0;
   /** 
+   * Add a centrality bin 
+   * 
+   * @param low  Low cut
+   * @param high High cut
+   */
+  void AddCentralityBin(UShort_t at, Short_t low, Short_t high);
+  /** 
    * Make a centrality bin 
    * 
    * @param name  Name used for histograms
@@ -391,11 +405,12 @@ protected:
   Bool_t          fCorrEmpty;    // Correct for empty bins 
   Double_t        fTriggerEff;   // Trigger efficiency for selected trigger(s)
   TH1*            fShapeCorr;    // Shape correction 
-  TList*          fListOfCentralities; // Centrality bins 
+  TObjArray*      fListOfCentralities; // Centrality bins 
   Bool_t          fUseShapeCorr; // Whether to use shape correction
   TNamed*         fSNNString;    // sqrt(s_NN) string 
   TNamed*         fSysString;    // Collision system string 
   TH1D*           fCent;         // Centrality distribution 
+  TAxis*          fCentAxis;     // Centrality axis
 
   ClassDef(AliBasedNdetaTask,2); // Determine multiplicity in base area
 };
index 0a0dd8d..ae79799 100644 (file)
@@ -72,6 +72,7 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
   //    name Name of object
   //
   fRingHistos.SetName(GetName());
+  fRingHistos.SetOwner();
   fRingHistos.Add(new RingHistos(1, 'I'));
   fRingHistos.Add(new RingHistos(2, 'I'));
   fRingHistos.Add(new RingHistos(2, 'O'));
@@ -205,7 +206,10 @@ AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
   }
-  if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
+  if (idx < 0 || idx >= fRingHistos.GetEntries()) {
+    AliWarning(Form("Index %d of FMD%d%c out of range", idx, d, r));
+    return 0;
+  }
   
   return static_cast<RingHistos*>(fRingHistos.At(idx));
 }
@@ -256,9 +260,12 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
       UShort_t    nt= (q == 0 ? 512 : 256);
       TH2D*       h = hists.Get(d,r);
       RingHistos* rh= GetRingHistos(d,r);
-      rh->fEmptyStrips->Reset();
-      rh->fTotalStrips->Reset();
-      rh->fBasicHits->Reset();
+      if (!rh) { 
+       AliError(Form("No ring histogram found for FMD%d%c", d, r));
+       fRingHistos.ls();
+       return false;
+      }
+      rh->ResetPoissonHistos();
       
       for (UShort_t s=0; s<ns; s++) { 
        for (UShort_t t=0; t<nt; t++) {
@@ -854,61 +861,7 @@ AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
   fEmptyVsTotal->SetDirectory(0);
   fEmptyVsTotal->SetXTitle("Total # strips");
   fEmptyVsTotal->SetYTitle("Empty # strips");
-  fEmptyVsTotal->SetZTitle("Correlation");
-  
-  //Inserted by HHD
-  
-  //Float_t nbinlimitsFMD3I[8] =  {-3.4,-3.2,-3,-2.8,-2.6,-2.4,-2.2,-2};
-  Int_t nbins = 40;
-  /*if(fName.Contains("FMD1I")) nbins = 7;
-  if(fName.Contains("FMD2I")) nbins = 8;
-  if(fName.Contains("FMD2O")) nbins = 4;
-  if(fName.Contains("FMD3I")) nbins = 8;
-  if(fName.Contains("FMD3O")) nbins = 4;*/
-  Float_t lowlimit = -4, highlimit = 6;
-  /*if(fName.Contains("FMD1I")) {lowlimit = 3.6; highlimit = 5;}
-  if(fName.Contains("FMD2I")) {lowlimit = 2.2; highlimit = 3.8;}
-  if(fName.Contains("FMD2O")) {lowlimit = 1.6; highlimit = 2.4;} 
-  if(fName.Contains("FMD3I")) {lowlimit = -2.4; highlimit = -1.6;}
-  if(fName.Contains("FMD3O")) {lowlimit = -3.5; highlimit = -2.1;} 
-  
-  std::cout<<nbins<<"   "<<lowlimit<<"    "<<highlimit<<std::endl;
-  */
-  fTotalStrips = new TH2D(Form("total%s", fName.Data()), 
-                         Form("Total number of strips in %s", fName.Data()), 
-                         nbins, 
-                         lowlimit,
-                         highlimit, 
-                         (fRing == 'I' || fRing == 'i' ? 5 : 10), 
-                         0, 2*TMath::Pi());
-  fEmptyStrips = new TH2D(Form("empty%s", fName.Data()), 
-                         Form("Empty number of strips in %s", fName.Data()), 
-                         nbins, 
-                         lowlimit,
-                         highlimit, 
-                         (fRing == 'I' || fRing == 'i' ? 5 : 10), 
-                         0, 2*TMath::Pi());
-  fBasicHits   = new TH2D(Form("basichits%s", fName.Data()), 
-                         Form("Basic number of hits in %s", fName.Data()), 
-                         200, 
-                         -4, 
-                         6, 
-                         (fRing == 'I' || fRing == 'i' ? 20 : 40), 
-                         0, 2*TMath::Pi());
-  
-  fTotalStrips->SetDirectory(0);
-  fEmptyStrips->SetDirectory(0);
-  fBasicHits->SetDirectory(0);
-  fTotalStrips->SetXTitle("#eta");
-  fEmptyStrips->SetXTitle("#eta");
-  fBasicHits->SetXTitle("#eta");
-  fTotalStrips->SetYTitle("#varphi [radians]");
-  fEmptyStrips->SetYTitle("#varphi [radians]");
-  fBasicHits->SetYTitle("#varphi [radians]");
-  fTotalStrips->Sumw2();
-  fEmptyStrips->Sumw2();
-  fBasicHits->Sumw2();
-  
+  fEmptyVsTotal->SetZTitle("Correlation");  
 }
 //____________________________________________________________________
 AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
@@ -987,7 +940,71 @@ AliFMDDensityCalculator::RingHistos::~RingHistos()
   if (fEmptyStrips)    delete fEmptyStrips;
 }
 
-  //____________________________________________________________________
+//____________________________________________________________________
+void
+AliFMDDensityCalculator::RingHistos::ResetPoissonHistos()
+{
+  if (fTotalStrips) { 
+    fTotalStrips->Reset();
+    fEmptyStrips->Reset();
+    fBasicHits->Reset();
+    return;
+  }
+  //Inserted by HHD
+  
+  //Float_t nbinlimitsFMD3I[8] =  {-3.4,-3.2,-3,-2.8,-2.6,-2.4,-2.2,-2};
+  Int_t nbins = 40;
+  /*if(fName.Contains("FMD1I")) nbins = 7;
+  if(fName.Contains("FMD2I")) nbins = 8;
+  if(fName.Contains("FMD2O")) nbins = 4;
+  if(fName.Contains("FMD3I")) nbins = 8;
+  if(fName.Contains("FMD3O")) nbins = 4;*/
+  Float_t lowlimit = -4, highlimit = 6;
+  /*if(fName.Contains("FMD1I")) {lowlimit = 3.6; highlimit = 5;}
+  if(fName.Contains("FMD2I")) {lowlimit = 2.2; highlimit = 3.8;}
+  if(fName.Contains("FMD2O")) {lowlimit = 1.6; highlimit = 2.4;} 
+  if(fName.Contains("FMD3I")) {lowlimit = -2.4; highlimit = -1.6;}
+  if(fName.Contains("FMD3O")) {lowlimit = -3.5; highlimit = -2.1;} 
+  
+  std::cout<<nbins<<"   "<<lowlimit<<"    "<<highlimit<<std::endl;
+  */
+  fTotalStrips = new TH2D(Form("total%s", fName.Data()), 
+                         Form("Total number of strips in %s", fName.Data()), 
+                         nbins, 
+                         lowlimit,
+                         highlimit, 
+                         (fRing == 'I' || fRing == 'i' ? 5 : 10), 
+                         0, 2*TMath::Pi());
+  fEmptyStrips = new TH2D(Form("empty%s", fName.Data()), 
+                         Form("Empty number of strips in %s", fName.Data()), 
+                         nbins, 
+                         lowlimit,
+                         highlimit, 
+                         (fRing == 'I' || fRing == 'i' ? 5 : 10), 
+                         0, 2*TMath::Pi());
+  fBasicHits   = new TH2D(Form("basichits%s", fName.Data()), 
+                         Form("Basic number of hits in %s", fName.Data()), 
+                         200, 
+                         -4, 
+                         6, 
+                         (fRing == 'I' || fRing == 'i' ? 20 : 40), 
+                         0, 2*TMath::Pi());
+  
+  fTotalStrips->SetDirectory(0);
+  fEmptyStrips->SetDirectory(0);
+  fBasicHits->SetDirectory(0);
+  fTotalStrips->SetXTitle("#eta");
+  fEmptyStrips->SetXTitle("#eta");
+  fBasicHits->SetXTitle("#eta");
+  fTotalStrips->SetYTitle("#varphi [radians]");
+  fEmptyStrips->SetYTitle("#varphi [radians]");
+  fBasicHits->SetYTitle("#varphi [radians]");
+  fTotalStrips->Sumw2();
+  fEmptyStrips->Sumw2();
+  fBasicHits->Sumw2();
+}
+
+//____________________________________________________________________
 void
 AliFMDDensityCalculator::RingHistos::Init(const TAxis& /*eAxis*/)
 {
index 291200e..2f78db6 100644 (file)
@@ -289,6 +289,10 @@ protected:
      * @param nEvents Number of events 
      */
     void ScaleHistograms(TList* dir, Int_t nEvents);
+    /** 
+     * Create Poisson histograms 
+     */
+    void ResetPoissonHistos();
     TH2D*     fEvsN;         // Correlation of Eloss vs uncorrected Nch
     TH2D*     fEvsM;         // Correlation of Eloss vs corrected Nch
     TProfile* fEtaVsN;       // Average uncorrected Nch vs eta
@@ -296,13 +300,13 @@ protected:
     TProfile* fCorr;         // Average correction vs eta
     TH2D*     fDensity;      // Distribution inclusive Nch
     TH2D*     fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
-    TH2D*     fTotalStrips;  //! Total number of strips in a region
-    TH2D*     fEmptyStrips;  //! Total number of strips in a region
-    TH2D*     fBasicHits  ;  //! Total number basic hits in a region
+    TH2D*     fTotalStrips;  // Total number of strips in a region
+    TH2D*     fEmptyStrips;  // Total number of strips in a region
+    TH2D*     fBasicHits  ;  // Total number basic hits in a region
     TH2D*     fEmptyVsTotal; // # of empty strips vs total number of strips 
     
     
-    ClassDef(RingHistos,1);
+    ClassDef(RingHistos,3);
   };
   /** 
    * Get the ring histogram container 
index 26bfb74..d6b2390 100644 (file)
@@ -41,13 +41,15 @@ AliFMDEventInspector::AliFMDEventInspector()
     fHType(0),
     fHWords(0),
     fHCent(0),
+    fHCentVsQual(0),
     fLowFluxCut(1000),
     fMaxVzErr(0.2),
     fList(0),
     fEnergy(0),
     fField(999), 
     fCollisionSystem(kUnknown),
-    fDebug(0)
+    fDebug(0),
+    fCentAxis(0)
 {
   // 
   // Constructor 
@@ -63,13 +65,15 @@ AliFMDEventInspector::AliFMDEventInspector(const char* name)
     fHType(0),
     fHWords(0),
     fHCent(0),
+    fHCentVsQual(0),
     fLowFluxCut(1000),
     fMaxVzErr(0.2),
     fList(0),
     fEnergy(0),
     fField(999), 
     fCollisionSystem(kUnknown),
-    fDebug(0)
+    fDebug(0),
+    fCentAxis(0)
 {
   // 
   // Constructor 
@@ -88,13 +92,15 @@ AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o)
     fHType(o.fHType),
     fHWords(o.fHWords),
     fHCent(o.fHCent),
+    fHCentVsQual(o.fHCentVsQual),
     fLowFluxCut(o.fLowFluxCut),
     fMaxVzErr(o.fMaxVzErr),
     fList(o.fList),
     fEnergy(o.fEnergy),
     fField(o.fField), 
     fCollisionSystem(o.fCollisionSystem),
-    fDebug(0)
+    fDebug(0),
+    fCentAxis(0)
 {
   // 
   // Copy constructor 
@@ -110,12 +116,6 @@ AliFMDEventInspector::~AliFMDEventInspector()
   // 
   // Destructor 
   //
-  if (fHEventsTr)    delete fHEventsTr;
-  if (fHEventsTrVtx) delete fHEventsTrVtx;
-  if (fHTriggers)    delete fHTriggers;  
-  if (fHType)        delete fHType;
-  if (fHWords)       delete fHWords;
-  if (fHCent)        delete fHCent;
   if (fList)         delete fList;
 }
 //____________________________________________________________________
@@ -138,6 +138,7 @@ AliFMDEventInspector::operator=(const AliFMDEventInspector& o)
   fHType             = o.fHType;
   fHWords            = o.fHWords;
   fHCent             = o.fHCent;
+  fHCentVsQual       = o.fHCentVsQual;
   fLowFluxCut        = o.fLowFluxCut;
   fMaxVzErr          = o.fMaxVzErr;
   fDebug             = o.fDebug;
@@ -153,6 +154,7 @@ AliFMDEventInspector::operator=(const AliFMDEventInspector& o)
     if (fHType)        fList->Add(fHType);
     if (fHWords)       fList->Add(fHWords);
     if (fHCent)        fList->Add(fHCent);
+    if (fHCentVsQual)  fList->Add(fHCentVsQual);
   }
   return *this;
 }
@@ -199,6 +201,13 @@ AliFMDEventInspector::Init(const TAxis& vtxAxis)
   // Parameters:
   //   vtxAxis Vertex axis in use 
   //
+  
+  // -1.5 -0.5 0.5 1.5 ... 89.5 ... 100.5
+  // ----- 92 number --------- ---- 1 ---
+  TArrayD limits(93);
+  for (Int_t i = 0; i < 92; i++) limits[i] = -1.5 + i;
+  
+  fCentAxis  = new TAxis(limits.GetSize()-1, limits.GetArray());
   fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger", 
                        vtxAxis.GetNbins(), 
                        vtxAxis.GetXmin(), 
@@ -261,7 +270,7 @@ AliFMDEventInspector::Init(const TAxis& vtxAxis)
   fHWords->SetBit(TH1::kCanRebin);
   fList->Add(fHWords);
 
-  fHCent = new TH1F("cent", "Centrality", 101, -1.5, 100.5);
+  fHCent = new TH1F("cent", "Centrality", limits.GetSize()-1,limits.GetArray());
   fHCent->SetFillColor(kBlue+1);
   fHCent->SetFillStyle(3001);
   fHCent->SetStats(0);
@@ -269,6 +278,18 @@ AliFMDEventInspector::Init(const TAxis& vtxAxis)
   fHCent->SetXTitle("Centrality [%]");
   fHCent->SetYTitle("Events");
   fList->Add(fHCent);
+
+  fHCentVsQual = new TH2F("centVsQuality", "Quality vs Centrality", 
+                         5, 0, 5, limits.GetSize()-1, limits.GetArray());
+  fHCentVsQual->SetXTitle("Quality");
+  fHCentVsQual->SetYTitle("Centrality [%]");
+  fHCentVsQual->SetZTitle("Events");
+  fHCentVsQual->GetXaxis()->SetBinLabel(1, "OK");
+  fHCentVsQual->GetXaxis()->SetBinLabel(2, "Outside v_{z} cut");
+  fHCentVsQual->GetXaxis()->SetBinLabel(3, "V0 vs SPD outlier");
+  fHCentVsQual->GetXaxis()->SetBinLabel(4, "V0 vs TPC outlier");
+  fHCentVsQual->GetXaxis()->SetBinLabel(5, "V0 vs ZDC outlier");
+  fList->Add(fHCentVsQual);
 }
 
 //____________________________________________________________________
@@ -338,11 +359,18 @@ AliFMDEventInspector::Process(const AliESDEvent* event,
   fHType->Fill(lowFlux ? 0 : 1);
   
   // --- Read centrality information 
-  cent = -10;
-  if (!ReadCentrality(event, cent)) {
+  cent          = -10;
+  UShort_t qual = 0;
+  if (!ReadCentrality(event, cent, qual)) {
     if (fDebug > 3) 
       AliWarning("Failed to get centrality");
   }
+  fHCent->Fill(cent);
+  if (qual == 0) fHCentVsQual->Fill(0., cent);
+  else { 
+    for (UShort_t i = 0; i < 4; i++) 
+      if (qual & (1 << i)) fHCentVsQual->Fill(Double_t(i+1), cent);
+  }
 
   // --- Get the vertex information ----------------------------------
   vz          = 0;
@@ -380,7 +408,9 @@ AliFMDEventInspector::Process(const AliESDEvent* event,
 
 //____________________________________________________________________
 Bool_t
-AliFMDEventInspector::ReadCentrality(const AliESDEvent* esd, Double_t& cent)
+AliFMDEventInspector::ReadCentrality(const AliESDEvent* esd, 
+                                    Double_t& cent, 
+                                    UShort_t& qual) const
 {
   // 
   // Read centrality from event 
@@ -392,15 +422,16 @@ AliFMDEventInspector::ReadCentrality(const AliESDEvent* esd, Double_t& cent)
   // Return:
   //    False on error, true otherwise 
   //
+  cent = -1;
+  qual = 0;
   AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
-  if (centObj) {
-    // AliInfo(Form("Got centrality object %p with quality %d", 
-    //              centObj, centObj->GetQuality()));
-    // centObj->Print();
-    cent = centObj->GetCentralityPercentile("V0M");  
-  }
-  // AliInfo(Form("Centrality is %f", cent));
-  fHCent->Fill(cent);
+  if (!centObj)  return true;
+
+  // AliInfo(Form("Got centrality object %p with quality %d", 
+  //              centObj, centObj->GetQuality()));
+  // centObj->Print();
+  cent = centObj->GetCentralityPercentile("V0M");  
+  qual = centObj->GetQuality();
 
   return true;
 }
index aaa3c8a..64a784d 100644 (file)
@@ -9,6 +9,7 @@ class TH2D;
 class TH1D;
 class TH1I;
 class TH1F;
+class TH2F;
 class TAxis;
 class TList;
 
@@ -234,7 +235,8 @@ protected:
    * 
    * @return False on error, true otherwise 
    */
-  virtual Bool_t ReadCentrality(const AliESDEvent* esd, Double_t& cent);
+  virtual Bool_t ReadCentrality(const AliESDEvent* esd, Double_t& cent,
+                               UShort_t& qual) const;
 
   TH1I*    fHEventsTr;    //! Histogram of events w/trigger
   TH1I*    fHEventsTrVtx; //! Events w/trigger and vertex 
@@ -242,6 +244,7 @@ protected:
   TH1I*    fHType;        //! Type (low/high flux) of event
   TH1I*    fHWords;       //! Trigger words 
   TH1F*    fHCent;        //! Centrality 
+  TH2F*    fHCentVsQual;  //! Centrality vs quality 
   Int_t    fLowFluxCut;   //  Low flux cut
   Double_t fMaxVzErr;     //  Maximum error on v_z
   TList*   fList;         //! Histogram container 
@@ -249,7 +252,8 @@ protected:
   Short_t  fField;        // L3 magnetic field [kG]
   UShort_t fCollisionSystem; //  Collision system
   Int_t    fDebug;        //  Debug level 
-  ClassDef(AliFMDEventInspector,2); // Inspect the event 
+  TAxis*   fCentAxis;     // Centrality axis used in histograms
+  ClassDef(AliFMDEventInspector,3); // Inspect the event 
 };
 
 #endif
index 74a9b01..17f2fe3 100644 (file)
 #include "AliGenGeVSimEventHeader.h"
 #include "AliHeader.h"
 #include <TList.h>
+#include <TH2F.h>
+
 //====================================================================
 AliFMDMCEventInspector::AliFMDMCEventInspector()
   : AliFMDEventInspector(), 
     fHVertex(0),
     fHPhiR(0), 
-    fHB(0)
+    fHB(0),
+    fHBvsPart(0),
+    fHBvsBin(0),
+    fHBvsCent(0),
+    fHVzComp(0),
+    fHCentVsPart(0),
+    fHCentVsBin(0)
 {
   // 
   // Constructor 
@@ -47,7 +55,13 @@ AliFMDMCEventInspector::AliFMDMCEventInspector(const char* /* name */)
   : AliFMDEventInspector("fmdEventInspector"), 
     fHVertex(0),
     fHPhiR(0), 
-    fHB(0)
+    fHB(0),
+    fHBvsPart(0),
+    fHBvsBin(0),
+    fHBvsCent(0),
+    fHVzComp(0),
+    fHCentVsPart(0),
+    fHCentVsBin(0)
 {
   // 
   // Constructor 
@@ -62,7 +76,13 @@ AliFMDMCEventInspector::AliFMDMCEventInspector(const AliFMDMCEventInspector& o)
   : AliFMDEventInspector(o), 
     fHVertex(0),
     fHPhiR(0), 
-    fHB(0)
+    fHB(0),
+    fHBvsPart(0),
+    fHBvsBin(0),
+    fHBvsCent(0),
+    fHVzComp(0),
+    fHCentVsPart(0),
+    fHCentVsBin(0)
 {
   // 
   // Copy constructor 
@@ -108,10 +128,13 @@ AliFMDMCEventInspector::Init(const TAxis& vtxAxis)
   //
   AliFMDEventInspector::Init(vtxAxis);
 
-  fHVertex = new TH1F("vertex", "True vertex distribution", 
-                     vtxAxis.GetNbins(), 
-                     vtxAxis.GetXmin(), 
-                     vtxAxis.GetXmax());
+  Int_t    maxPart = 450;
+  Int_t    maxBin  = 225;
+  Int_t    maxB    = 25;
+  Int_t    nVtx = vtxAxis.GetNbins();
+  Double_t lVtx = vtxAxis.GetXmin();
+  Double_t hVtx = vtxAxis.GetXmax();
+  fHVertex = new TH1F("vertex", "True vertex distribution", nVtx, lVtx, hVtx);
   fHVertex->SetXTitle("v_{z} [cm]");
   fHVertex->SetYTitle("# of events");
   fHVertex->SetFillColor(kGreen+1);
@@ -128,14 +151,67 @@ AliFMDMCEventInspector::Init(const TAxis& vtxAxis)
   fHPhiR->SetDirectory(0);
   fList->Add(fHPhiR);
 
-  fHB = new TH1F("b", "Impact parameter", 125, 0, 25);
+  fHB = new TH1F("b", "Impact parameter", 5*maxB, 0, maxB);
   fHB->SetXTitle("b [fm]");
   fHB->SetYTitle("# of events");
   fHB->SetFillColor(kGreen+1);
   fHB->SetFillStyle(3001);
   fHB->SetDirectory(0);
   fList->Add(fHB);
+
+  fHBvsPart = new TH2F("bVsParticipants", "Impact parameter vs Participants",
+                      5*maxB, 0, maxB, maxPart, -.5, maxPart-.5);
+  fHBvsPart->SetXTitle("b [fm]");
+  fHBvsPart->SetYTitle("# of participants");
+  fHBvsPart->SetZTitle("Events");
+  fHBvsPart->SetDirectory(0);
+  fList->Add(fHBvsPart);
+
+  fHBvsBin = new TH2F("bVsBinary", "Impact parameter vs Binary Collisions",
+                      5*maxB, 0, maxB, maxBin, -.5, maxBin-.5);
+  fHBvsBin->SetXTitle("b [fm]");
+  fHBvsBin->SetYTitle("# of binary collisions");
+  fHBvsBin->SetZTitle("Events");
+  fHBvsBin->SetDirectory(0);
+  fList->Add(fHBvsBin);
   
+  fHBvsCent = new TH2F("bVsCentrality", "Impact parameter vs Centrality",
+                      5*maxB, 0, maxB, fCentAxis->GetNbins(), 
+                      fCentAxis->GetXbins()->GetArray());
+  fHBvsCent->SetXTitle("b [fm]");
+  fHBvsCent->SetYTitle("Centrality [%]");
+  fHBvsCent->SetZTitle("Event");
+  fHBvsCent->SetDirectory(0);
+  fList->Add(fHBvsCent);
+  
+  
+  fHVzComp = new TH2F("vzComparison", "v_{z} truth vs reconstructed",
+                     10*nVtx, lVtx, hVtx, 10*nVtx, lVtx, hVtx);
+  fHVzComp->SetXTitle("True v_{z} [cm]");
+  fHVzComp->SetYTitle("Reconstructed v_{z} [cm]");
+  fHVzComp->SetZTitle("Events");
+  fHVzComp->SetDirectory(0);
+  fList->Add(fHVzComp);
+
+  fHCentVsPart = new TH2F("centralityVsParticipans", 
+                         "# of participants vs Centrality",
+                         maxPart, -.5, maxPart-.5, fCentAxis->GetNbins(), 
+                         fCentAxis->GetXbins()->GetArray());
+  fHCentVsPart->SetXTitle("Participants");
+  fHCentVsPart->SetYTitle("Centrality [%]");
+  fHCentVsPart->SetZTitle("Event");
+  fHCentVsPart->SetDirectory(0);
+  fList->Add(fHCentVsPart);
+
+  fHCentVsBin = new TH2F("centralityVsBinary", 
+                        "# of binary collisions vs Centrality",
+                        maxBin, -.5, maxBin-.5, fCentAxis->GetNbins(), 
+                        fCentAxis->GetXbins()->GetArray());
+  fHCentVsBin->SetXTitle("Binary collisions");
+  fHCentVsBin->SetYTitle("Centrality [%]");
+  fHCentVsBin->SetZTitle("Event");
+  fHCentVsBin->SetDirectory(0);
+  fList->Add(fHCentVsBin);
 }
 
 //____________________________________________________________________
@@ -145,6 +221,8 @@ AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event,
                                  UShort_t&         ivz, 
                                  Double_t&         vz,
                                  Double_t&         b,
+                                 Int_t&            npart, 
+                                 Int_t&            nbin,
                                  Double_t&         phiR)
 {
   // 
@@ -172,57 +250,69 @@ AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event,
   //Assign MC only triggers : True NSD etc.
   AliHeader*               header          = event->Header();
   AliGenEventHeader*       genHeader       = header->GenEventHeader();
+  AliCollisionGeometry*    colGeometry     = 
+    dynamic_cast<AliCollisionGeometry*>(genHeader);
   AliGenPythiaEventHeader* pythiaHeader    = 
     dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
   AliGenDPMjetEventHeader* dpmHeader       = 
     dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
   AliGenGeVSimEventHeader* gevHeader       = 
     dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
-  AliGenHijingEventHeader* hijingHeader    = 
-    dynamic_cast<AliGenHijingEventHeader*>(genHeader);
   AliGenHerwigEventHeader* herwigHeader    = 
     dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
+  // AliGenHijingEventHeader* hijingHeader    = 
+  //   dynamic_cast<AliGenHijingEventHeader*>(genHeader);
   // AliGenHydjetEventHeader* hydjetHeader    = 
   //   dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
-  AliGenEposEventHeader*   eposHeader      = 
-    dynamic_cast<AliGenEposEventHeader*>(genHeader);
+  // AliGenEposEventHeader*   eposHeader      = 
+  //   dynamic_cast<AliGenEposEventHeader*>(genHeader);
   
   // Check if this is a single diffractive event 
-  Bool_t   sd = kFALSE;
-  Double_t phi = -1111;
-  b            = -1;
+  Bool_t   sd    = kFALSE;
+  Double_t phi   = -1111;
+  npart          = 0;
+  nbin           = 0;
+  b              = -1;
+  if (colGeometry) { 
+    b     = colGeometry->ImpactParameter();
+    phi   = colGeometry->ReactionPlaneAngle();
+    npart = (colGeometry->ProjectileParticipants() + 
+            colGeometry->TargetParticipants());
+    nbin  = colGeometry->NN();
+  }
   if(pythiaHeader) {
     Int_t pythiaType = pythiaHeader->ProcessType();
     if (pythiaType==92 || pythiaType==93) sd = kTRUE;
-    b = pythiaHeader->GetImpactParameter();
+    b     = pythiaHeader->GetImpactParameter();
+    npart = 2; // Always 2 protons
+    nbin  = 1; // Always 1 binary collision 
   }
-  if(dpmHeader) {
+  if(dpmHeader) { // Also an AliCollisionGeometry 
     Int_t processType = dpmHeader->ProcessType();
     if (processType == 5 || processType == 6)  sd = kTRUE;
-    b    = dpmHeader->ImpactParameter();
-    phi  = dpmHeader->ReactionPlaneAngle();
-
   }
   if (gevHeader) { 
     phi  = gevHeader->GetEventPlane();
   }
-  if (hijingHeader) { 
-    b    = hijingHeader->ImpactParameter();
-    phi  = hijingHeader->ReactionPlaneAngle();
-  }
   if (herwigHeader) {
     Int_t processType = herwigHeader->ProcessType();
     // This is a guess 
     if (processType == 5 || processType == 6)  sd = kTRUE;
+    npart = 2; // Always 2 protons
+    nbin  = 1; // Always 1 binary collision 
   }
+  // if (hijingHeader) { 
+  // b    = hijingHeader->ImpactParameter();
+  // phi  = hijingHeader->ReactionPlaneAngle();
+  // }
   // if (hydjetHeader) {
   //   b    = hydjetHeader->ImpactParameter();
   //   phi  = hydjetHeader->ReactionPlaneAngle();
   // }
-  if (eposHeader) {
-    b    = eposHeader->ImpactParameter();
-    phi  = eposHeader->ReactionPlaneAngle();
-  }
+  // if (eposHeader) {
+  //   b    = eposHeader->ImpactParameter();
+  // phi  = eposHeader->ReactionPlaneAngle();
+  // }
 
   // Normalize event plane angle to [0,2pi]
   if (phi <= -1111) phiR = -1;
@@ -248,6 +338,8 @@ AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event,
   fHVertex->Fill(vz);
   fHPhiR->Fill(phiR);
   fHB->Fill(b);
+  fHBvsPart->Fill(b, npart);
+  fHBvsBin->Fill(b, nbin);
 
   // Check for the vertex bin 
   ivz = fHEventsTr->GetXaxis()->FindBin(vz);
@@ -265,7 +357,9 @@ AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event,
 }
 //____________________________________________________________________
 Bool_t
-AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd, Double_t& cent)
+AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd, 
+                                      Double_t& cent, 
+                                      UShort_t& qual) const
 {
   // 
   // Read centrality from event 
@@ -277,23 +371,33 @@ AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd, Double_t& cent)
   // Return:
   //    False on error, true otherwise 
   //
+  cent = -1;
+  qual = 0;
   AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
-  if (centObj) {
-    // AliInfo(Form("Got centrality object %p with quality %d", 
-    //              centObj, centObj->GetQuality()));
-    // centObj->Print();
-    if (centObj->GetQuality() == 0x8) 
-      cent = centObj->GetCentralityPercentileUnchecked("V0M");  
-    else
-      cent = centObj->GetCentralityPercentile("V0M");        
-  }
-  // AliInfo(Form("Centrality is %f", cent));
-  fHCent->Fill(cent);
+  if (!centObj) return true;
+
+  qual = centObj->GetQuality();
+  if (qual == 0x8) // Ignore ZDC outliers 
+    cent = centObj->GetCentralityPercentileUnchecked("V0M");  
+  else
+    cent = centObj->GetCentralityPercentile("V0M");        
 
   return true;
 }
 
-  
+//____________________________________________________________________
+Bool_t
+AliFMDMCEventInspector::CompareResults(Double_t vz,    Double_t trueVz, 
+                                      Double_t cent,  Double_t b,
+                                      Int_t    npart, Int_t    nbin)
+{
+  fHVzComp->Fill(trueVz, vz);
+  fHBvsCent->Fill(b, cent);
+  fHCentVsPart->Fill(npart, cent);
+  fHCentVsBin->Fill(nbin, cent);
+
+  return true;
+}  
 //
 // EOF
 //
index 50812f3..2994988 100644 (file)
@@ -5,6 +5,7 @@
 #define ALIFMDMCEVENTINSPECTOR_H
 #include "AliFMDEventInspector.h"
 class AliMCEvent;
+class TH2F;
 
 /** 
  * This class inspects the event 
@@ -83,7 +84,22 @@ public:
                   UShort_t&         ivz, 
                   Double_t&         vz,
                   Double_t&         b,
+                  Int_t&            npart, 
+                  Int_t&            nbin,
                   Double_t&         phiR);
+  /** 
+   * Compare the result of analysing the ESD for 
+   * the inclusive charged particle density to analysing 
+   * MC truth 
+   * 
+   * @param esd 
+   * @param mc 
+   * 
+   * @return 
+   */
+  virtual Bool_t CompareResults(Double_t vz,    Double_t trueVz, 
+                               Double_t cent,  Double_t b,
+                               Int_t    npart, Int_t    nbin);
 protected:
   /** 
    * Read centrality from event 
@@ -93,12 +109,19 @@ protected:
    * 
    * @return False on error, true otherwise 
    */
-  virtual Bool_t ReadCentrality(const AliESDEvent* esd, Double_t& cent);
+  virtual Bool_t ReadCentrality(const AliESDEvent* esd, Double_t& cent,
+                               UShort_t& qual) const;
 
-  TH1F* fHVertex; // Histogram of vertex 
-  TH1F* fHPhiR;   // Histogram of event plane 
-  TH1F* fHB;      // Histogram of impact parameter 
-  ClassDef(AliFMDMCEventInspector,1); // Inspect the event 
+  TH1F* fHVertex;  // Histogram of vertex 
+  TH1F* fHPhiR;    // Histogram of event plane 
+  TH1F* fHB;       // Histogram of impact parameter 
+  TH2F* fHBvsPart; // Impact parameter vs # participants 
+  TH2F* fHBvsBin;  // Impact parameter vs # participants 
+  TH2F* fHBvsCent; // Impact parameter vs centrality
+  TH2F* fHVzComp;  // True vs reconstructed vz
+  TH2F* fHCentVsPart; // Centrality versus # participants 
+  TH2F* fHCentVsBin;  // Centrality versus # binary collisions 
+  ClassDef(AliFMDMCEventInspector,2); // Inspect the event 
 };
 
 #endif
index ec5a90a..639bc2f 100644 (file)
@@ -267,9 +267,12 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   Double_t vzMC     = 0;
   Double_t phiR     = 0;
   Double_t b        = 0;
+  Int_t    npart    = 0;
+  Int_t    nbin     = 0;
   // UInt_t   foundMC  = 
-  fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, phiR);
-  
+  fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, 
+                           npart, nbin, phiR);
+  fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
   
   //Store all events
   MarkEventForStore();
@@ -280,8 +283,14 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   //MarkEventForStore();
   // Set trigger bits, and mark this event for storage 
   fAODFMD.SetTriggerBits(triggers);
-  fMCAODFMD.SetTriggerBits(triggers);
+  fAODFMD.SetSNN(fEventInspector.GetEnergy());
+  fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
   fAODFMD.SetCentrality(cent);
+
+  fMCAODFMD.SetTriggerBits(triggers);
+  fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
+  fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
+  fMCAODFMD.SetCentrality(cent);
   
   //All events should be stored - HHD
   //if (isAccepted) MarkEventForStore();
index 937a566..6650af9 100644 (file)
@@ -195,29 +195,94 @@ struct dNdetaDrawer
    * 
    * @param filename  File containing the data 
    */
-  void Run(const char* filename="forward_dndeta.C") 
+  void Run(const char* filename="forward_dndeta.root") 
   {
-    Double_t max = 0;
+    Double_t max = 0, rmax=0, amax=0;
 
-    if (!Open(filename, max)) return;
-    Info("Run", "Got data from %s with maximum %f", filename, max);
+    gStyle->SetPalette(1);
 
-    // Create our stack of results
-    THStack* results = fResults; // StackResults(max);
+    // --- Open input file -------------------------------------------
+    TFile* file = TFile::Open(filename, "READ");
+    if (!file) { 
+      Error("Open", "Cannot open %s", filename);
+      return;
+    }
+    // --- Get forward list ------------------------------------------
+    TList* forward = static_cast<TList*>(file->Get("ForwardResults"));
+    if (!forward) { 
+      Error("Open", "Couldn't find list ForwardResults");
+      return;
+    }
+    // --- Get information on the run --------------------------------
+    FetchInformation(forward);
+    // --- Set the macro pathand load other data script --------------
+    gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
+                            gROOT->GetMacroPath()));
+    gROOT->LoadMacro("OtherData.C");
+
+    // --- Get the central results -----------------------------------
+    TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
+    if (!clusters) Warning("Open", "Couldn't find list CentralResults");
 
-    // Create our stack of other results 
-    TMultiGraph* other = 0;
-    if (fShowOthers || fShowAlice) other = StackOther(max);
+    // --- Make our containtes ---------------------------------------
+    fResults   = new THStack("results", "Results");
+    fRatios    = new THStack("ratios",  "Ratios");
+    fLeftRight = new THStack("asymmetry", "Left-right asymmetry");
+    fOthers    = new TMultiGraph();
     
-    Double_t smax = 0;
-    THStack* ratios = 0;
-    if (fShowRatios) ratios = StackRatios(other, smax);
+    // --- Loop over input data --------------------------------------
+    FetchResults(forward,  "Forward", max, rmax, amax);
+    FetchResults(clusters, "Central", max, rmax, amax);
 
-    Double_t amax = 0;
-    THStack* leftright = 0;
-    if (fShowLeftRight) leftright = StackLeftRight(amax);
+    // --- Get trigger information -----------------------------------
+    TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
+    if (sums) {
+      TList* all = static_cast<TList*>(sums->FindObject("all"));
+      if (all) {
+       fTriggers = FetchResult(all, "triggers");
+       if (!fTriggers) all->ls();
+      }
+      else  {
+       Warning("Run", "List all not found in ForwardSums");
+       sums->ls();
+      }
+    }
+    else { 
+      Warning("Run", "No ForwardSums directory found in %s", file->GetName());
+      file->ls();
+    }
+    
+    // --- Check our stacks ------------------------------------------
+    if (!fResults->GetHists() || 
+       fResults->GetHists()->GetEntries() <= 0) { 
+      Error("Run", "No histograms in result stack!");
+      return;
+    }
+    if (!fOthers->GetListOfGraphs() || 
+       fOthers->GetListOfGraphs()->GetEntries() <= 0) { 
+      Warning("Run", "No other data found - disabling that");
+      fShowOthers = false;
+    }
+    if (!fRatios->GetHists() || 
+       fRatios->GetHists()->GetEntries() <= 0) { 
+      Warning("Run", "No ratio data found - disabling that");
+      // fRatios->ls();
+      fShowRatios = false;
+    }
+    if (!fLeftRight->GetHists() || 
+       fLeftRight->GetHists()->GetEntries() <= 0) { 
+      Warning("Run", "No left/right data found - disabling that");
+      // fLeftRight->ls();
+      fShowLeftRight = false;
+    }
+    
+    // --- Close the input file --------------------------------------
+    file->Close();
+
+    
 
-    Plot(results, other, max, ratios, smax, leftright, amax);
+    // --- Plot results ----------------------------------------------
+    Plot(max, rmax, amax);
   }
 
   //__________________________________________________________________
@@ -236,6 +301,8 @@ struct dNdetaDrawer
       fSysString  = static_cast<TNamed*>(results->FindObject("sys"));
     if (!fVtxAxis)
       fVtxAxis    = static_cast<TAxis*>(results->FindObject("vtxAxis"));
+    if (!fCentAxis) 
+      fCentAxis   = static_cast<TAxis*>(results->FindObject("centAxis"));
     if (!fTrigString) fTrigString = new TNamed("trigString", "unknown");
     if (!fSNNString)  fSNNString  = new TNamed("sNN", "unknown");
     if (!fSysString)  fSysString  = new TNamed("sys", "unknown");
@@ -245,246 +312,258 @@ struct dNdetaDrawer
       fVtxAxis->SetTitle("v_{z} range unspecified");
     }
 
+    TString centTxt;
+    if (fCentAxis) { 
+      Int_t nCent = fCentAxis->GetNbins();
+      centTxt = Form("\n   Centrality: %d bins", nCent);
+      for (Int_t i = 0; i <= nCent; i++) 
+       centTxt.Append(Form("%c%d", i == 0 ? ' ' : ',', 
+                           fCentAxis->GetXbins()->At(i)));
+    }
     Info("FetchInformation", 
         "Initialized for\n"
         "   Trigger:    %s  (%d)\n"
         "   sqrt(sNN):  %s  (%dGeV)\n"
         "   System:     %s  (%d)\n"
-        "   Vz range:   %s  (%f,%f)",
+        "   Vz range:   %s  (%f,%f)%s",
         fTrigString->GetTitle(), fTrigString->GetUniqueID(), 
         fSNNString->GetTitle(),  fSNNString->GetUniqueID(), 
         fSysString->GetTitle(),  fSysString->GetUniqueID(), 
-        fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
+        fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax(),
+        centTxt.Data());
   }
+  //__________________________________________________________________
+  TMultiGraph* FetchOthers(UShort_t centLow, UShort_t centHigh)
+  {
+    TMultiGraph* thisOther = 0;
+    if (!fShowOthers && !fShowAlice) return 0;
+
+    Bool_t   onlya = (fShowOthers ? false : true);
+    UShort_t sys   = (fSysString  ? fSysString->GetUniqueID() : 0);
+    UShort_t trg   = (fTrigString ? fTrigString->GetUniqueID() : 0);
+    UShort_t snn   = (fSNNString  ? fSNNString->GetUniqueID() : 0);
+    Long_t   ret   = gROOT->ProcessLine(Form("GetData(%d,%d,%d,%d,%d,%d);",
+                                            sys,snn,trg,
+                                            centLow,centHigh,onlya));
+    if (!ret) { 
+      Warning("FetchResults", "No other data found for sys=%d, sNN=%d, "
+             "trigger=%d %d%%-%d%% central %s", 
+             sys, snn, trg, centLow, centHigh, 
+             onlya ? " (ALICE results)" : "all");
+      return 0;
+    }
+    thisOther = reinterpret_cast<TMultiGraph*>(ret);
     
+    return thisOther;
+  }
   //__________________________________________________________________
   /** 
-   * Open input file, and find data 
+   * Get the results from the top-level list 
    * 
-   * @param filename File name
-   * 
-   * @return true on success 
+   * @param list  List 
+   * @param name  name 
+   * @param max   On return, maximum of data 
+   * @param rmax  On return, maximum of ratios
+   * @param amax  On return, maximum of left-right comparisons
    */
-  Bool_t Open(const char* filename, Double_t& max)
+  void FetchResults(const TList* list, 
+                   const char*  name, 
+                   Double_t&    max,
+                   Double_t&    rmax,
+                   Double_t&    amax)
   {
-    TFile* file = TFile::Open(filename, "READ");
-    if (!file) { 
-      Error("Open", "Cannot open %s", filename);
-      return false;
+    UShort_t n = fCentAxis ? fCentAxis->GetNbins() : 0;
+    if (n == 0) {
+      TList* all = static_cast<TList*>(list->FindObject("all"));
+      if (!all)
+       Error("FetchResults", "Couldn't find list 'all' in %s", 
+             list->GetName());
+      else 
+       FetchResults(all, name, FetchOthers(0,0), -1, 0, max, rmax, amax);
+      return;
     }
     
-    TList* results = static_cast<TList*>(file->Get("ForwardResults"));
-    if (!results) { 
-      Error("Open", "Couldn't find list ForwardResults");
-      return false;
+    Int_t   nCol = gStyle->GetNumberOfColors();
+    for (UShort_t i = 0; i < n; i++) { 
+      UShort_t centLow  = fCentAxis->GetXbins()->At(i);
+      UShort_t centHigh = fCentAxis->GetXbins()->At(i+1);
+      TString  lname    = Form("cent%03d_%03d", centLow, centHigh);
+      TList*   thisCent = static_cast<TList*>(list->FindObject(lname));
+
+      Float_t fc   = (centLow+double(centHigh-centLow)/2) / 100;
+      Int_t   icol = TMath::Min(nCol-1,int(fc * nCol + .5));
+      Int_t   col  = gStyle->GetColorPalette(icol);
+      Info("FetchResults","Centrality %d-%d color index %d (=%d*%f) -> %d", 
+          centLow, centHigh, icol, nCol, fc, col);
+
+      TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
+      if (!thisCent)
+       Error("FetchResults", "Couldn't find list '%s' in %s", 
+             lname.Data(), list->GetName());
+      else 
+       FetchResults(thisCent, name, FetchOthers(centLow, centHigh), 
+                    col, centTxt.Data(), max, rmax, amax);
     }
-
-    // Get information on the run 
-    FetchInformation(results);
-
-    TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
-    if (!clusters) Warning("Open", "Couldn't find list CentralResults");
-
-    fResults   = new THStack("results", "Results");
-    FetchResults(fResults, results,  "Forward", max);
-    FetchResults(fResults, clusters, "Central", max);
-
-    TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
-    if (sums) fTriggers = FetchResult(sums, "triggers");
-    
-    file->Close();
-
-    return true;
+  } 
+  //__________________________________________________________________
+  void SetAttributes(TH1* h, Int_t color)
+  {
+    if (!h) return;
+    if (color < 0) return;
+    h->SetLineColor(color);
+    h->SetMarkerColor(color);
+    h->SetFillColor(color);
+  }
+  //__________________________________________________________________
+  void SetAttributes(TGraph* g, Int_t color)
+  {
+    if (!g) return;
+    if (color < 0) return;
+    g->SetLineColor(color);
+    g->SetMarkerColor(color);
+    g->SetFillColor(color);
   }
   //__________________________________________________________________
-  void FetchResults(THStack* stack, const TList* list, const char* name, 
-                   Double_t& max)
+  void ModifyTitle(TNamed* h, const char* centTxt)
+  {
+    if (!centTxt || !h) return;
+    h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
+  }
+
+  //__________________________________________________________________
+  /** 
+   * Fetch results for a particular centrality bin
+   * 
+   * @param list       List 
+   * @param name       Name 
+   * @param centLow    Low cut on centrality 
+   * @param centHigh   high cut on centrality 
+   * @param max        On return, data maximum
+   * @param rmax       On return, ratio maximum 
+   * @param amax       On return, left-right maximum 
+   */
+  void FetchResults(const TList* list, 
+                   const char*  name, 
+                   TMultiGraph* thisOther,
+                   Int_t        color,
+                   const char*  centTxt,
+                   Double_t&    max,
+                   Double_t&    rmax,
+                   Double_t&    amax)
   {
     TH1* dndeta      = FetchResult(list, Form("dndeta%s", name));
     TH1* dndetaMC    = FetchResult(list, Form("dndeta%sMC", name));
     TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
     TH1* dndetaSym   = 0;
     TH1* dndetaMCSym = 0;
+    TH1* tracks      = FetchResult(list, "tracks");
+    SetAttributes(dndeta,     color);
+    SetAttributes(dndetaMC,   color+2);
+    SetAttributes(dndetaTruth,color);
+    SetAttributes(dndetaSym,  color);
+    SetAttributes(dndetaMCSym,color+2);
+    SetAttributes(tracks,     color+3);
+    ModifyTitle(dndeta,     centTxt);
+    ModifyTitle(dndetaMC,   centTxt);
+    ModifyTitle(dndetaTruth,centTxt);
+    ModifyTitle(dndetaSym,  centTxt);
+    ModifyTitle(dndetaMCSym,centTxt);
+    ModifyTitle(tracks,     centTxt);
+      
 
-    max = TMath::Max(max, AddHistogram(stack, dndetaTruth, "e5 p"));
-    max = TMath::Max(max, AddHistogram(stack, dndetaMC,    "", dndetaMCSym));
-    max = TMath::Max(max, AddHistogram(stack, dndeta,      "", dndetaSym));
+    max = TMath::Max(max, AddHistogram(fResults, dndetaTruth, "e5 p"));
+    max = TMath::Max(max, AddHistogram(fResults, dndetaMC,    dndetaMCSym));
+    max = TMath::Max(max, AddHistogram(fResults, dndeta,      dndetaSym));
+    max = TMath::Max(max, AddHistogram(fResults, tracks));
     
-    Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f", 
-        dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
-
-    if (dndetaMCSym) fNumerators.Add(dndetaMCSym);
-    if (dndetaMC)    fNumerators.Add(dndetaMC);
-    if (dndetaSym)   fNumerators.Add(dndetaSym);
-    if (dndeta)      fNumerators.Add(dndeta);
-    if (dndetaTruth) fDenominators.Add(dndetaTruth);
-  }
-  //__________________________________________________________________
-  /** 
-   * Make a histogram stack of results 
-   * 
-   * @param max On return, the maximum value in the stack 
-   * 
-   * @return Newly allocated stack
-   */
-  TMultiGraph* StackOther(Double_t& max)
-  {
-    gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/OtherData.C");
-    Int_t    error = 0;
-    Bool_t   onlya = (fShowOthers ? false : true);
-    Int_t    trg   = (fTrigString ? fTrigString->GetUniqueID() : 0);
-    UShort_t snn   = (fSNNString  ? fSNNString->GetUniqueID() : 0);
-    Long_t   ret   = gROOT->ProcessLine(Form("GetData(%d,%d,%d);",
-                                            snn,trg,onlya));
-    if (error) { 
-      Error("StackOther", "Failed to execute GetData(%d,%d,%d)", 
-           snn, trg, onlya);
-      return 0;
-    }
-    if (!ret) { 
-      Warning("StackOther", "No other data found for sNN=%d, trigger=%d", 
-             snn, trg);
-      return 0;
-    }
-    TMultiGraph* other = reinterpret_cast<TMultiGraph*>(ret);
+    // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f", 
+    //      dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
 
-    TGraphAsymmErrors* o      = 0;
-    TIter              next(other->GetListOfGraphs());
-    while ((o = static_cast<TGraphAsymmErrors*>(next()))) {
-      max = TMath::Max(max, TMath::MaxElement(o->GetN(), o->GetY()));
-      fDenominators.Add(o);
+    if (fShowLeftRight) {
+      fLeftRight->Add(Asymmetry(dndeta,    amax));
+      fLeftRight->Add(Asymmetry(dndetaMC,  amax));
+      fLeftRight->Add(Asymmetry(tracks,    amax));
     }
 
-    return other;
-  }
-  //__________________________________________________________________
-  /** 
-   * Make a histogram stack of ratios of results to other data
-   * 
-   * @param max On return, the maximum value in the stack 
-   * 
-   * @return Newly allocated stack
-   */
-  THStack* StackRatios(TMultiGraph* others, Double_t& max) 
-  {
-    THStack* ratios = new THStack("ratios", "Ratios");
-
-    if (others) {
-      TObject* ua5_1  = 0;
-      TObject* ua5_2  = 0;
-      TObject* alice  = 0;
-      TObject* cms    = 0;
-      TObject* denom  = 0;
-      TIter    nextD(&fDenominators);
-      while ((denom = nextD())) {
-       TIter nextN(&fNumerators);
-       TObject* numer = 0;
-       while ((numer = nextN())) { 
-         ratios->Add(Ratio(numer, denom, max));
-       }
-       TString oName(denom->GetName());
-       oName.ToLower();
-       if (oName.Contains("ua5"))  { 
-         if (ua5_1) ua5_2 = denom; 
-         else       ua5_1 = denom; 
-       }
-       if (oName.Contains("alice")) alice = denom;
-       if (oName.Contains("cms"))   cms   = denom;
+    if (thisOther) {
+      TIter next(thisOther->GetListOfGraphs());
+      TGraph* g = 0;
+      while ((g = static_cast<TGraph*>(next()))) {
+       fRatios->Add(Ratio(dndeta,    g, rmax));
+       fRatios->Add(Ratio(dndetaSym, g, rmax));
+       SetAttributes(g, color);
+       ModifyTitle(g, centTxt);
+       if (!fOthers->GetListOfGraphs() || 
+           !fOthers->GetListOfGraphs()->FindObject(g->GetName()))
+         fOthers->Add(g);
       }
-      if (ua5_1 && alice) ratios->Add(Ratio(alice, ua5_1, max));
-      if (ua5_2 && alice) ratios->Add(Ratio(alice, ua5_2, max));
-      if (cms   && alice) ratios->Add(Ratio(alice, cms,   max));
+      // fOthers->Add(thisOther);
     }
-    // Check if we have ratios 
-    if (!ratios->GetHists() || 
-       (ratios->GetHists()->GetEntries() <= 0)) { 
-      delete ratios; 
-      ratios = 0; 
+    if (tracks) { 
+      if (!fRatios->GetHists()->FindObject(tracks->GetName()))
+       fRatios->Add(Ratio(dndeta, tracks, rmax));
     }
-    return ratios;
-  }
-  //__________________________________________________________________
-  /** 
-   * Make a histogram stack of the left-right asymmetry 
-   * 
-   * @param max On return, the maximum value in the stack 
-   * 
-   * @return Newly allocated stack
-   */
-  THStack* StackLeftRight(Double_t& max)
-  {
-    THStack* ret = new THStack("leftright", "Left-right asymmetry");
-    // ret->Add(Asymmetry(fForward,    max));
-    // ret->Add(Asymmetry(fForwardMC,  max));
 
-    if (!ret->GetHists() || 
-       (ret->GetHists()->GetEntries() <= 0)) { 
-      delete ret; 
-      ret = 0; 
+    if (dndetaTruth) { 
+      fRatios->Add(Ratio(dndeta,      dndetaTruth, rmax));
+      fRatios->Add(Ratio(dndetaSym,   dndetaTruth, rmax));
+      fRatios->Add(Ratio(dndetaMC,    dndetaTruth, rmax));
+      fRatios->Add(Ratio(dndetaMCSym, dndetaTruth, rmax));
     }
-    return ret;
   }
   //__________________________________________________________________
   /** 
    * Plot the results
-   * 
-   * @param results    Results
-   * @param others     Other data
    * @param max        Max value 
-   * @param ratios     Stack of ratios (optional)
    * @param rmax       Maximum diviation from 1 of ratios 
-   * @param leftright  Stack of left-right asymmetry (optional)         
    * @param amax       Maximum diviation from 1 of asymmetries 
    */
-  void Plot(THStack*     results,    
-           TMultiGraph* others, 
-           Double_t     max, 
-           THStack*     ratios,     
+  void Plot(Double_t     max, 
            Double_t     rmax,
-           THStack*     leftright, 
            Double_t     amax)
   {
     gStyle->SetOptTitle(0);
     gStyle->SetTitleFont(132, "xyz");
     gStyle->SetLabelFont(132, "xyz");
     
-    Int_t    h = 800;
-    Int_t    w = 800; // h / TMath::Sqrt(2);
-    if (!ratios) w *= 1.4;
-    if (!leftright) w *= 1.4;
+    Int_t    h  = 800;
+    Int_t    w  = 800; // h / TMath::Sqrt(2);
+    Double_t y1 = 0;
+    Double_t y2 = 0;
+    Double_t y3 = 0;
+    if (!fShowRatios)    w  *= 1.4;
+    else                 y1 =  0.3;
+    if (!fShowLeftRight) w  *= 1.4;
+    else { 
+      Double_t y11 = y1;
+      y1 = (y11 > 0.0001 ? 0.4 : 0.2);
+      y2 = (y11 > 0.0001 ? y11 : 0.2);
+    }
     TCanvas* c = new TCanvas("Results", "Results", w, h);
     c->SetFillColor(0);
     c->SetBorderSize(0);
     c->SetBorderMode(0);
 
-    Double_t y1 = 0;
-    Double_t y2 = 0;
-    Double_t y3 = 0;
-    if (ratios)    y1 = 0.3;
-    if (leftright) { 
-      if (y1 > 0.0001) {
-       y2 = 0.2;
-       y1 = 0.4;
-      }
-      else {
-       y1 = 0.2;
-       y2 = y1;
-      }
-    }
-    PlotResults(results, others, max, y1);
+#if 0
+    Info("Plot", "y1=%f, y2=%f, y3=%f extra: %s %s", 
+        y1, y2, y2, (fShowRatios ? "ratios" : ""), 
+        (fShowLeftRight ? "right/left" : ""));
+#endif
+    PlotResults(max, y1);
     c->cd();
 
-    PlotRatios(ratios, rmax, y2, y1);
+    PlotRatios(rmax, y2, y1);
     c->cd( );
 
-    PlotLeftRight(leftright, amax, y3, y2);
+    PlotLeftRight(amax, y3, y2);
     c->cd();
 
     
     Int_t   vMin = fVtxAxis->GetXmin();
     Int_t   vMax = fVtxAxis->GetXmax();    
     TString trg(fTrigString->GetTitle());
-    Int_t   nev  = fTriggers->GetBinContent(fTriggers->GetNbinsX());
+    Int_t   nev  = 0;
+    if (fTriggers) nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
     trg          = trg.Strip(TString::kBoth);
     TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev",
                      fSysString->GetTitle(), 
@@ -499,15 +578,54 @@ struct dNdetaDrawer
   }
   //__________________________________________________________________
   /** 
+   * Build main legend 
+   * 
+   * @param x1 
+   * @param y1 
+   * @param x2 
+   * @param y2 
+   */
+  void BuildLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
+  {
+    TLegend* l = new TLegend(x1,y1,x2,y2);
+    l->SetNColumns(2);
+    l->SetFillColor(0);
+    l->SetFillStyle(0);
+    l->SetBorderSize(0);
+    l->SetTextFont(132);
+
+    TIter    next(fResults->GetHists());
+    TObject* hist = 0;
+    while ((hist = next())) { 
+      TString n(hist->GetTitle());
+      if (n.Contains("mirrored")) continue;
+      l->AddEntry(hist, hist->GetTitle(), "pl");
+    }
+    TIter nexto(fOthers->GetListOfGraphs());
+    while ((hist = nexto())) { 
+      TString n(hist->GetTitle());
+      if (n.Contains("mirrored")) continue;
+      l->AddEntry(hist, hist->GetTitle(), "pl");
+    }
+    TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
+    d1->SetLineColor(kBlack);
+    d1->SetMarkerColor(kBlack);
+    d1->SetMarkerStyle(20);
+    TLegendEntry* d2 = l->AddEntry("d2", "Mirrored data", "lp");
+    d2->SetLineColor(kBlack);
+    d2->SetMarkerColor(kBlack);
+    d2->SetMarkerStyle(24);
+    
+    l->Draw();
+  }
+  //__________________________________________________________________
+  /** 
    * Plot the results
    *    
-   * @param results   Results
-   * @param others    Other data
    * @param max       Maximum 
    * @param yd        Bottom position of pad 
    */
-  void PlotResults(THStack* results, TMultiGraph* others, 
-                  Double_t max, Double_t yd) 
+  void PlotResults(Double_t max, Double_t yd) 
   {
     // Make a sub-pad for the result itself
     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
@@ -522,34 +640,36 @@ struct dNdetaDrawer
     p1->Draw();
     p1->cd();
     
-    Info("PlotResults", "Plotting results with max=%f", max);
-    results->ls();
-    results->SetMaximum(1.15*max);
-    results->SetMinimum(yd > 0.00001 ? -0.1 : 0);
+    // Info("PlotResults", "Plotting results with max=%f", max);
+    fResults->SetMaximum(1.15*max);
+    fResults->SetMinimum(yd > 0.00001 ? -0.1 : 0);
 
-    FixAxis(results, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
+    FixAxis(fResults, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
 
     p1->Clear();
-    results->DrawClone("nostack e1");
+    fResults->DrawClone("nostack e1");
 
-    fRangeParam->fSlave1Axis = results->GetXaxis();
+    fRangeParam->fSlave1Axis = fResults->GetXaxis();
     fRangeParam->fSlave1Pad  = p1;
 
     // Draw other data
-    if (others) {
+    if (fShowOthers || fShowAlice) {
       TGraphAsymmErrors* o      = 0;
-      TIter              nextG(others->GetListOfGraphs());
+      TIter              nextG(fOthers->GetListOfGraphs());
       while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
         o->DrawClone("same p");
     }
 
     // Make a legend in the main result pad
+    BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
+#if 0
     TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
     l->SetNColumns(2);
     l->SetFillColor(0);
     l->SetFillStyle(0);
     l->SetBorderSize(0);
     l->SetTextFont(132);
+#endif
 
     // Put a title on top
     TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
@@ -574,7 +694,8 @@ struct dNdetaDrawer
     tt->Draw();
 
     // Put number of accepted events on the plot
-    Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
+    Int_t nev = 0;
+    if (fTriggers) nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
     TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
     et->SetNDC();
     et->SetTextFont(132);
@@ -591,7 +712,7 @@ struct dNdetaDrawer
     }
     // results->Draw("nostack e1 same");
 
-    fRangeParam->fSlave1Axis = FindXAxis(p1, results->GetName());
+    fRangeParam->fSlave1Axis = FindXAxis(p1, fResults->GetName());
     fRangeParam->fSlave1Pad  = p1;
 
 
@@ -619,14 +740,14 @@ struct dNdetaDrawer
   /** 
    * Plot the ratios 
    * 
-   * @param ratios  Ratios to plot (if any)
    * @param max     Maximum diviation from 1 
    * @param y1      Lower y coordinate of pad
    * @param y2      Upper y coordinate of pad
    */
-  void PlotRatios(THStack* ratios, Double_t max, Double_t y1, Double_t y2) 
+  void PlotRatios(Double_t max, Double_t y1, Double_t y2) 
   {
-    if (!ratios) return;
+    if (!fShowRatios) return;
+
     bool isBottom = (y1 < 0.0001);
     Double_t yd = y2 - y1;
     // Make a sub-pad for the result itself
@@ -641,12 +762,12 @@ struct dNdetaDrawer
     p2->cd();
 
     // Fix up axis
-    FixAxis(ratios, 1/yd/1.7, "Ratios", 7);
+    FixAxis(fRatios, 1/yd/1.7, "Ratios", 7);
 
-    ratios->SetMaximum(1+TMath::Max(.22,1.05*max));
-    ratios->SetMinimum(1-TMath::Max(.32,1.05*max));
+    fRatios->SetMaximum(1+TMath::Max(.22,1.05*max));
+    fRatios->SetMinimum(1-TMath::Max(.32,1.05*max));
     p2->Clear();
-    ratios->DrawClone("nostack e1");
+    fRatios->DrawClone("nostack e1");
 
     
     // Make a legend
@@ -672,15 +793,15 @@ struct dNdetaDrawer
     band->DrawClone("X L same");
     
     // Replot the ratios on top
-    ratios->DrawClone("nostack e1 same");
+    fRatios->DrawClone("nostack e1 same");
 
     if (isBottom) {
-      fRangeParam->fMasterAxis = FindXAxis(p2, ratios->GetName());
+      fRangeParam->fMasterAxis = FindXAxis(p2, fRatios->GetName());
       p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)", 
                                fRangeParam));
     }
     else { 
-      fRangeParam->fSlave2Axis = FindXAxis(p2, ratios->GetName());
+      fRangeParam->fSlave2Axis = FindXAxis(p2, fRatios->GetName());
       fRangeParam->fSlave2Pad  = p2;
     }
   }
@@ -688,15 +809,14 @@ struct dNdetaDrawer
   /** 
    * Plot the asymmetries
    * 
-   * @param ratios  Asymmetries to plot (if any)
    * @param max     Maximum diviation from 1 
    * @param y1      Lower y coordinate of pad
    * @param y2      Upper y coordinate of pad
    */
-  void PlotLeftRight(THStack* leftright, Double_t max, 
-                    Double_t y1, Double_t y2) 
+  void PlotLeftRight(Double_t max, Double_t y1, Double_t y2) 
   {
-    if (!leftright) return;
+    if (!fShowLeftRight) return;
+
     bool isBottom = (y1 < 0.0001);
     Double_t yd = y2 - y1;
     // Make a sub-pad for the result itself
@@ -711,22 +831,22 @@ struct dNdetaDrawer
     p3->cd();
 
     TH1* dummy = 0;
-    if (leftright->GetHists()->GetEntries() == 1) { 
+    if (fLeftRight->GetHists()->GetEntries() == 1) { 
       // Add dummy histogram
       dummy = new TH1F("dummy","", 10, -6, 6);
       dummy->SetLineColor(0);
       dummy->SetFillColor(0);
       dummy->SetMarkerColor(0);
-      leftright->Add(dummy);
+      fLeftRight->Add(dummy);
     }
 
     // Fix up axis
-    FixAxis(leftright, 1/yd/1.7, "Right/Left", 4);
+    FixAxis(fLeftRight, 1/yd/1.7, "Right/Left", 4);
 
-    leftright->SetMaximum(1+TMath::Max(.12,1.05*max));
-    leftright->SetMinimum(1-TMath::Max(.15,1.05*max));
+    fLeftRight->SetMaximum(1+TMath::Max(.12,1.05*max));
+    fLeftRight->SetMinimum(1-TMath::Max(.15,1.05*max));
     p3->Clear();
-    leftright->DrawClone("nostack e1");
+    fLeftRight->DrawClone("nostack e1");
 
     
     // Make a legend
@@ -762,14 +882,14 @@ struct dNdetaDrawer
     band->Draw("3 same");
     band->DrawClone("X L same");
 
-    leftright->DrawClone("nostack e1 same");
+    fLeftRight->DrawClone("nostack e1 same");
     if (isBottom) {
-      fRangeParam->fMasterAxis = FindXAxis(p3, leftright->GetName());
+      fRangeParam->fMasterAxis = FindXAxis(p3, fLeftRight->GetName());
       p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)", 
                                fRangeParam));
     }
     else { 
-      fRangeParam->fSlave2Axis = FindXAxis(p3, leftright->GetName());
+      fRangeParam->fSlave2Axis = FindXAxis(p3, fLeftRight->GetName());
       fRangeParam->fSlave2Pad  = p3;
     }
   }
@@ -791,14 +911,7 @@ struct dNdetaDrawer
   {
     if (!list) return 0;
     
-    TList* all = static_cast<TList*>(list->FindObject("all"));
-    if (!all) { 
-      Warning("GetResult", "No 'all' list find in %s", list->GetName());
-      // list->ls();
-      return 0;
-    }
-
-    TH1* ret = static_cast<TH1*>(all->FindObject(name));
+    TH1* ret = static_cast<TH1*>(list->FindObject(name));
     if (!ret) {
       // all->ls();
       Warning("GetResult", "Histogram %s not found", name);
@@ -816,7 +929,7 @@ struct dNdetaDrawer
    * 
    * @return Maximum of histogram 
    */
-  Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option) const 
+  Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option="") const 
   {
     // Check if we have input 
     if (!hist) return 0;
@@ -824,7 +937,10 @@ struct dNdetaDrawer
     // Rebin if needed 
     Rebin(hist);
 
+    // Info("AddHistogram", "Adding %s to %s", 
+    //      hist->GetName(), stack->GetName());
     stack->Add(hist, option);
+    // stack->ls();
     return hist->GetMaximum();
   }
   //__________________________________________________________________
@@ -838,8 +954,8 @@ struct dNdetaDrawer
    * 
    * @return Maximum of histogram 
    */
-  Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option, 
-                       TH1*& sym) const 
+  Double_t AddHistogram(THStack* stack, TH1* hist, TH1*& sym,
+                       Option_t* option="") const 
   {
     // Check if we have input 
     if (!hist) return 0;
@@ -852,6 +968,8 @@ struct dNdetaDrawer
     sym = Symmetrice(hist);
     stack->Add(sym, option);
 
+    // Info("AddHistogram", "Adding %s and %s to %s", 
+    //      hist->GetName(), sym->GetName(), stack->GetName());
     return hist->GetMaximum();
   }
 
@@ -990,6 +1108,7 @@ struct dNdetaDrawer
     // Double_t dBin  = (high - low) / oBins;
     // Int_t    tBins = Int_t(2*high/dBin+.5);
     // ret->SetBins(tBins, -high, high);
+    ret->SetDirectory(0);
     ret->Reset();
     ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
     ret->SetYTitle("Right/Left");
@@ -1134,6 +1253,7 @@ struct dNdetaDrawer
     ret->SetMarkerStyle(g->GetMarkerStyle());
     ret->SetMarkerColor(h->GetMarkerColor());
     ret->SetMarkerSize(0.9*g->GetMarkerSize());
+    ret->SetDirectory(0);
     Double_t xlow  = g->GetX()[0];
     Double_t xhigh = g->GetX()[g->GetN()-1];
     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
@@ -1178,6 +1298,7 @@ struct dNdetaDrawer
     t1->Divide(h2);
     t1->SetMarkerColor(h1->GetMarkerColor());
     t1->SetMarkerStyle(h2->GetMarkerStyle());
+    t1->SetDirectory(0);
     max = TMath::Max(RatioMax(t1), max);
     return t1;
   }
@@ -1218,6 +1339,7 @@ struct dNdetaDrawer
     h->SetMarkerStyle(g2->GetMarkerStyle());
     h->SetMarkerColor(g1->GetMarkerColor());
     h->SetMarkerSize(0.9*g2->GetMarkerSize());
+    h->SetDirectory(0);
 
     Double_t low  = g2->GetX()[0];
     Double_t high = g2->GetX()[g2->GetN()-1];
@@ -1282,11 +1404,18 @@ struct dNdetaDrawer
   void FixAxis(THStack* stack, Double_t s, const char* ytitle,
                Int_t ynDiv=210, Bool_t force=true)
   {
+    if (!stack) { 
+      Warning("FixAxis", "No stack passed for %s!", ytitle);
+      return;
+    }
     if (force) stack->Draw("nostack e1");
 
     TH1* h = stack->GetHistogram();
-    if (!h) return;
-
+    if (!h) { 
+      Warning("FixAxis", "Stack %s has no histogram", stack->GetName());
+      return;
+    }
+    
     h->SetXTitle("#eta");
     h->SetYTitle(ytitle);
     TAxis* xa = h->GetXaxis();
@@ -1297,6 +1426,11 @@ struct dNdetaDrawer
       xa->SetTitleSize(s*xa->GetTitleSize());
       xa->SetLabelSize(s*xa->GetLabelSize());
       xa->SetTickLength(s*xa->GetTickLength());
+
+      if (stack != fResults) {
+       TAxis* rxa = fResults->GetXaxis();
+       xa->Set(rxa->GetNbins(), rxa->GetXmin(), rxa->GetXmax());
+      }
     }
     if (ya) {
       ya->SetTitle(ytitle);
@@ -1313,31 +1447,24 @@ struct dNdetaDrawer
 
 
   //__________________________________________________________________
-  Bool_t      fShowOthers;   // Show other data
-  Bool_t      fShowAlice;    // Show ALICE published data
-  Bool_t      fShowRatios;   // Show ratios 
-  Bool_t      fShowLeftRight;// Show asymmetry 
-  UShort_t    fRebin;        // Rebinning factor 
-  Bool_t      fCutEdges;     // Whether to cut edges
-  TString     fTitle;        // Title on plot
-  TNamed*     fTrigString;   // Trigger string (read, or set)
-  TNamed*     fSNNString;    // Energy string (read, or set)
-  TNamed*     fSysString;    // Collision system string (read or set)
-  TAxis*      fVtxAxis;      // Vertex cuts (read or set)
-  TList       fNumerators;   // Numerators in ratios 
-  TList       fDenominators; // Denominators in ratios 
-  THStack*    fResults;      // Stack of results 
-  THStack*    fRatios;       // Stack of ratios 
-  // TH1*        fForward;      // Results
-  // TH1*        fForwardMC;    // MC results
-  // TH1*        fForwardHHD;   // Old results
-  // TH1*        fTruth;        // MC truth
-  // TH1*        fCentral;      // Central region data
-  // TH1*        fForwardSym;   // Symmetric extension
-  // TH1*        fForwardMCSym; // Symmetric extension
-  // TH1*        fForwardHHDSym;// Symmetric extension
-  TH1*        fTriggers;     // Number of triggers
-  RangeParam* fRangeParam;   // Parameter object for range zoom 
+  Bool_t       fShowOthers;   // Show other data
+  Bool_t       fShowAlice;    // Show ALICE published data
+  Bool_t       fShowRatios;   // Show ratios 
+  Bool_t       fShowLeftRight;// Show asymmetry 
+  UShort_t     fRebin;        // Rebinning factor 
+  Bool_t       fCutEdges;     // Whether to cut edges
+  TString      fTitle;        // Title on plot
+  TNamed*      fTrigString;   // Trigger string (read, or set)
+  TNamed*      fSNNString;    // Energy string (read, or set)
+  TNamed*      fSysString;    // Collision system string (read or set)
+  TAxis*       fVtxAxis;      // Vertex cuts (read or set)
+  TAxis*       fCentAxis;     // Centrality axis
+  THStack*     fResults;      // Stack of results 
+  THStack*     fRatios;       // Stack of ratios 
+  THStack*     fLeftRight;    // Left-right asymmetry
+  TMultiGraph* fOthers;       // Older data 
+  TH1*         fTriggers;     // Number of triggers
+  RangeParam*  fRangeParam;   // Parameter object for range zoom 
   
 };
 
index f251c08..fa585af 100644 (file)
@@ -95,6 +95,7 @@ void MakeAOD(const char* esddir,
   // --- AOD output handler ------------------------------------------
   AliAODHandler* aodHandler   = new AliAODHandler();
   mgr->SetOutputEventHandler(aodHandler);
+  aodHandler->SetNeedsHeaderReplication();
   aodHandler->SetOutputFileName("AliAOD.root");
 
   // --- Add tasks ---------------------------------------------------
@@ -133,7 +134,7 @@ void MakeAOD(const char* esddir,
   mgr->SetSkipTerminate(false);
   // Some informative output 
   mgr->PrintStatus();
-  // mgr->SetDebugLevel(3);
+  if (proof) mgr->SetDebugLevel(3);
   if (mgr->GetDebugLevel() < 1 && !proof) 
     mgr->SetUseProgressBar(kTRUE,100);
 
index 9b7536d..bd10808 100644 (file)
  *
  * @ingroup pwg2_forward_scripts
  */
-void MakedNdeta(const char* aoddir=".", 
-               Int_t       nEvents=-1, 
-               const char* trig="INEL",
-               Float_t     centLow = 0,
-               Float_t     centHigh = 100,
-               Double_t    vzMin=-10,
-               Double_t    vzMax=10,
-               Int_t       proof=0)
+void MakedNdeta(const char* aoddir   = ".", 
+               Int_t       nEvents  = -1, 
+               const char* trig     = "INEL",
+               Bool_t      useCent  = false,
+               Double_t    vzMin    = -10,
+               Double_t    vzMax    = +10,
+               Int_t       proof    = 0)
 {
   // --- Libraries to load -------------------------------------------
   gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
@@ -61,10 +60,10 @@ void MakedNdeta(const char* aoddir=".",
   // --- Add tasks ---------------------------------------------------
   // Forward 
   gROOT->LoadMacro("AddTaskForwarddNdeta.C");
-  AddTaskForwarddNdeta(trig, vzMin, vzMax, centLow, centHigh, true);
+  AddTaskForwarddNdeta(trig, vzMin, vzMax, useCent, true);
   // Central
   gROOT->LoadMacro("AddTaskCentraldNdeta.C");
-  AddTaskCentraldNdeta(trig, vzMin, vzMax, centLow, centHigh);
+  AddTaskCentraldNdeta(trig, vzMin, vzMax, useCent);
 
   
   // --- Run the analysis --------------------------------------------
index e0a4863..f7fa944 100644 (file)
@@ -585,6 +585,7 @@ TGraphAsymmErrors* CMSNsd7000()
 /** 
  * Get a multi graph of data for a given energy and trigger type 
  * 
+ * @param sys    Collision system (1: pp, 2: PbPb)
  * @param energy Energy in GeV (900, 2360, 7000)
  * @param type   Bit pattern of trigger type 
  *   - 0x1 INEL 
@@ -596,63 +597,94 @@ TGraphAsymmErrors* CMSNsd7000()
  * @ingroup pwg2_forward_otherdata
  */
 TMultiGraph* 
-GetData(Int_t energy, Int_t type, bool aliceOnly=false)
+GetData(UShort_t sys, 
+       UShort_t energy,
+       UShort_t type=0x1, 
+       UShort_t centLow=0, 
+       UShort_t centHigh=0, 
+       bool     aliceOnly=false)
 {
-  TMultiGraph* mp = new TMultiGraph(Form("dndeta_%dGeV_%d", energy, type),"");
+  TMultiGraph* mp = new TMultiGraph(Form("dndeta_%dGeV_%d_%03d_%03d", 
+                                        energy, type, centLow, centHigh),"");
   TString tn;
   TString en;
-  if (TMath::Abs(energy-900) < 10) {
-    en.Append(" #sqrt{s}=900GeV");
-    if (type & 0x1) { 
-      tn.Append(" INEL");
-      if (!aliceOnly) mp->Add(UA5Inel(false));
-      if (!aliceOnly) mp->Add(UA5Inel(true));
-      mp->Add(AliceCentralInel900());
-    }      
-    if (type & 0x4) { 
-      tn.Append(" NSD");
-      if (!aliceOnly) mp->Add(UA5Nsd(false));
-      if (!aliceOnly) mp->Add(UA5Nsd(true));
-      mp->Add(AliceCentralNsd900());
-      if (!aliceOnly) mp->Add(CMSNsd900());
+  TString sn;
+  TString cn;
+  if (sys == 1) { 
+    sn = ", pp(p#bar{p})";
+    if (energy < 1000) 
+      en = Form(", #sqrt{s}=%dGeV", energy);
+    else 
+      en = Form(", #sqrt{s}=%f4.2TeV", float(energy)/1000);
+    if (!(type & 0x7)) 
+      Warning("GetData", "Unknown trigger mask 0x%x", type);
+
+    if (TMath::Abs(energy-900) < 10) {
+      if (type & 0x1) { 
+       tn.Append(" INEL");
+       if (!aliceOnly) mp->Add(UA5Inel(false));
+       if (!aliceOnly) mp->Add(UA5Inel(true));
+       mp->Add(AliceCentralInel900());
+      }      
+      if (type & 0x4) { 
+       tn.Append(" NSD");
+       if (!aliceOnly) mp->Add(UA5Nsd(false));
+       if (!aliceOnly) mp->Add(UA5Nsd(true));
+       mp->Add(AliceCentralNsd900());
+       if (!aliceOnly) mp->Add(CMSNsd900());
+      }
+      if (type & 0x2) { 
+       tn.Append(" INEL>0");
+       mp->Add(AliceCentralInelGt900());
+      }
     }
-    if (type & 0x2) { 
-      tn.Append(" INEL>0");
-      mp->Add(AliceCentralInelGt900());
+    else if (TMath::Abs(energy-2360) < 10) {
+      if (type & 0x1) { 
+       tn.Append(" INEL");
+       mp->Add(AliceCentralInel2360());
+      }
+      if (type & 0x4) { 
+       tn.Append(" NSD");
+       mp->Add(AliceCentralNsd2360());
+       if (!aliceOnly) mp->Add(CMSNsd2360());
+      }
+      if (type & 0x2) { 
+       tn.Append(" INEL>0");
+       mp->Add(AliceCentralInelGt2360());
+      }
     }
-  }
-  if (TMath::Abs(energy-2360) < 10) {
-    en.Append(" #sqrt{s}=2.36TeV");
-    if (type & 0x1) { 
-      tn.Append(" INEL");
-      mp->Add(AliceCentralInel2360());
-    }
-    if (type & 0x4) { 
-      tn.Append(" NSD");
-      mp->Add(AliceCentralNsd2360());
-      if (!aliceOnly) mp->Add(CMSNsd2360());
-    }
-    if (type & 0x1) { 
-      tn.Append(" INEL>0");
-      mp->Add(AliceCentralInelGt2360());
+    else if (TMath::Abs(energy-7000) < 10) {
+      if (type & 0x1) { 
+       tn.Append(" INEL");
+      }
+      if (type & 0x4) { 
+       tn.Append(" NSD");
+       if (!aliceOnly) mp->Add(CMSNsd7000());
+      }
+      if (type & 0x2) { 
+       tn.Append(" INEL>0");
+       mp->Add(AliceCentralInelGt7000());
+      }
     }
+    else 
+      Warning("GetData", "No other results for sys=%d, energy=%d",
+             sys, energy);
   }
-  if (TMath::Abs(energy-7000) < 10) {
-    en.Append(" #sqrt{s}=7TeV");
-    if (type & 0x1) { 
-      tn.Append(" INEL");
-    }
-    if (type & 0x4) { 
-      tn.Append(" NSD");
-      if (!aliceOnly) mp->Add(CMSNsd7000());
-    }
-    if (type & 0x1) { 
-      tn.Append(" INEL>0");
-      mp->Add(AliceCentralInelGt7000());
-    }
+  else if (sys == 2) { 
+    // Nothing for PbPb so far 
+    cn = Form(", %d%%-%d%% central", centLow, centHigh);
+    sn = ", PbPb";
+    if (energy < 1000) 
+      en = Form(", #sqrt{s_{NN}}=%dGeV", energy);
+    else 
+      en = Form(", #sqrt{s_{NN}}=%f4.2TeV", float(energy)/1000);
+    Warning("GetData", "No other data for PbP b yet");
   }
-  mp->SetTitle(Form("1/N dN_{ch}/d#eta, pp(p#bar{p}), %s, %s", 
-                   en.Data(), tn.Data()));
+  else 
+    Warning("GetData", "Unknown system %d", sys);
+  TString tit(Form("1/N dN_{ch}/d#eta%s%s%s%s", 
+                  sn.Data(), en.Data(), tn.Data(), cn.Data()));
+  mp->SetTitle(tit.Data());
   if (!mp->GetListOfGraphs() || mp->GetListOfGraphs()->GetEntries() <= 0) {
     delete mp;
     mp = 0;
@@ -671,9 +703,14 @@ GetData(Int_t energy, Int_t type, bool aliceOnly=false)
  * @ingroup pwg2_forward_otherdata
  */
 void
-OtherData(Int_t energy=900, Int_t type=0x1, bool aliceOnly=false)
+OtherData(UShort_t sys=1, 
+             UShort_t energy=900, 
+             UShort_t type=0x1, 
+             UShort_t centLow=0, 
+             UShort_t centHigh=5, 
+             bool     aliceOnly=false)
 {
-  TMultiGraph* mp = GetData(energy, type, aliceOnly);
+  TMultiGraph* mp = GetData(sys, energy, type, centLow, centHigh, aliceOnly);
   if (!mp) return;
 
   gStyle->SetTitleX(0.1);
index 97ca8e8..68c5abd 100755 (executable)
@@ -5,8 +5,6 @@ nev=-1
 rebin=1
 vzmin=-10
 vzmax=10
-centlow=0
-centhigh=100
 batch=0
 gdb=0
 proof=0
@@ -177,7 +175,7 @@ fi
 if test $dopass2 -gt 0 ; then
     rotate ${output2}
 
-    args=(\(\".\",$nev,\"$type\",$centlow,$centhigh,$vzmin,$vzmax,$proof\))
+    args=(\(\".\",$nev,\"$type\",$cent,$vzmin,$vzmax,$proof\))
     if test "x$pass1" = "xMakeELossFits.C" ; then 
        args=(\(\"${output1}\"\))
     fi