Code update for handling centrality in the dN/deta analysis.
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 15 Mar 2011 15:10:21 +0000 (15:10 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 15 Mar 2011 15:10:21 +0000 (15:10 +0000)
New MC specific event inspector
Various updates to run script
More flags for the histogram collector
and a few other things

21 files changed:
PWG2/FORWARD/analysis2/AddTaskCentraldNdeta.C
PWG2/FORWARD/analysis2/AddTaskFMD.C
PWG2/FORWARD/analysis2/AddTaskForwarddNdeta.C
PWG2/FORWARD/analysis2/AliBasedNdetaTask.cxx
PWG2/FORWARD/analysis2/AliBasedNdetaTask.h
PWG2/FORWARD/analysis2/AliFMDEventInspector.cxx
PWG2/FORWARD/analysis2/AliFMDEventInspector.h
PWG2/FORWARD/analysis2/AliFMDHistCollector.cxx
PWG2/FORWARD/analysis2/AliFMDHistCollector.h
PWG2/FORWARD/analysis2/AliFMDMCEventInspector.cxx [new file with mode: 0644]
PWG2/FORWARD/analysis2/AliFMDMCEventInspector.h [new file with mode: 0644]
PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx
PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.h
PWG2/FORWARD/analysis2/AliForwarddNdetaTask.cxx
PWG2/FORWARD/analysis2/AliForwarddNdetaTask.h
PWG2/FORWARD/analysis2/DrawdNdeta.C
PWG2/FORWARD/analysis2/ForwardAODConfig.C
PWG2/FORWARD/analysis2/MakeAOD.C
PWG2/FORWARD/analysis2/MakeELossFits.C
PWG2/FORWARD/analysis2/MakedNdeta.C
PWG2/FORWARD/analysis2/Run.sh

index 4b87813..116ebfa 100644 (file)
@@ -9,7 +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, 
+                    Float_t centLow=0, 
+                    Float_t centHigh=100)
 {
   // analysis manager
   AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
@@ -19,18 +23,16 @@ AddTaskCentraldNdeta(const char* trig="INEL", Double_t vzMin=-10, Double_t vzMax
   AliCentraldNdetaTask* task = new AliCentraldNdetaTask("Central");
   task->SetVertexRange(vzMin, vzMax);
   task->SetTriggerMask(trig);
-  task->SetCentLow(centLow);
-  task->SetCentHigh(centHigh);
-  
-  /*
-  if(trig == "INEL") task->SetTriggerEff(0.95);
-  if(trig == "NSD")  task->SetTriggerEff(1.04);
-  TFile f("/home/canute/ALICE/FMDanalysis/productionData/normalizationHists900GeV.root", "READ");
-  //TFile f("/home/canute/ALICE/FMDanalysis/BackgroundCorrection/normalizationHists.root", "READ");
-  TH2D* hnorm = 0 ;
-  if(trig == "INEL") hnorm = (TH2D*)f.Get("hInelNormalization");
-  if(trig == "NSD") hnorm = (TH2D*)f.Get("hNSDNormalization");
-  task->SetShapeCorrection(hnorm);*/
+  // 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);
   mgr->AddTask(task);
 
   // create containers for input/output
index 08cc413..6ce04bb 100644 (file)
@@ -35,8 +35,6 @@ AddTaskFMD(Bool_t mc, UShort_t sys=0, UShort_t sNN=0, Short_t field=0)
     AliForwardCorrectionManager::Instance().Init(sys,sNN,field);
 
   // --- Configure the task ------------------------------------------
-  gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2", 
-                          gROOT->GetMacroPath()));
   const char* config = gSystem->Which(gROOT->GetMacroPath(),
                                      "ForwardAODConfig.C");
   if (!config) 
@@ -44,8 +42,8 @@ AddTaskFMD(Bool_t mc, UShort_t sys=0, UShort_t sNN=0, Short_t field=0)
            gROOT->GetMacroPath());
   else {
     Info("AddTaskFMD", 
-        "Loading configuration of AliForwardMultiplicityTask from %s",
-        config);
+        "Loading configuration of '%s' from %s",
+        task->ClassName(), config);
     gROOT->Macro(Form("%s((AliForwardMultiplicityBase*)%p)", config, task));
     delete config;
   }
index 27c0d58..3a41490 100644 (file)
@@ -9,7 +9,12 @@
  * 
  */
 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, 
+                    Float_t centlow = 0, 
+                    Float_t centhigh = 100, 
+                    Bool_t cutEdges = false)
 {
   // analysis manager
   AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
@@ -20,18 +25,16 @@ AddTaskForwarddNdeta(const char* trig="INEL", Double_t vzMin=-10, Double_t vzMax
   task->SetVertexRange(vzMin, vzMax);
   task->SetTriggerMask(trig);
   task->SetCutEdges(cutEdges);
-  task->SetCentLow(centlow);
-  task->SetCentHigh(centhigh);
-   
-  //task->SetTriggerEff(0.93);
-  /*if(trig == "INEL") task->SetTriggerEff(0.95);
-  if(trig == "NSD")  task->SetTriggerEff(1.04);
-  TFile f("/home/canute/ALICE/FMDanalysis/productionData/normalizationHists900GeV.root", "READ");
-  //TFile f("/home/canute/ALICE/FMDanalysis/BackgroundCorrection/normalizationHists.root", "READ");
-  TH2D* hnorm = 0 ;
-  if(trig == "INEL") hnorm = (TH2D*)f.Get("hInelNormalization");
-  if(trig == "NSD") hnorm = (TH2D*)f.Get("hNSDNormalization");
-  task->SetShapeCorrection(hnorm);*/
+  // 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);
   mgr->AddTask(task);
 
   // create containers for input/output
index b382518..3913591 100644 (file)
 //____________________________________________________________________
 AliBasedNdetaTask::AliBasedNdetaTask()
   : AliAnalysisTaskSE(), 
-    fSum(0),           //  Sum of histograms 
-    fSumMC(0),         //  Sum of MC histograms (if any)
     fSums(0),          // Container of sums 
-    fOutput(0),                // Container of outputs 
-    fTriggers(0),      // Histogram of triggers 
+    fOutput(0),         // Container of output 
     fVtxMin(0),                // Minimum v_z
     fVtxMax(0),                // Maximum v_z
     fTriggerMask(0),    // Trigger mask 
@@ -31,8 +28,11 @@ AliBasedNdetaTask::AliBasedNdetaTask()
     fCorrEmpty(true), 
     fTriggerEff(1),
     fShapeCorr(0),
-    fCentLow(0),
-    fCentHigh(100)
+    fListOfCentralities(0),
+    fUseShapeCorr(true),
+    fSNNString(0),
+    fSysString(0),
+    fCent(0)
 {
   // 
   // Constructor
@@ -42,11 +42,8 @@ AliBasedNdetaTask::AliBasedNdetaTask()
 //____________________________________________________________________
 AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
   : AliAnalysisTaskSE(name), 
-    fSum(0),           // Sum of histograms 
-    fSumMC(0),         // Sum of MC histograms (if any)
     fSums(0),          // Container of sums 
-    fOutput(0),                // Container of outputs 
-    fTriggers(0),      // Histogram of triggers 
+    fOutput(0),         // Container of output 
     fVtxMin(-10),      // Minimum v_z
     fVtxMax(10),       // Maximum v_z
     fTriggerMask(AliAODForwardMult::kInel), 
@@ -56,12 +53,16 @@ AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
     fCorrEmpty(true), 
     fTriggerEff(1),
     fShapeCorr(0),
-    fCentLow(0),
-    fCentHigh(100)
+    fListOfCentralities(0),
+    fUseShapeCorr(true),
+    fSNNString(0),
+    fSysString(0),
+    fCent(0)
 {
   // 
   // Constructor
   // 
+  fListOfCentralities = new TList;
 
   // Output slot #1 writes into a TH1 container
   DefineOutput(1, TList::Class()); 
@@ -71,11 +72,8 @@ AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
 //____________________________________________________________________
 AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
   : AliAnalysisTaskSE(o), 
-    fSum(o.fSum),      // TH2D* -  Sum of histograms 
-    fSumMC(o.fSumMC),// TH2D* -  Sum of MC histograms (if any)
     fSums(o.fSums),            // TList* - Container of sums 
-    fOutput(o.fOutput),                // TList* - Container of outputs 
-    fTriggers(o.fTriggers),    // TH1D* - Histogram of triggers 
+    fOutput(o.fOutput),         // Container of output 
     fVtxMin(o.fVtxMin),                // Double_t - Minimum v_z
     fVtxMax(o.fVtxMax),                // Double_t - Maximum v_z
     fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask 
@@ -85,8 +83,11 @@ AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
     fCorrEmpty(true), 
     fTriggerEff(o.fTriggerEff),
     fShapeCorr(o.fShapeCorr),
-    fCentLow(o.fCentLow),
-    fCentHigh(o.fCentHigh)
+    fListOfCentralities(o.fListOfCentralities),
+    fUseShapeCorr(o.fUseShapeCorr),
+    fSNNString(o.fSNNString),
+    fSysString(o.fSysString),
+    fCent(o.fCent)
 {}
 
 //____________________________________________________________________
@@ -100,15 +101,47 @@ AliBasedNdetaTask::~AliBasedNdetaTask()
     delete fSums;
     fSums = 0;
   }
-  if (fOutputs) { 
-    fOutputs->Delete();
-    delete fOutputs;
-    fOutputs = 0;
+  if (fOutput) { 
+    fOutput->Delete();
+    delete fOutput;
+    fOutput = 0;
   }
 }
 
 //________________________________________________________________________
 void 
+AliBasedNdetaTask::AddCentralityBin(Short_t low, Short_t high)
+{
+  // 
+  // Add a centrality bin 
+  // 
+  // Parameters:
+  //    low  Low cut
+  //    high High cut
+  //
+  fListOfCentralities->Add(MakeCentralityBin(GetName(), low, high));
+}
+
+//________________________________________________________________________
+AliBasedNdetaTask::CentralityBin*
+AliBasedNdetaTask::MakeCentralityBin(const char* name, 
+                                    Short_t low, Short_t high) const
+{
+  // 
+  // Make a centrality bin 
+  // 
+  // Parameters:
+  //    name  Name used for histograms
+  //    low   Low cut in percent
+  //    high  High cut in percent
+  // 
+  // Return:
+  //    A newly created centrality bin 
+  //
+  return new CentralityBin(name, low, high);
+}
+//________________________________________________________________________
+void 
 AliBasedNdetaTask::SetTriggerMask(const char* mask)
 {
   // 
@@ -160,33 +193,24 @@ AliBasedNdetaTask::UserCreateOutputObjects()
   //
   // This is called once per slave process 
   //
-  fOutput = new TList;
-  fOutput->SetName(Form("%s_result", GetName()));
-  fOutput->SetOwner();
-
   fSums = new TList;
   fSums->SetName(Form("%s_sums", GetName()));
   fSums->SetOwner();
 
+  // Automatically add 'all' centrality bin if nothing has been defined. 
+  if (fListOfCentralities->GetEntries() <= 0) AddCentralityBin(0,0); 
 
-  fTriggers = new TH1D("triggers", "Number of triggers", 
-                      kAccepted, 1, kAccepted);
-  fTriggers->SetYTitle("# of events");
-  fTriggers->GetXaxis()->SetBinLabel(kAll,         "All events");
-  fTriggers->GetXaxis()->SetBinLabel(kB,           "w/B trigger");
-  fTriggers->GetXaxis()->SetBinLabel(kA,           "w/A trigger");
-  fTriggers->GetXaxis()->SetBinLabel(kC,           "w/C trigger");
-  fTriggers->GetXaxis()->SetBinLabel(kE,           "w/E trigger");
-  fTriggers->GetXaxis()->SetBinLabel(kMB,          "w/Collision trigger");
-  fTriggers->GetXaxis()->SetBinLabel(kWithVertex,  "w/Vertex");
-  fTriggers->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger");
-  fTriggers->GetXaxis()->SetBinLabel(kPileUp,      "with pileup");
-  fTriggers->GetXaxis()->SetBinLabel(kAccepted,    "Accepted by cut");
-  fTriggers->GetXaxis()->SetNdivisions(kAccepted, false);
-  fTriggers->SetFillColor(kRed+1);
-  fTriggers->SetFillStyle(3001);
-  fTriggers->SetStats(0);
-  fSums->Add(fTriggers);
+  // Centrality histogram 
+  fCent = new TH1D("cent", "Centrality", 100, 0, 100);
+  fCent->SetDirectory(0);
+  fCent->SetXTitle(0);
+  fSums->Add(fCent);
+
+  // Loop over centrality bins 
+  TIter next(fListOfCentralities);
+  CentralityBin* bin = 0;
+  while ((bin = static_cast<CentralityBin*>(next()))) 
+    bin->CreateOutputObjects(fSums);
 
   // Check that we have an AOD input handler 
   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
@@ -197,87 +221,6 @@ AliBasedNdetaTask::UserCreateOutputObjects()
   // Post data for ALL output slots >0 here, to get at least an empty histogram
   PostData(1, fSums); 
 }
-
-//____________________________________________________________________
-TH2D*
-AliBasedNdetaTask::CloneHist(const TH2D* in, const char* name) 
-{
-  // 
-  // Clone a 2D histogram
-  // 
-  // Parameters:
-  //    in    Histogram to clone.
-  //    name  New name of clone.
-  // 
-  // Return:
-  //    The clone
-  //
-  if (!in) return 0;
-  TH2D* ret = static_cast<TH2D*>(in->Clone(name));
-  ret->SetDirectory(0);
-  ret->Sumw2();
-  fSums->Add(ret);
-
-  return ret;
-}
-
-//____________________________________________________________________
-Bool_t
-AliBasedNdetaTask::CheckEvent(const AliAODEvent* aod)
-{
-  // 
-  // Check the trigger and vertex 
-  // 
-  // Parameters:
-  //    aod 
-  // 
-  // Return:
-  //    
-  //
-  TObject*           oForward   = aod->FindListObject("Forward");
-  if (!oForward) { 
-    AliWarning("No forward object found");
-    return false;
-  }
-  AliAODForwardMult* forward   = static_cast<AliAODForwardMult*>(oForward);
-  
-  // Count event 
-  fTriggers->AddBinContent(kAll);
-  if (forward->IsTriggerBits(AliAODForwardMult::kB)) 
-    fTriggers->AddBinContent(kB);
-  if (forward->IsTriggerBits(AliAODForwardMult::kA)) 
-    fTriggers->AddBinContent(kA);
-  if (forward->IsTriggerBits(AliAODForwardMult::kC)) 
-    fTriggers->AddBinContent(kC);
-  if (forward->IsTriggerBits(AliAODForwardMult::kE)) 
-    fTriggers->AddBinContent(kE);
-  if (forward->IsTriggerBits(AliAODForwardMult::kInel)) 
-    fTriggers->AddBinContent(kMB);
-  if (forward->IsTriggerBits(AliAODForwardMult::kMCNSD)) 
-    fTriggers->AddBinContent(kMCNSD);
-  if (forward->IsTriggerBits(AliAODForwardMult::kPileUp)) 
-    fTriggers->AddBinContent(kPileUp);
-  
-  
-  
-  // Check if we have an event of interest. 
-  if (!forward->IsTriggerBits(fTriggerMask)) return false;
-  
-  //Check for pileup
-  if (forward->IsTriggerBits(AliAODForwardMult::kPileUp)) return false;
-  
-  fTriggers->AddBinContent(kWithTrigger);
-  
-  // Check that we have a valid vertex
-  if (!forward->HasIpZ()) return false;
-  fTriggers->AddBinContent(kWithVertex);
-
-  // Check that vertex is within cuts 
-  if (!forward->InRange(fVtxMin, fVtxMax)) return false;
-  fTriggers->AddBinContent(kAccepted);
-
-  return true;
-}
 //____________________________________________________________________
 void 
 AliBasedNdetaTask::UserExec(Option_t *) 
@@ -295,47 +238,40 @@ AliBasedNdetaTask::UserExec(Option_t *)
     return;
   }  
   
-  TObject* obj = 0;
-  obj = aod->FindListObject("Forward");
-  
-  AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
-  Float_t cent = forward->GetCentrality();
-  
-  if(cent > 0) {
-    if( cent < fCentLow || cent > fCentHigh ) return;
-    //else std::cout<<"selecting event with cent "<<cent<<std::endl;
+  TObject* obj = aod->FindListObject("Forward");
+  if (!obj) { 
+    AliWarning("No forward object found");
+    return;
   }
-  
+  AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
+
+  // Fill centrality histogram 
+  fCent->Fill(forward->GetCentrality());
+
   // Get the histogram(s) 
   TH2D* data   = GetHistogram(aod, false);
   TH2D* dataMC = GetHistogram(aod, true);
 
-  // We should have a base object at least 
-  if (!data) { 
-    AliWarning("No data object found in AOD");
-    return;
-  }
-  
-  
+  // Loop over centrality bins 
+  TIter next(fListOfCentralities);
+  CentralityBin* bin = 0;
+  while ((bin = static_cast<CentralityBin*>(next()))) 
+    bin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
 
-  // Check the event 
-  if (!CheckEvent(aod)) return;
+  // Here, we get the update 
+  if (!fSNNString) { 
+    UShort_t sNN = forward->GetSNN();
+    fSNNString = new TNamed("sNN", "");
+    fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
+    fSNNString->SetUniqueID(sNN);
+    fSums->Add(fSNNString);
   
-  // Create our sum histograms 
-  if (!fSum)  { 
-    fSum   = CloneHist(data,  GetName()); 
-    if(forward->GetSystem()  < 2)
-      LoadNormalizationData(forward->GetSystem(),forward->GetSNN());
+    UShort_t sys = forward->GetSystem();
+    fSysString = new TNamed("sys", "");
+    fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
+    fSysString->SetUniqueID(sys);
+    fSums->Add(fSysString);
   }
-  else fSum->Add((data));
-  
-  //for(Int_t i = 1; i<=data->GetNbinsX(); i++) {
-  //  std::cout<<fSum->GetBinContent(i,0) <<"   "<<data->GetBinContent(i,0)<<std::endl;
-  // }
-  
-  if (!fSumMC && dataMC) fSumMC = CloneHist(dataMC,Form("%sMC", GetName()));
-  else if (fSumMC) fSumMC->Add((dataMC));
-  // Add contribution 
   
   PostData(1, fSums);
 }
@@ -378,7 +314,7 @@ AliBasedNdetaTask::ProjectX(const TH2D* h,
                            Int_t firstbin, 
                            Int_t lastbin, 
                            bool  corr,
-                           bool  error) const
+                           bool  error)
 {
   // 
   // Project onto the X axis 
@@ -394,9 +330,9 @@ AliBasedNdetaTask::ProjectX(const TH2D* h,
   //    Newly created histogram or null
   //
   if (!h) return 0;
-  //#if USE_ROOT_PROJECT
-  //return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
-  //#endif
+#if USE_ROOT_PROJECT
+  return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
+#endif
   
   TAxis* xaxis = h->GetXaxis();
   TAxis* yaxis = h->GetYaxis();
@@ -414,7 +350,8 @@ AliBasedNdetaTask::ProjectX(const TH2D* h,
   if      (last  < 0)                    last  = yaxis->GetNbins();
   else if (last  >  yaxis->GetNbins()+1) last  = yaxis->GetNbins();
   if (last-first < 0) { 
-    AliWarning(Form("Nothing to project [%d,%d]", first, last));
+    AliWarningGeneral("AliBasedNdetaTask", 
+                     Form("Nothing to project [%d,%d]", first, last));
     return 0;
     
   }
@@ -479,185 +416,45 @@ AliBasedNdetaTask::Terminate(Option_t *)
     return; 
   }
   
-  if (!fOutput) { 
-    fOutput = new TList;
-    fOutput->SetName(Form("%s_result", GetName()));
-    fOutput->SetOwner();
-  }
-
-  fSum     = static_cast<TH2D*>(fSums->FindObject(GetName()));
-  fSumMC   = static_cast<TH2D*>(fSums->FindObject(Form("%sMC", GetName())));
-  fTriggers       = static_cast<TH1D*>(fSums->FindObject("triggers"));
-
-  if (!fTriggers) { 
-    AliError("Couldn't find histogram 'triggers' in list");
-    return;
-  }
-  if (!fSum) { 
-    AliError("Couldn't find histogram 'base' in list");
-    return;
-  }
-  
-  Int_t nAll        = Int_t(fTriggers->GetBinContent(kAll));
-  Int_t nB          = Int_t(fTriggers->GetBinContent(kB));
-  Int_t nA          = Int_t(fTriggers->GetBinContent(kA));
-  Int_t nC          = Int_t(fTriggers->GetBinContent(kC));
-  Int_t nE          = Int_t(fTriggers->GetBinContent(kE));
-  Int_t nMB         = Int_t(fTriggers->GetBinContent(kMB));
-  Int_t nTriggered  = Int_t(fTriggers->GetBinContent(kWithTrigger));
-  Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
-  Int_t nAccepted   = Int_t(fTriggers->GetBinContent(kAccepted));
-  Int_t nGood       = nB - nA - nC + 2 * nE;
-  
-  //  Int_t nBg         = nA + nC -nE;
-  //std::cout<<"nBackground "<<nBg<<"   , with vertex in range "<<fNbgValid<<" ,  total "<<fNbgVertex<<"  "<<std::endl<<std::endl;
-  Float_t alpha            = ((Float_t)nAccepted) / ((Float_t)nWithVertex);
-  //Float_t alphaBG  = 1;
-  //if(fNbgVertex > 0) alphaBG= (Float_t)fNbgValid / (Float_t)fNbgVertex;
-  //Float_t nNormalizationJF = nAccepted + alpha*(nMB - nAccepted - nBg);
-  Float_t nNormalizationJF = nAccepted + alpha*(nTriggered - nWithVertex) ; // -alphaBG*(nBg-fNbgVertex);
-  
-  //Float_t nNormalizationBg = fNbgValid + alphaBG*nBg;
-  
-  //std::cout<<nTriggered -nWithVertex<<std::endl;
-  
-  //Double_t vtxEff   = Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
-  Double_t vtxEff   =  (Float_t)nAccepted / nNormalizationJF;
-  vtxEff = vtxEff*fTriggerEff;
-  
-  //if ( fTriggerMask == AliAODForwardMult::kNSD ) vtxEff = 0.97; //From paper...
-  //std::cout<<
-
-  //Double_t vtxEffBg   =  (Float_t)fNbgValid / nNormalizationBg;
-  //Double_t vtxEff   =  1.;//(Float_t)nAccepted / nNormalizationJF;
-  
-  
-  
-  
-  if (nTriggered <= 0) { 
-    AliError("Number of triggered events <= 0");
-    return;
-  }
-  if (nGood <= 0) { 
-    AliWarning(Form("Number of good events=%d=%d-%d-%d+2*%d<=0",
-                   nGood, nB, nA, nC, nE));
-    nGood = nMB;
-  }
-  //Double_t vtxEff   =  Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
-  
-  Double_t vNorm    = Double_t(nAccepted) / nGood;
-  AliInfo(Form("Total of %9d events\n"
-              "                   of these %9d are 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"
-              "                   Triggers by type:\n"
-              "                     B   = %9d\n"
-              "                     A|C = %9d (%9d+%-9d)\n"
-              "                     E   = %9d\n"
-              "                   Implies %9d good triggers\n"
-              "                   Vertex efficiency: %f (%f)",
-              nAll, nMB, nTriggered, 
-              AliAODForwardMult::GetTriggerString(fTriggerMask),
-              nWithVertex, nAccepted,
-              fVtxMin, fVtxMax, 
-              nB, nA+nC, nA, nC, nE, nGood, vtxEff, vNorm));
-  
-  const char* name = GetName();
-  
-  // Get acceptance normalisation from underflow bins 
-  TH1D* norm   = ProjectX(fSum, Form("norm%s",name), 0, 0, fCorrEmpty, false);
-  
-    
-  
+  fOutput = new TList;
+  fOutput->SetName(Form("%s_result", GetName()));
+  fOutput->SetOwner();
 
-  //
-  //std::cout<<norm->GetMaximumBin()<<"    "<< (Float_t)nAccepted / norm->GetBinContent((Float_t)norm->GetMaximumBin()) <<std::endl;
-  Float_t scaleForNorm =  (Float_t)nAccepted / (Float_t)norm->GetBinContent(norm->GetMaximumBin()) ;
-  //Float_t scaleForNorm =  norm->Integral() ;
-  std::cout<<norm->Integral()<<std::endl;
-  
-  norm->Scale(scaleForNorm);
-  //norm->Scale(1., "width");
-  //norm->Scale(norm->GetNbinsX() / (norm->GetXaxis()->GetXmax() - norm->GetXaxis()->GetXmin() ));
-  
-  //norm->Scale(1.,"width");
-  //
-  // Project onto eta axis - _ignoring_underflow_bins_!
-  
-  TH2D* tmpNorm = (TH2D*)fSum->Clone("tmpNorm");
-  
-  if(fShapeCorr)
-    fSum->Divide(fShapeCorr);
-  
-  TH1D* dndeta = ProjectX(fSum, Form("dndeta%s",name), 1, fSum->GetNbinsY(),
-                         fCorrEmpty);
-  
-  // Normalize to the acceptance 
-  //dndeta->Divide(norm);
-  
-  for(Int_t i = 1; i<=tmpNorm->GetNbinsX(); i++) {
-    
-    if( tmpNorm->GetBinContent(i,0) > 0 ) {
-      //std::cout<<tmpNorm->GetBinContent(i,0) <<std::endl;
-      dndeta->SetBinContent(i,dndeta->GetBinContent(i) / tmpNorm->GetBinContent(i,0));
-      dndeta->SetBinError(i,dndeta->GetBinError(i) / tmpNorm->GetBinContent(i,0));
-    }
-  }
-  
-  // Scale by the vertex efficiency 
-  //dndeta->Scale(vNorm, "width");
-    
-  dndeta->Scale(vtxEff, "width");
-  
-  //dndeta->Scale(1./(Float_t)nAccepted);
-  //dndeta->Scale(dndeta->GetNbinsX() / (dndeta->GetXaxis()->GetXmax() - dndeta->GetXaxis()->GetXmin() ));
-  
-  //norm->Scale(1. / nAccepted);
-  //norm->Scale(1. / nNormalizationJF);
-  SetHistogramAttributes(dndeta, kRed+1, 20, Form("ALICE %s", name));
-  SetHistogramAttributes(norm, kRed+1,20,Form("ALICE %s normalisation", name)); 
+  fSNNString = static_cast<TNamed*>(fSums->FindObject("sNN"));
+  fSysString = static_cast<TNamed*>(fSums->FindObject("sys"));
 
-  fOutput->Add(fTriggers->Clone());
-  if (fSymmetrice)   fOutput->Add(Symmetrice(dndeta));
-  fOutput->Add(dndeta);
-  fOutput->Add(norm);
-  fOutput->Add(Rebin(dndeta));
-  if (fSymmetrice)   fOutput->Add(Symmetrice(Rebin(dndeta)));
+  if(fSysString && fSNNString && 
+     fSysString->GetUniqueID() == AliForwardUtil::kPP)
+    LoadNormalizationData(fSysString->GetUniqueID(),
+                         fSNNString->GetUniqueID());
 
-  if (fSumMC) { 
-    // Get acceptance normalisation from underflow bins 
-    TH1D* normMC   = ProjectX(fSumMC,Form("norm%sMC", name), 0, 1, 
-                             fCorrEmpty, false);
-    // Project onto eta axis - _ignoring_underflow_bins_!
-    TH1D* dndetaMC = ProjectX(fSumMC,Form("dndeta%sMC", name),1,
-                             fSum->GetNbinsY(), fCorrEmpty);
-    // Normalize to the acceptance 
-    dndetaMC->Divide(normMC);
-    // Scale by the vertex efficiency 
-    dndetaMC->Scale(vNorm, "width");
-    normMC->Scale(1. / nAccepted);
+  // Loop over centrality bins 
+  TIter next(fListOfCentralities);
+  CentralityBin* bin = 0;
+  while ((bin = static_cast<CentralityBin*>(next()))) 
+    bin->End(fSums, fOutput, fUseShapeCorr ? fShapeCorr : 0, fTriggerEff,
+            fSymmetrice, fRebin, fCorrEmpty, fCutEdges, 
+            fVtxMin, fVtxMax, fTriggerMask);
 
-    SetHistogramAttributes(dndetaMC, kRed+3, 21, Form("ALICE %s (MC)",name));
-    SetHistogramAttributes(normMC, kRed+3, 21, Form("ALICE %s (MC) "
-                                                   "normalisation",name)); 
+  // Output collision energy string 
+  if (fSNNString) fOutput->Add(fSNNString->Clone());
 
-    fOutput->Add(dndetaMC);
-    if (fSymmetrice)   fOutput->Add(Symmetrice(dndetaMC));    
-    fOutput->Add(normMC);
-    fOutput->Add(Rebin(dndetaMC));
-    if (fSymmetrice)   fOutput->Add(Symmetrice(Rebin(dndetaMC)));
-  }
+  // Output collision system string 
+  if (fSysString) fOutput->Add(fSysString->Clone());
 
+  // Output trigger string 
   TNamed* trigString = 
     new TNamed("trigString", AliAODForwardMult::GetTriggerString(fTriggerMask));
   trigString->SetUniqueID(fTriggerMask);
   fOutput->Add(trigString);
 
+  // Output vertex axis 
   TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
   vtxAxis->SetName("vtxAxis");
   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
   fOutput->Add(vtxAxis);
+
+  // Output trigger efficiency and shape correction 
   fOutput->Add(new TParameter<Double_t>("triggerEff", fTriggerEff)); 
   if (fShapeCorr) fOutput->Add(fShapeCorr);
   
@@ -667,49 +464,49 @@ AliBasedNdetaTask::Terminate(Option_t *)
 void
 AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
 {
-  sys = 1;
-  energy = 900;
+  // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD
   TString type("pp");
-  if(sys == 2) type.Form("PbPb");
   TString snn("900");
   if(energy == 7000) snn.Form("7000");
   if(energy == 2750) snn.Form("2750"); 
   
-  if(fShapeCorr && (fTriggerEff != 1)) {AliInfo("Objects already set for normalization - no action taken"); return; }
-  
-  TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWG2/FORWARD/corrections/Normalization/normalizationHists_%s_%s.root",type.Data(),snn.Data()));
-  if(!fin) AliWarning("no file for normalization");
-  
-  if ( fTriggerMask == AliAODForwardMult::kInel ) {
-    TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get("hInelNormalization"));
-    if(shapeCor) SetShapeCorrection(shapeCor);
-    TParameter<float>* ineleff = (TParameter<float>*)fin->Get("inelTriggerEff");
-  if ( ineleff) SetTriggerEff(ineleff->GetVal());
-  if (shapeCor && ineleff) AliInfo("Loaded objects for normalization.");
+  if(fShapeCorr && (fTriggerEff != 1)) {
+    AliInfo("Objects already set for normalization - no action taken"); 
+    return; 
+  }
   
-}
+  TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWG2/FORWARD/corrections/"
+                               "Normalization/normalizationHists_%s_%s.root",
+                               type.Data(),snn.Data()));
+  if(!fin) {
+    AliWarning(Form("no file for normalization of %d/%d", sys, energy));
+    return;
+  }
 
-if ( fTriggerMask == AliAODForwardMult::kNSD ) {
-  TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get("hNSDNormalization"));
-  if(shapeCor) SetShapeCorrection(shapeCor);
-  TParameter<float>* nsdeff = (TParameter<float>*)fin->Get("nsdTriggerEff");
-if ( nsdeff) SetTriggerEff(nsdeff->GetVal());
-if (shapeCor && nsdeff) AliInfo("Loaded objects for normalization.");
-}
+  // Shape correction 
+  TString shapeCorName(Form("h%sNormalization", 
+                           fTriggerMask == AliAODForwardMult::kInel ? "Inel" :
+                           fTriggerMask == AliAODForwardMult::kNSD ? "NSD" :
+                           fTriggerMask == AliAODForwardMult::kInelGt0 ?
+                           "InelGt0" : "All"));
+  TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName));
+  if (shapeCor) SetShapeCorrection(shapeCor);
 
-if ( fTriggerMask == AliAODForwardMult::kInelGt0 ) {
-  TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get("hInelGt0Normalization"));
-    if(shapeCor) SetShapeCorrection(shapeCor);
-    TParameter<float>* inelgt0eff = (TParameter<float>*)fin->Get("inelgt0TriggerEff");
-if ( inelgt0eff) SetTriggerEff(inelgt0eff->GetVal());
-if (shapeCor && inelgt0eff) AliInfo("Loaded objects for normalization.");
-  }
-  
+  // Trigger efficiency
+  TString effName(Form("%sTriggerEff", 
+                      fTriggerMask == AliAODForwardMult::kInel ? "inel" :
+                      fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
+                      fTriggerMask == AliAODForwardMult::kInelGt0 ?
+                      "inelgt0" : "all"));
+  TParameter<float>* eff = static_cast<TParameter<float>*>(fin->Get(effName));
+  if (eff) SetTriggerEff(eff->GetVal());
 
+  // Print - out
+  if (shapeCor && eff) AliInfo("Loaded objects for normalization.");
 }
 //________________________________________________________________________
 TH1D*
-AliBasedNdetaTask::Rebin(const TH1D* h) const
+AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges)
 {
   // 
   // Make a copy of the input histogram and rebin that histogram
@@ -720,40 +517,45 @@ AliBasedNdetaTask::Rebin(const TH1D* h) const
   // Return:
   //    New (rebinned) histogram
   //
-  if (fRebin <= 1) return 0;
+  if (rebin <= 1) return 0;
 
   Int_t nBins = h->GetNbinsX();
-  if(nBins % fRebin != 0) {
-    Warning("Rebin", "Rebin factor %d is not a devisor of current number "
-           "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
+  if(nBins % rebin != 0) {
+    AliWarningGeneral("AliBasedNdetaTask", 
+                     Form("Rebin factor %d is not a devisor of current number "
+                          "of bins %d in the histogram %s", 
+                          rebin, nBins, h->GetName()));
     return 0;
   }
     
   // Make a copy 
   TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d", 
-                                              h->GetName(), fRebin)));
-  tmp->Rebin(fRebin);
+                                              h->GetName(), rebin)));
+  tmp->Rebin(rebin);
   tmp->SetDirectory(0);
 
   // The new number of bins 
-  Int_t nBinsNew = nBins / fRebin;
+  Int_t nBinsNew = nBins / rebin;
   for(Int_t i = 1;i<= nBinsNew; i++) {
     Double_t content = 0;
     Double_t sumw    = 0;
     Double_t wsum    = 0;
     Int_t    nbins   = 0;
-    for(Int_t j = 1; j<=fRebin;j++) {
-      Int_t    bin = (i-1)*fRebin + j;
+    for(Int_t j = 1; j<=rebin;j++) {
+      Int_t    bin = (i-1)*rebin + j;
       Double_t c   =  h->GetBinContent(bin);
       if (c <= 0) continue;
       
-      if (fCutEdges) {
+      if (cutEdges) {
        if (h->GetBinContent(bin+1)<=0 || 
            h->GetBinContent(bin-1)<=0) {
-         Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)", 
-                 bin, c, h->GetName(), 
-                 bin+1, h->GetBinContent(bin+1), 
-                 bin-1, h->GetBinContent(bin-1));
+#if 0
+         AliWarningGeneral("AliBasedNdetaTask", 
+                           Form("removing bin %d=%f of %s (%d=%f,%d=%f)", 
+                                bin, c, h->GetName(), 
+                                bin+1, h->GetBinContent(bin+1), 
+                                bin-1, h->GetBinContent(bin-1)));
+#endif
          continue;
        }       
       }
@@ -776,7 +578,7 @@ AliBasedNdetaTask::Rebin(const TH1D* h) const
 
 //__________________________________________________________________
 TH1* 
-AliBasedNdetaTask::Symmetrice(const TH1* h) const
+AliBasedNdetaTask::Symmetrice(const TH1* h)
 {
   // 
   // Make an extension of @a h to make it symmetric about 0 
@@ -820,6 +622,419 @@ AliBasedNdetaTask::Symmetrice(const TH1* h) const
   s->SetBinError(l2+1, h->GetBinError(first));
   return s;
 }
+
+//====================================================================
+AliBasedNdetaTask::CentralityBin::CentralityBin()
+  : TNamed("", ""), 
+    fSums(0), 
+    fOutput(0),
+    fSum(0), 
+    fSumMC(0), 
+    fTriggers(0), 
+    fLow(0), 
+    fHigh(0)
+{
+  // 
+  // Constructor 
+  //
+}
+//____________________________________________________________________
+AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name, 
+                                               Short_t low, Short_t high)
+  : TNamed(name, ""), 
+    fSums(0), 
+    fOutput(0),
+    fSum(0), 
+    fSumMC(0), 
+    fTriggers(0),
+    fLow(low), 
+    fHigh(high)
+{
+  // 
+  // Constructor 
+  // 
+  // Parameters:
+  //    name Name used for histograms (e.g., Forward)
+  //    low  Lower centrality cut in percent 
+  //    high Upper centrality cut in percent 
+  //
+  if (low <= 0 && high <= 0) { 
+    fLow  = 0; 
+    fHigh = 0;
+    SetTitle("All centralities");
+  }
+  else {
+    fLow  = low;
+    fHigh = high;
+    SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high));
+  }
+}
+//____________________________________________________________________
+AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o)
+  : TNamed(o), 
+    fSums(o.fSums), 
+    fOutput(o.fOutput),
+    fSum(o.fSum), 
+    fSumMC(o.fSumMC), 
+    fTriggers(o.fTriggers), 
+    fLow(o.fLow), 
+    fHigh(o.fHigh)
+{
+  // 
+  // Copy constructor 
+  // 
+  // Parameters:
+  //    other Object to copy from 
+  //
+}
+//____________________________________________________________________
+AliBasedNdetaTask::CentralityBin::~CentralityBin()
+{
+  // 
+  // Destructor 
+  //
+  if (fSums) fSums->Delete();
+  if (fOutput) fOutput->Delete();
+}
+
+//____________________________________________________________________
+AliBasedNdetaTask::CentralityBin&
+AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
+{
+  // 
+  // Assignment operator 
+  // 
+  // Parameters:
+  //    other Object to assign from 
+  // 
+  // Return:
+  //    Reference to this 
+  //
+  SetName(o.GetName());
+  SetTitle(o.GetTitle());
+  fSums      = o.fSums;
+  fOutput    = o.fOutput;
+  fSum       = o.fSum;
+  fSumMC     = o.fSumMC;
+  fTriggers  = o.fTriggers;
+  fLow       = o.fLow;
+  fHigh      = o.fHigh;
+
+  return *this;
+}
+//____________________________________________________________________
+const char* 
+AliBasedNdetaTask::CentralityBin::GetListName() const
+{
+  // 
+  // Get the list name 
+  // 
+  // Return:
+  //    List Name 
+  //
+  if (IsAllBin()) return "all"; 
+  return Form("cent%03d_%03d", fLow, fHigh);
+}
+//____________________________________________________________________
+void
+AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir)
+{
+  // 
+  // Create output objects 
+  // 
+  // Parameters:
+  //    dir   Parent list
+  //
+  fSums = new TList;
+  fSums->SetName(GetListName());
+  fSums->SetOwner();
+  dir->Add(fSums);
+
+  fTriggers = new TH1D("triggers", "Number of triggers", 
+                      kAccepted, 1, kAccepted);
+  fTriggers->SetYTitle("# of events");
+  fTriggers->GetXaxis()->SetBinLabel(kAll,         "All events");
+  fTriggers->GetXaxis()->SetBinLabel(kB,           "w/B trigger");
+  fTriggers->GetXaxis()->SetBinLabel(kA,           "w/A trigger");
+  fTriggers->GetXaxis()->SetBinLabel(kC,           "w/C trigger");
+  fTriggers->GetXaxis()->SetBinLabel(kE,           "w/E trigger");
+  fTriggers->GetXaxis()->SetBinLabel(kMB,          "w/Collision trigger");
+  fTriggers->GetXaxis()->SetBinLabel(kWithVertex,  "w/Vertex");
+  fTriggers->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger");
+  fTriggers->GetXaxis()->SetBinLabel(kPileUp,      "w/Pileup");
+  fTriggers->GetXaxis()->SetBinLabel(kAccepted,    "Accepted by cut");
+  fTriggers->GetXaxis()->SetNdivisions(kAccepted, false);
+  fTriggers->SetFillColor(kRed+1);
+  fTriggers->SetFillStyle(3001);
+  fTriggers->SetStats(0);
+  fSums->Add(fTriggers);
+}
+//____________________________________________________________________
+void
+AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc)
+{
+  // 
+  // Create sum histogram 
+  // 
+  // Parameters:
+  //    data  Data histogram to clone 
+  //    mc    (optional) MC histogram to clone 
+  //
+  fSum = static_cast<TH2D*>(data->Clone(GetName()));
+  fSum->SetDirectory(0);
+  fSum->Reset();
+  fSums->Add(fSum);
+  
+  // If no MC data is given, then do not create MC sum histogram 
+  if (!mc) return;
+
+  fSumMC = static_cast<TH2D*>(mc->Clone(Form("%sMC", GetName())));
+  fSumMC->SetDirectory(0);
+  fSumMC->Reset();
+  fSums->Add(fSumMC);
+}
+
+//____________________________________________________________________
+Bool_t
+AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
+                                            Int_t triggerMask,
+                                            Double_t vzMin, Double_t vzMax)
+{
+  // 
+  // Check the trigger, vertex, and centrality
+  // 
+  // Parameters:
+  //    aod Event input 
+  // 
+  // Return:
+  //    true if the event is to be used 
+  //
+  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;
+  }
+  
+  fTriggers->AddBinContent(kAll);
+  if (forward->IsTriggerBits(AliAODForwardMult::kB)) 
+    fTriggers->AddBinContent(kB);
+  if (forward->IsTriggerBits(AliAODForwardMult::kA)) 
+    fTriggers->AddBinContent(kA);
+  if (forward->IsTriggerBits(AliAODForwardMult::kC)) 
+    fTriggers->AddBinContent(kC);
+  if (forward->IsTriggerBits(AliAODForwardMult::kE)) 
+    fTriggers->AddBinContent(kE);
+  if (forward->IsTriggerBits(AliAODForwardMult::kInel)) 
+    fTriggers->AddBinContent(kMB);
+  if (forward->IsTriggerBits(AliAODForwardMult::kMCNSD)) 
+    fTriggers->AddBinContent(kMCNSD);
+  if (forward->IsTriggerBits(AliAODForwardMult::kPileUp)) 
+    fTriggers->AddBinContent(kPileUp);
+  
+  // Check if we have an event of interest. 
+  if (!forward->IsTriggerBits(triggerMask)) return false;
+  
+  //Check for pileup
+  if (forward->IsTriggerBits(AliAODForwardMult::kPileUp)) return false;
+  
+  fTriggers->AddBinContent(kWithTrigger);
+  
+  // Check that we have a valid vertex
+  if (!forward->HasIpZ()) return false;
+  fTriggers->AddBinContent(kWithVertex);
+
+  // Check that vertex is within cuts 
+  if (!forward->InRange(vzMin, vzMax)) return false;
+  fTriggers->AddBinContent(kAccepted);
+  
+  return true;
+}
+  
+  
+//____________________________________________________________________
+void
+AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward,
+                                              Int_t triggerMask,
+                                              Double_t vzMin, Double_t vzMax,
+                                              const TH2D* data, const TH2D* mc)
+{
+  // 
+  // Process an event
+  // 
+  // Parameters:
+  //    forward     Forward data (for trigger, vertex, & centrality)
+  //    triggerMask Trigger mask 
+  //    vzMin       Minimum IP z coordinate
+  //    vzMax       Maximum IP z coordinate
+  //    data        Data histogram 
+  //    mc          MC histogram
+  //
+  if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return;
+  if (!data) return;
+  if (!fSum) CreateSums(data, mc);
+  
+  fSum->Add(data);
+  if (mc) fSumMC->Add(mc);
+}
+
+//________________________________________________________________________
+void 
+AliBasedNdetaTask::CentralityBin::End(TList*      sums, 
+                                     TList*      results,
+                                     const TH1*  shapeCorr, 
+                                     Double_t    trigEff,
+                                     Bool_t      symmetrice,
+                                     Int_t       rebin, 
+                                     Bool_t      corrEmpty, 
+                                     Bool_t      cutEdges, 
+                                     Double_t    vzMin, 
+                                     Double_t    vzMax, 
+                                     Int_t       triggerMask) 
+{
+  // 
+  // End of processing 
+  // 
+  // Parameters:
+  //    sums        List of sums
+  //    results     Output list of results
+  //    shapeCorr   Shape correction or nil
+  //    trigEff     Trigger efficiency 
+  //    symmetrice  Whether to symmetrice the results
+  //    rebin       Whether to rebin the results
+  //    corrEmpty   Whether to correct for empty bins
+  //    cutEdges    Whether to cut edges when rebinning
+  //    vzMin       Minimum IP z coordinate
+  //    vzMax    Maximum IP z coordinate
+  //    triggerMask Trigger mask 
+  //
+
+  fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
+  if(!fSums) {
+    AliError("Could not retrieve TList fSums"); 
+    return; 
+  }
+  
+  fOutput = new TList;
+  fOutput->SetName(GetListName());
+  fOutput->SetOwner();
+  results->Add(fOutput);
+
+  fSum      = static_cast<TH2D*>(fSums->FindObject(GetName()));
+  fSumMC    = static_cast<TH2D*>(fSums->FindObject(Form("%sMC", GetName())));
+  fTriggers = static_cast<TH1D*>(fSums->FindObject("triggers"));
+
+  if (!fTriggers) { 
+    AliError("Couldn't find histogram 'triggers' in list");
+    return;
+  }
+  if (!fSum) { 
+    AliError("Couldn't find histogram 'base' in list");
+    return;
+  }
+  
+  Int_t nAll        = Int_t(fTriggers->GetBinContent(kAll));
+  Int_t nB          = Int_t(fTriggers->GetBinContent(kB));
+  Int_t nA          = Int_t(fTriggers->GetBinContent(kA));
+  Int_t nC          = Int_t(fTriggers->GetBinContent(kC));
+  Int_t nE          = Int_t(fTriggers->GetBinContent(kE));
+  Int_t nMB         = Int_t(fTriggers->GetBinContent(kMB));
+  Int_t nTriggered  = Int_t(fTriggers->GetBinContent(kWithTrigger));
+  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) { 
+    AliError("Number of triggered events <= 0");
+    return;
+  }
+  if (nGood <= 0) { 
+    AliWarning(Form("Number of good events=%d=%d-%d-%d+2*%d<=0",
+                   nGood, nB, nA, nC, nE));
+    nGood = nMB;
+  }
+  AliInfo(Form("Total of %9d events for %s\n"
+              "                   of these %9d are 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"
+              "                   Triggers by type:\n"
+              "                     B   = %9d\n"
+              "                     A|C = %9d (%9d+%-9d)\n"
+              "                     E   = %9d\n"
+              "                   Implies %9d good triggers\n"
+              "                   Vertex efficiency: %f (%f)",
+              nAll, GetTitle(), nMB, nTriggered, 
+              AliAODForwardMult::GetTriggerString(triggerMask),
+              nWithVertex, nAccepted,
+              vzMin, vzMax, 
+              nB, nA+nC, nA, nC, nE, nGood, vtxEff, vNorm));
+  
+  const char* name = GetName();
+  
+  // Get acceptance normalisation from underflow bins 
+  TH1D* norm = ProjectX(fSum, Form("norm%s",name), 0, 0, corrEmpty, false);
+  norm->SetDirectory(0);
+
+  // Scale by shape correction 
+  if(shapeCorr) fSum->Divide(shapeCorr);
+  else AliInfo("No shape correction specified, or disabled");
+  
+  // Project on X axis 
+  TH1D* dndeta = ProjectX(fSum, Form("dndeta%s",name), 1, fSum->GetNbinsY(),
+                         corrEmpty);
+  dndeta->SetDirectory(0);
+
+  // Normalize to the acceptance -
+  dndeta->Divide(norm);
+  
+  // Scale by the vertex efficiency 
+  dndeta->Scale(vtxEff, "width");
+  
+  SetHistogramAttributes(dndeta, kRed+1, 20, Form("ALICE %s", name));
+  SetHistogramAttributes(norm, kRed+1,20,Form("ALICE %s normalisation", name)); 
+
+  fOutput->Add(fTriggers->Clone());
+  if (symmetrice)   fOutput->Add(Symmetrice(dndeta));
+  fOutput->Add(dndeta);
+  fOutput->Add(norm);
+  fOutput->Add(Rebin(dndeta, rebin, cutEdges));
+  if (symmetrice)   fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges)));
+
+  if (fSumMC) { 
+    // Get acceptance normalisation from underflow bins 
+    TH1D* normMC   = ProjectX(fSumMC,Form("norm%sMC", name), 0, 1, 
+                             corrEmpty, false);
+    // Project onto eta axis - _ignoring_underflow_bins_!
+    TH1D* dndetaMC = ProjectX(fSumMC,Form("dndeta%sMC", name),1,
+                             fSum->GetNbinsY(), corrEmpty);
+    // Normalize to the acceptance 
+    dndetaMC->Divide(normMC);
+    // Scale by the vertex efficiency 
+    dndetaMC->Scale(vNorm, "width");
+    normMC->Scale(1. / nAccepted);
+
+    SetHistogramAttributes(dndetaMC, kRed+3, 21, Form("ALICE %s (MC)",name));
+    SetHistogramAttributes(normMC, kRed+3, 21, Form("ALICE %s (MC) "
+                                                   "normalisation",name)); 
+
+    fOutput->Add(dndetaMC);
+    if (symmetrice)   fOutput->Add(Symmetrice(dndetaMC));    
+
+    fOutput->Add(normMC);
+    fOutput->Add(Rebin(dndetaMC, rebin, cutEdges));
+
+    if (symmetrice)   
+      fOutput->Add(Symmetrice(Rebin(dndetaMC, rebin, cutEdges)));
+  }
+}
+
 //
 // EOF
 //
index e0cbea3..fd7bea7 100644 (file)
@@ -4,11 +4,11 @@
 #ifndef ALIBASEDNDETATASK_H
 #define ALIBASEDNDETATASK_H
 #include <AliAnalysisTaskSE.h>
-// #include <AliAODBaseMult.h>
 class TList;
 class TH2D;
 class TH1D;
 class AliAODEvent;
+class AliAODForwardMult;
 
 /**
  * Task to determine the 
@@ -28,8 +28,24 @@ public:
    * @param maxVtx  Set @f$v_z@f$ range
    */
   AliBasedNdetaTask(const char* name);
+  /**
+   * Destructor
+   * 
+   */
+  virtual ~AliBasedNdetaTask();
 
   /** 
+   * @{ 
+   * @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)
@@ -74,30 +90,23 @@ public:
    */
   void SetShapeCorrection(const TH1* h);
   /** 
-   * Load the normalization data - done automatically if not set from outside
-   * 
-   * @param sys system
-   * @param energy energy
-   */
-  void LoadNormalizationData(UShort_t sys, UShort_t energy);
-  /** 
-   * Load the normalization data - done automatically if not set from outside
-   * 
-   * @param centlow low cent cut
+   * Set whether to use the shape correction 
+   *
+   * @param use  whether to use the shape correction 
    */
-  void SetCentLow(Float_t centlow)   { fCentLow = centlow; }
+  void SetUseShapeCorrection(Bool_t use) { fUseShapeCorr = use; }
   /** 
    * Load the normalization data - done automatically if not set from outside
    * 
-   * @param centhigh high cent cut
+   * @param sys system
+   * @param energy energy
    */
-  void SetCentHigh(Float_t centhigh) { fCentHigh = centhigh; }
-  
-  /**
-   * Destructor
-   * 
+  void LoadNormalizationData(UShort_t sys, UShort_t energy);  
+  /** @} */
+
+  /** @{ 
+   *  @name Task interface 
    */
-  virtual ~AliBasedNdetaTask();
   /** 
    * Initialise on master - does nothing
    * 
@@ -123,36 +132,12 @@ public:
    * @param option Not used 
    */
   virtual void Terminate(Option_t* option);
-protected:
-  AliBasedNdetaTask(const AliBasedNdetaTask&);
-  AliBasedNdetaTask& operator=(const AliBasedNdetaTask&) { return *this; }
+  /** @} */
 
   /** 
-   * Retrieve the histogram 
-   * 
-   * @param aod AOD event 
-   * @param mc  Whether to get the MC histogram or not
-   * 
-   * @return Retrieved histogram or null
-   */
-  virtual TH2D* GetHistogram(const AliAODEvent* aod, Bool_t mc=false) = 0;
-  /** 
-   * Check the trigger and vertex 
-   * 
-   * @param aod 
-   * 
-   * @return 
+   * @{ 
+   * @name Services member functions 
    */
-  Bool_t CheckEvent(const AliAODEvent* aod);
-  /** 
-   * Clone a 2D histogram
-   * 
-   * @param in    Histogram to clone.
-   * @param name  New name of clone.
-   * 
-   * @return The clone
-   */
-  TH2D* CloneHist(const TH2D* in, const char* name);
   /** 
    * Make a copy of the input histogram and rebin that histogram
    * 
@@ -160,7 +145,7 @@ protected:
    * 
    * @return New (rebinned) histogram
    */
-  TH1D* Rebin(const TH1D* h) const;
+  static TH1D* Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges=false);
   /** 
    * Make an extension of @a h to make it symmetric about 0 
    * 
@@ -168,7 +153,7 @@ protected:
    * 
    * @return Symmetric extension of @a h 
    */
-  TH1* Symmetrice(const TH1* h) const;
+  static TH1* Symmetrice(const TH1* h);
   /** 
    * Project onto the X axis 
    * 
@@ -180,12 +165,12 @@ protected:
    * 
    * @return Newly created histogram or null
    */
-  TH1D* ProjectX(const TH2D* h, 
-                const char* name,
-                Int_t firstbin, 
-                Int_t lastbin, 
-                bool  corr=true,
-                bool  error=true) const;
+  static TH1D* ProjectX(const TH2D* h, 
+                       const char* name,
+                       Int_t firstbin, 
+                       Int_t lastbin, 
+                       bool  corr=true,
+                       bool  error=true);
   /** 
    * Set histogram graphical options, etc. 
    * 
@@ -195,9 +180,35 @@ protected:
    * @param title   Title of histogram
    * @param ytitle  Title on y-axis. 
    */
-  void  SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker, 
-                              const char* title, 
-                              const char* ytitle="#frac{1}{N} #frac{dN_{ch}}{d#eta}");
+  static void SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker, 
+                                    const char* title, 
+                                    const char* ytitle="#frac{1}{N} #frac{dN_{ch}}{d#eta}");
+  /** @} */
+protected:
+  AliBasedNdetaTask(const AliBasedNdetaTask&);
+  AliBasedNdetaTask& operator=(const AliBasedNdetaTask&) { return *this; }
+  class CentralityBin;
+
+  /** 
+   * Retrieve the histogram 
+   * 
+   * @param aod AOD event 
+   * @param mc  Whether to get the MC histogram or not
+   * 
+   * @return Retrieved histogram or null
+   */
+  virtual TH2D* GetHistogram(const AliAODEvent* aod, Bool_t mc=false) = 0;
+  /** 
+   * Make a centrality bin 
+   * 
+   * @param name  Name used for histograms
+   * @param low   Low cut in percent
+   * @param high  High cut in percent
+   * 
+   * @return A newly created centrality bin 
+   */
+  virtual CentralityBin* MakeCentralityBin(const char* name, Short_t low, 
+                                          Short_t high) const;
   /** 
    * Trigger histogram bins 
    */
@@ -214,15 +225,163 @@ protected:
     kAccepted   = 10,
     kMCNSD      = 11
   };
+  /**
+   * Calculations done per centrality 
+   * 
+   */
+  struct CentralityBin : public TNamed
+  {
+    /** 
+     * Constructor 
+     */
+    CentralityBin();
+    /** 
+     * Constructor 
+     * 
+     * @param name Name used for histograms (e.g., Forward)
+     * @param low  Lower centrality cut in percent 
+     * @param high Upper centrality cut in percent 
+     */
+    CentralityBin(const char* name, Short_t low, Short_t high);
+    /** 
+     * Copy constructor 
+     * 
+     * @param other Object to copy from 
+     */
+    CentralityBin(const CentralityBin& other);
+    /** 
+     * Destructor 
+     */
+    virtual ~CentralityBin();
+    /** 
+     * Assignment operator 
+     * 
+     * @param other Object to assign from 
+     * 
+     * @return Reference to this 
+     */
+    CentralityBin& operator=(const CentralityBin& other);
+    /** 
+     * Check if this is the 'all' bin 
+     * 
+     * @return true if low and high cuts are both zero
+     */    
+    Bool_t IsAllBin() const { return fLow == 0 && fHigh == 0; }
+    /** 
+     * Get the list name 
+     * 
+     * @return List Name 
+     */
+    const char* GetListName() const;
+    /** 
+     * Create output objects 
+     * 
+     * @param dir   Parent list
+     */
+    virtual void CreateOutputObjects(TList* dir);
+    /** 
+     * Process an event
+     * 
+     * @param forward     Forward data (for trigger, vertex, & centrality)
+     * @param triggerMask Trigger mask 
+     * @param vzMin       Minimum IP z coordinate
+     * @param vzMax       Maximum IP z coordinate
+     * @param data        Data histogram 
+     * @param mc          MC histogram
+     */
+    virtual void ProcessEvent(const AliAODForwardMult* forward, 
+                             Int_t triggerMask,
+                             Double_t vzMin, Double_t vzMax, 
+                             const TH2D* data, const TH2D* mc);
+    /** 
+     * End of processing 
+     * 
+     * @param sums        List of sums
+     * @param results     Output list of results
+     * @param shapeCorr   Shape correction or nil
+     * @param trigEff     Trigger efficiency 
+     * @param symmetrice  Whether to symmetrice the results
+     * @param rebin       Whether to rebin the results
+     * @param corrEmpty   Whether to correct for empty bins
+     * @param cutEdges    Whether to cut edges when rebinning
+     * @param vzMin       Minimum IP z coordinate
+     * @param vzMax      Maximum IP z coordinate
+     * @param triggerMask Trigger mask 
+     */
+    virtual void End(TList*      sums, 
+                    TList*      results,
+                    const TH1*  shapeCorr, 
+                    Double_t    trigEff,
+                    Bool_t      symmetrice,
+                    Int_t       rebin, 
+                    Bool_t      corrEmpty, 
+                    Bool_t      cutEdges, 
+                    Double_t    vzMin, 
+                    Double_t    vzMax, 
+                    Int_t       triggerMask);
+    /**
+     * @{
+     * @name Access histograms
+     */
+    /** 
+     * Get sum histogram 
+     * 
+     * @param mc If true, return MC histogram 
+     * 
+     * @return Sum histogram
+     */
+    const TH2D* GetSum(Bool_t mc=false) const { return mc ? fSumMC : fSum; }
+    /** 
+     * Get sum histogram 
+     * 
+     * @param mc If true, return MC histogram 
+     * 
+     * @return Sum histogram
+     */
+    TH2D* GetSum(Bool_t mc=false) { return mc ? fSumMC : fSum; }
+    /** 
+     * Get trigger histogram
+     * 
+     * @return Trigger histogram
+     */
+    const TH1D* GetTriggers() const { return fTriggers; } 
+    /** 
+     * Get trigger histogram
+     * 
+     * @return Trigger histogram 
+     */
+    TH1D* GetTrigggers() { return fTriggers; }
+    /** @} */
+  protected:
+    /** 
+     * Create sum histogram 
+     * 
+     * @param data  Data histogram to clone 
+     * @param mc    (optional) MC histogram to clone 
+     */
+    virtual void CreateSums(const TH2D* data, const TH2D* mc);
+    /** 
+     * Check the trigger, vertex, and centrality
+     * 
+     * @param aod Event input 
+     * 
+     * @return true if the event is to be used 
+     */
+    virtual Bool_t CheckEvent(const AliAODForwardMult* forward, 
+                             Int_t triggerMask,
+                             Double_t vzMin, Double_t vzMax);
+    TList*   fSums;      // Output list 
+    TList*   fOutput;    // Output list 
+    TH2D*    fSum;       // Sum histogram
+    TH2D*    fSumMC;     // MC sum histogram
+    TH1D*    fTriggers;  // Trigger histogram 
+    UShort_t fLow;       // Lower limit (inclusive)
+    UShort_t fHigh;      // Upper limit (exclusive)
 
-  TH2D*           fSum;          // Sum of histograms 
-  TH2D*           fSumMC;        // Sum of MC histograms (if any)
-
+    ClassDef(CentralityBin,1); // A centrality bin 
+  };
   TList*          fSums;         // Container of sums 
   TList*          fOutput;       // Container of outputs 
-
-  TH1D*           fTriggers;     // Histogram of triggers 
-   
   Double_t        fVtxMin;       // Minimum v_z
   Double_t        fVtxMax;       // Maximum v_z
   Int_t           fTriggerMask;  // Trigger mask 
@@ -232,10 +391,13 @@ protected:
   Bool_t          fCorrEmpty;    // Correct for empty bins 
   Double_t        fTriggerEff;   // Trigger efficiency for selected trigger(s)
   TH1*            fShapeCorr;    // Shape correction 
-  Float_t         fCentLow;      // Low centrality cut
-  Float_t         fCentHigh;      // High centrality cut
-  
-  ClassDef(AliBasedNdetaTask,1); // Determine multiplicity in base area
+  TList*          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 
+
+  ClassDef(AliBasedNdetaTask,2); // Determine multiplicity in base area
 };
 
 #endif
index a51d8b5..4adcaeb 100644 (file)
 #include <TROOT.h>
 #include <iostream>
 #include <iomanip>
-#include "AliMCEvent.h"
-#include "AliGenPythiaEventHeader.h"
-#include "AliGenDPMjetEventHeader.h"
-#include "AliHeader.h"
-#include "AliMCEventHandler.h"
+
 //====================================================================
 AliFMDEventInspector::AliFMDEventInspector()
   : TNamed(),
@@ -317,49 +313,38 @@ AliFMDEventInspector::Process(const AliESDEvent* event,
   //    0 (or kOk) on success, otherwise a bit mask of error codes 
   //
 
-  // Check that we have an event 
+  // --- Check that we have an event ---------------------------------
   if (!event) { 
     AliWarning("No ESD event found for input event");
     return kNoEvent;
   }
 
-  // Read trigger information from the ESD and store in AOD object
+  // --- Read trigger information from the ESD and store in AOD object
   if (!ReadTriggers(event, triggers)) { 
     if (fDebug > 2) {
       AliWarning("Failed to read triggers from ESD"); }
     return kNoTriggers;
   }
 
-  // Check if this is a high-flux event 
+  // --- Check if this is a high-flux event --------------------------
   const AliMultiplicity* testmult = event->GetMultiplicity();
   if (!testmult) {
     if (fDebug > 3) {
       AliWarning("No central multiplicity object found"); }
-    return kNoSPD;
   }
-  lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
+  else 
+    lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
 
   fHType->Fill(lowFlux ? 0 : 1);
   
+  // --- Read centrality information 
   cent = -10;
-  AliCentrality* centObj = const_cast<AliESDEvent*>(event)->GetCentrality();
-  if (centObj) {
-    // AliInfo(Form("Got centrality object %p with quality %d", 
-    //              centObj, centObj->GetQuality()));
-    // centObj->Print();
-    cent = centObj->GetCentralityPercentileUnchecked("V0M");  
-  }
-  // AliInfo(Form("Centrality is %f", cent));
-  fHCent->Fill(cent);
-
-  // Check the FMD ESD data 
-  if (!event->GetFMDData()) { 
-    if (fDebug > 3) {
-      AliWarning("No FMD data found in ESD"); }
-    return kNoFMD;
+  if (!ReadCentrality(event, cent)) {
+    if (fDebug > 3) 
+      AliWarning("Failed to get centrality");
   }
 
-  // Get the vertex information 
+  // --- Get the vertex information ----------------------------------
   vz          = 0;
   Bool_t vzOk = ReadVertex(event, vz);
 
@@ -371,7 +356,7 @@ AliFMDEventInspector::Process(const AliESDEvent* event,
   }
   fHEventsTrVtx->Fill(vz);
 
-  // Get the vertex bin 
+  // --- Get the vertex bin ------------------------------------------
   ivz = fHEventsTr->GetXaxis()->FindBin(vz);
   if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) { 
     if (fDebug > 3) {
@@ -382,12 +367,46 @@ AliFMDEventInspector::Process(const AliESDEvent* event,
     return kBadVertex;
   }
   
+  // --- Check the FMD ESD data --------------------------------------
+  if (!event->GetFMDData()) { 
+    if (fDebug > 3) {
+      AliWarning("No FMD data found in ESD"); }
+    return kNoFMD;
+  }
+
   
   return kOk;
 }
 
 //____________________________________________________________________
 Bool_t
+AliFMDEventInspector::ReadCentrality(const AliESDEvent* esd, Double_t& cent)
+{
+  // 
+  // Read centrality from event 
+  // 
+  // Parameters:
+  //    esd  Event 
+  //    cent On return, the centrality or negative if not found
+  // 
+  // Return:
+  //    False on error, true otherwise 
+  //
+  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->GetCentralityPercentileUnchecked("V0M");  
+  }
+  // AliInfo(Form("Centrality is %f", cent));
+  fHCent->Fill(cent);
+
+  return true;
+}
+
+//____________________________________________________________________
+Bool_t
 AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
 {
   // 
@@ -416,39 +435,8 @@ AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
     AliWarning("No input handler");
     return kFALSE;
   }
-  
-  AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-  AliMCEvent* mcEvent = 0;
-  if(mcHandler)
-    mcEvent = mcHandler->MCEvent();
-  
-  if(mcEvent) {
-    //Assign MC only triggers : True NSD etc.
-    AliHeader* header            = mcEvent->Header();
-    AliGenEventHeader* genHeader = header->GenEventHeader();
-    
-    AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
-    AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(header->GenEventHeader());
-    Bool_t sd = kFALSE;
-    if(pythiaGenHeader) {
-      Int_t pythiaType = pythiaGenHeader->ProcessType();
-      if(pythiaType==92 || pythiaType==93)
-       sd = kTRUE;
-    }
-    if(dpmHeader) {
-      Int_t processType = dpmHeader->ProcessType();
-      if(processType == 5 || processType == 6)  
-       sd = kTRUE;
-      
-    }
-    if(!sd) {
-      triggers |= AliAODForwardMult::kMCNSD;
-      fHTriggers->Fill(kMCNSD+0.5);
-    }
     
-  }
-    
-  // Check if this is a collision candidate (INEL)
+  // Check if this is a collision candidate (MB)
   // Note, that we should use the value cached in the input 
   // handler rather than calling IsCollisionCandiate directly 
   // on the AliPhysicsSelection obejct.  If we called the latter
@@ -574,16 +562,6 @@ AliFMDEventInspector::ReadVertex(const AliESDEvent* esd, Double_t& vz)
     return kFALSE;
   }
 
-  // HHD test
-  /*
-  if (vertex->IsFromVertexerZ()) { 
-    return kFALSE;
-  }
-  if (TMath::Sqrt(TMath::Power(vertex->GetX(),2) + TMath::Power(vertex->GetY(),2)) > 3 ) { 
-      return kFALSE;
-  }
-  */
-  
   // Get the z coordiante 
   vz = vertex->GetZ();
   return kTRUE;
@@ -642,7 +620,7 @@ AliFMDEventInspector::Print(Option_t*) const
   field.ReplaceAll("m",  "-");
   field.ReplaceAll("kG", " kG");
   
-  std::cout << ind << "AliFMDEventInspector: " << GetName() << '\n'
+  std::cout << ind << ClassName() << ": " << GetName() << '\n'
            << ind << " Low flux cut:           " << fLowFluxCut << '\n'
            << ind << " Max(delta v_z):         " << fMaxVzErr << " cm\n"
            << ind << " System:                 " 
index 6cdb3aa..aaa3c8a 100644 (file)
@@ -111,7 +111,7 @@ public:
    * 
    * @param vtxAxis Vertex axis in use 
    */
-  void Init(const TAxis& vtxAxis);
+  virtual void Init(const TAxis& vtxAxis);
   /** 
    * Process the event 
    * 
@@ -226,6 +226,15 @@ protected:
    * @return @c true on success, @c false otherwise 
    */
   Bool_t ReadVertex(const AliESDEvent* esd, Double_t& vz);
+  /** 
+   * Read centrality from event 
+   * 
+   * @param esd  Event 
+   * @param cent On return, the centrality or negative if not found
+   * 
+   * @return False on error, true otherwise 
+   */
+  virtual Bool_t ReadCentrality(const AliESDEvent* esd, Double_t& cent);
 
   TH1I*    fHEventsTr;    //! Histogram of events w/trigger
   TH1I*    fHEventsTrVtx; //! Events w/trigger and vertex 
index 614dfe8..37fdb80 100644 (file)
@@ -61,7 +61,8 @@ AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
   fList           = o.fList;
   fSumRings       = o.fSumRings;
   fCoverage       = o.fCoverage;
-
+  fMergeMethod    = o.fMergeMethod;
+  fFiducialMethod = o.fFiducialMethod;
   return *this;
 }
 
@@ -93,8 +94,8 @@ AliFMDHistCollector::Init(const TAxis& vtxAxis,
   fList->Add(fSumRings);
 
   fCoverage = new TH2D("coverage", "#eta coverage per v_{z}", 
-                      etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
-                      vtxAxis.GetNbins(), vtxAxis.GetXmin(), vtxAxis.GetXmax());
+                      etaAxis.GetNbins(),etaAxis.GetXmin(),etaAxis.GetXmax(),
+                      vtxAxis.GetNbins(),vtxAxis.GetXmin(),vtxAxis.GetXmax());
   fCoverage->SetDirectory(0);
   fCoverage->SetXTitle("#eta");
   fCoverage->SetYTitle("v_{z} [cm]");
@@ -107,6 +108,14 @@ AliFMDHistCollector::Init(const TAxis& vtxAxis,
 
   // Find the eta bin ranges 
   for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
+    TList*   vtxList = new TList;
+    Double_t vMin    = vtxAxis.GetBinLowEdge(iVz);
+    Double_t vMax    = vtxAxis.GetBinUpEdge(iVz);
+    vtxList->SetName(Form("%c%02d_%c%02d", 
+                         vMin < 0 ? 'm' : 'p', int(TMath::Abs(vMin)),
+                         vMax < 0 ? 'm' : 'p', int(TMath::Abs(vMax))));
+    fList->Add(vtxList);
+
     // Find the first and last eta bin to use for each ring for 
     // each vertex bin.   This is instead of using the methods 
     // provided by AliFMDAnaParameters 
@@ -128,7 +137,7 @@ AliFMDHistCollector::Init(const TAxis& vtxAxis,
        // do not have holes in the coverage 
        bool ok = true;
        for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) { 
-         if (bg->GetBinContent(ie,ip) < fCorrectionCut) {
+         if (!CheckCorrection(bg, ie, ip)) {
            ok = false;
            continue;
          }
@@ -142,15 +151,54 @@ AliFMDHistCollector::Init(const TAxis& vtxAxis,
       // Store the result for later use 
       fFirstBins[(iVz-1)*5+iIdx] = first;
       fLastBins[(iVz-1)*5+iIdx]  = last;
-      for (Int_t iC = first+fNCutBins; iC <= last-fNCutBins; iC++) {
-       Double_t old = fCoverage->GetBinContent(iC, iVz);
-       fCoverage->SetBinContent(iC, iVz, old+1);
+      TH2D* obg = static_cast<TH2D*>(bg->Clone());
+      obg->SetDirectory(0);
+      obg->Reset();
+      vtxList->Add(obg);
+
+      // Fill diagnostics histograms 
+      for (Int_t ie = first+fNCutBins; ie <= last-fNCutBins; ie++) {
+       Double_t old = fCoverage->GetBinContent(ie, iVz);
+       fCoverage->SetBinContent(ie, iVz, old+1);
+       for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
+         obg->SetBinContent(ie, ip, bg->GetBinContent(ie, ip));
+         obg->SetBinError(ie, ip, bg->GetBinError(ie, ip));
+       }
       }
     } // for j 
   }
 }
 
 //____________________________________________________________________
+Bool_t
+AliFMDHistCollector::CheckCorrection(const TH2D* bg, Int_t ie, Int_t ip) const
+{
+  // 
+  // Check if we should include the bin in the data range 
+  // 
+  // Parameters:
+  //    bg Secondary map histogram
+  //    ie Eta bin
+  //    ip Phi bin
+  // 
+  // Return:
+  //    True if to be used
+  //
+  Double_t c = bg->GetBinContent(ie,ip);
+  switch (fFiducialMethod) { 
+  case kByCut:
+    return c >= fCorrectionCut;
+  case kDistance: 
+    if (2 * c < bg->GetBinContent(ie+1,ip) ||
+       2 * c < bg->GetBinContent(ie-1,ip)) return false;
+    return true;
+  default: 
+    AliError("No fiducal cut method defined");
+  }
+  return false;
+}
+    
+//____________________________________________________________________
 void
 AliFMDHistCollector::DefineOutput(TList* dir)
 {
@@ -342,6 +390,92 @@ AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, UShort_t vtxbin) const
   
   
 //____________________________________________________________________
+void
+AliFMDHistCollector::MergeBins(Double_t c,   Double_t e, 
+                              Double_t oc,  Double_t oe,
+                              Double_t& rc, Double_t& re) const
+{
+  // 
+  // Merge bins accoring to set method
+  // 
+  // Parameters:
+  //    c   Current content
+  //    e   Current error
+  //    oc  Old content
+  //    oe  Old error
+  //    rc  On return, the new content
+  //    re  On return, tne new error
+  //
+  rc = re = 0;
+  switch (fMergeMethod) { 
+  case kStraightMean:
+    // calculate the average of old value (half the original), 
+    // and this value, as well as the summed squared errors 
+    // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2) 
+    // and half the error of this.   
+    //
+    // So, on the first overlapping histogram we get 
+    // 
+    //    c = c_1 / 2
+    //    e = sqrt((e_1 / 2)^2) = e_1/2
+    // 
+    // On the second we get 
+    // 
+    //    c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2 
+    //    e' = sqrt(e^2 + (e_2/2)^2) 
+    //       = sqrt(e_1^2/4 + e_2^2/4) 
+    //       = sqrt(1/4 * (e_1^2+e_2^2)) 
+    //       = 1/2 * sqrt(e_1^2 + e_2^2)
+    rc = oc + c/2;
+    re = TMath::Sqrt(oe*oe+(e*e)/4);
+    break;
+  case kStraightMeanNoZero:
+    // If there's data in the overlapping histogram, 
+    // calculate the average and add the errors in 
+    // quadrature.  
+    // 
+    // If there's no data in the overlapping histogram, 
+    // then just fill in the data 
+    if (oe > 0) {
+      rc = (oc + c)/2;
+      re = TMath::Sqrt(oe*oe + e*e)/2;
+    }
+    else {
+      rc = c;
+      re = e;
+    }      
+    break;
+  case kWeightedMean: {
+    // Calculate the weighted mean 
+    Double_t w  = 1/(e*e);
+    Double_t sc = w * c;
+    Double_t sw = w;
+    if (oe > 0) {
+      Double_t ow =  1/(oe*oe);
+      sc          += ow * oc;
+      sw          += ow;
+    }
+    rc = sc / sw;
+    re = TMath::Sqrt(1 / sw);
+  }
+    break;
+  case kLeastError:
+    if (e < oe) {
+      rc = c;
+      re = e;
+    }
+    else {
+      rc = oc;
+      re = oe;
+    }
+    break;
+  default:
+    AliError("No method for defining content of overlapping bins defined");
+    return;
+  }
+}
+
+//____________________________________________________________________
 Bool_t
 AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
                             UShort_t                vtxbin, 
@@ -406,82 +540,10 @@ AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
          Double_t oc = out.GetBinContent(iEta,iPhi);
          Double_t oe = out.GetBinError(iEta,iPhi);
 
-#define USE_STRAIGHT_MEAN
-// #define USE_STRAIGHT_MEAN_NONZERO
-// #define USE_WEIGHTED_MEAN
-// #define USE_MOST_CERTAIN
-#if defined(USE_STRAIGHT_MEAN)
-         // calculate the average of old value (half the original), 
-         // and this value, as well as the summed squared errors 
-         // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2) 
-         // and half the error of this.   
-         //
-         // So, on the first overlapping histogram we get 
-         // 
-         //    c = c_1 / 2
-         //    e = sqrt((e_1 / 2)^2) = e_1/2
-         // 
-         // On the second we get 
-         // 
-         //    c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2 
-         //    e' = sqrt(e^2 + (e_2/2)^2) 
-         //       = sqrt(e_1^2/4 + e_2^2/4) 
-         //       = sqrt(1/4 * (e_1^2+e_2^2)) 
-         //       = 1/2 * sqrt(e_1^2 + e_2^2)
-         out.SetBinContent(iEta,iPhi,oc + c/2);
-         out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe+(e*e)/4));
-#elif defined(USE_STRAIGHT_MEAN_NONZERO)
-# define ZERO_OTHER
-         // If there's data in the overlapping histogram, 
-         // calculate the average and add the errors in 
-         // quadrature.  
-         // 
-         // If there's no data in the overlapping histogram, 
-         // then just fill in the data 
-         if (oe > 0) {
-           out.SetBinContent(iEta,iPhi,(oc + c)/2);
-           out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe + e*e)/2);
-         }
-         else {
-           out.SetBinContent(iEta,iPhi,c);
-           out.SetBinError(iEta,iPhi,e);
-         }         
-#elif defined(USE_WEIGHTED_MEAN) 
-         // Calculate the weighted mean 
-         Double_t w  = 1/(e*e);
-         Double_t sc = w * c;
-         Double_t sw = w;
-         if (oe > 0) {
-           Double_t ow = 1/(oe*oe);
-           sc          += ow * oc;
-           sw          += ow;
-         }
-         Double_t nc = sc / sw;
-         Double_t ne = TMath::Sqrt(1 / sw);
-         out.SetBinContent(iEta,iPhi,nc,ne);
-#elif defined(USE_MOST_CERTAIN)
-# define ZERO_OTHER
-         if (e < oe) {
-           out.SetBinContent(iEta,iPhi,c);
-           out.SetBinError(iEta,iPhi,e);
-         }
-         else {
-           out.SetBinContent(iEta,iPhi,oc);
-           out.SetBinError(iEta,iPhi,oe);
-         }
-#else 
-#         error No method for defining content of overlapping bins defined
-#endif
-#if defined(ZERO_OTHER)
-         // Get the content of the overlapping histogram, 
-         // and zero the content so that we won't use it 
-         // again 
-         UShort_t od; Char_t oR; 
-         GetDetRing(overlap, od, oR);
-         TH2D* other = hists.Get(od,oR);
-         other->SetBinContent(iEta,iPhi,0);
-         other->SetBinError(iEta,iPhi,0);
-#endif
+         Double_t rc, re;
+         MergeBins(c, e, oc, oe, rc, re);
+         out.SetBinContent(iEta,iPhi, rc);
+         out.SetBinError(iEta,iPhi, re);
        }
       }
       // Remove temporary histogram 
@@ -507,7 +569,15 @@ AliFMDHistCollector::Print(Option_t* /* option */) const
   std::cout << ind << "AliFMDHistCollector: " << GetName() << '\n'
            << ind << " # of cut bins:          " << fNCutBins << '\n'
            << ind << " Correction cut:         " << fCorrectionCut << '\n'
-           << ind << " Bin ranges:\n" << ind << "  v_z bin";
+           << ind << " Merge method:           ";
+  switch (fMergeMethod) {
+  case kStraightMean:       std::cout << "straight mean\n"; break;
+  case kStraightMeanNoZero: std::cout << "straight mean (no zeros)\n"; break;
+  case kWeightedMean:       std::cout << "weighted mean\n"; break;
+  case kLeastError:         std::cout << "least error\n"; break;
+  }
+    
+  std::cout << ind << " Bin ranges:\n" << ind << "  v_z bin";
   Int_t nVz = fFirstBins.fN / 5;
   for (UShort_t iVz = 1; iVz <= nVz; iVz++) 
     std::cout << " | " << std::setw(7) << iVz;
index e259bb8..34fe9b1 100644 (file)
@@ -31,6 +31,59 @@ class AliFMDHistCollector : public TNamed
 {
 public:
   /** 
+   * Methods to use when merging overlapping bins @f$b_1@f$, @f$b_2@f$
+   * with content @f$c_1@f$, @f$c_2@f$, and errors @f$e_1@f$,
+   * @f$e_2@f$ into bin @f$ b@f$ with content @f$c@f$ and error @f$e@f$ 
+   */
+  enum MergeMethod {
+    /**
+     * @f[
+     *   c = \frac{1}{2}(c_1+c_2) 
+     * @f]
+     * @f[
+     *   e = \sqrt{e_1^2+e_2^2} 
+     * @f]
+     */
+    kStraightMean,       
+    /**
+     * As above, exept zero's are ignored 
+     */
+    kStraightMeanNoZero, 
+    /** 
+     * @f[ 
+     *   c = \frac{\frac{c_1}{e_1^2}+\frac{c_2}{e_2^2}}{
+     *             \frac{1}{e_1^2}+\frac{1}{e_2^2}}
+     * @f]
+     * @f[
+     *   e = \sqrt{\frac{1}{\frac{1}{e_1^2}+\frac{1}{e_2^2}}
+     * @f]
+     */
+    kWeightedMean, 
+    /** 
+     * @f[
+     *     c = \left\{\begin{array}{cl}
+     *          c_1 & \text{if $e_1 < e_2} \\
+     *          c_2 & \text{otherwise}\end{array}\right.
+     * @f]
+     */
+    kLeastError
+  };
+  /**
+   * How to obtain the fiducial cuts 
+   */
+  enum FiducialMethod { 
+    /**
+     * Select bins by fixed cut.  Bins with a secondary correction
+     * less than the cut is considered as non-valid
+     */
+    kByCut, 
+    /**
+     * A bin is considered non-valid, if it is less then twice as
+     * large as it's neighbors (in eta)
+     */
+    kDistance 
+  };
+  /** 
    * Constructor 
    */
   AliFMDHistCollector() 
@@ -41,7 +94,9 @@ public:
       fDebug(0),
       fList(0),
       fSumRings(0),
-      fCoverage(0)
+      fCoverage(0),
+      fMergeMethod(kStraightMean),
+      fFiducialMethod(kByCut)
   {}
   /** 
    * Constructor 
@@ -57,7 +112,9 @@ public:
       fDebug(0),
       fList(0),
       fSumRings(0),
-      fCoverage(0)
+      fCoverage(0),
+      fMergeMethod(kStraightMean),
+      fFiducialMethod(kByCut)
   {}
   /** 
    * Copy constructor 
@@ -73,7 +130,9 @@ public:
       fDebug(o.fDebug),
       fList(o.fList),
       fSumRings(o.fSumRings),
-      fCoverage(o.fCoverage)
+      fCoverage(o.fCoverage),
+      fMergeMethod(o.fMergeMethod),
+      fFiducialMethod(o.fFiducialMethod)
   {}
 
   /** 
@@ -113,6 +172,18 @@ public:
    */  
   virtual void DefineOutput(TList* dir);
   /** 
+   * Set the merge method 
+   * 
+   * @param m Method
+   */
+  void SetMergeMethod(MergeMethod m) { fMergeMethod = m; }
+  /** 
+   * Set the method for finding the fidicual area of the secondary maps 
+   * 
+   * @param m Method
+   */
+  void SetFiducialMethod(FiducialMethod m) { fFiducialMethod = m; }
+  /** 
    * Set the number of extra bins (beyond the secondary map border) 
    * to cut away. 
    * 
@@ -262,7 +333,29 @@ protected:
    * @return True if there's an overlapping histogram 
    */
   Bool_t HasOverlap(Int_t i, Int_t e, UShort_t v) const;
-
+  /** 
+   * Check if we should include the bin in the data range 
+   * 
+   * @param bg Secondary map histogram
+   * @param ie Eta bin
+   * @param ip Phi bin
+   * 
+   * @return True if to be used
+   */
+  Bool_t CheckCorrection(const TH2D* bg, Int_t ie, Int_t ip) const;
+  /** 
+   * Merge bins accoring to set method
+   * 
+   * @param c   Current content
+   * @param e   Current error
+   * @param oc  Old content
+   * @param oe  Old error
+   * @param rc  On return, the new content
+   * @param re  On return, tne new error
+   */
+  void MergeBins(Double_t c,   Double_t e, 
+                Double_t oc,  Double_t oe,
+                Double_t& rc, Double_t& re) const;
 
   Int_t       fNCutBins;        // Number of additional bins to cut away
   Float_t     fCorrectionCut;   // Cut-off on secondary corrections 
@@ -272,7 +365,8 @@ protected:
   TList*      fList;           // Output list
   TH2D*       fSumRings;        // Sum per ring (on y-axis)
   TH2D*       fCoverage;        // Sum per ring (on y-axis)
-
+  MergeMethod fMergeMethod;     // Merge methiod for overlapping bins 
+  FiducialMethod fFiducialMethod; // Fidicual method
   ClassDef(AliFMDHistCollector,1); // Calculate Nch density 
 };
 
diff --git a/PWG2/FORWARD/analysis2/AliFMDMCEventInspector.cxx b/PWG2/FORWARD/analysis2/AliFMDMCEventInspector.cxx
new file mode 100644 (file)
index 0000000..74a9b01
--- /dev/null
@@ -0,0 +1,300 @@
+// 
+// This class inspects the event 
+//
+// Input:
+//   - AliESDFMD object possibly corrected for sharing
+//
+// Output:
+//   - A histogram of v_z of events with triggers. 
+//   - A histogram of v_z of events with vertex and triggers 
+//   - A histogram of trigger counters 
+// 
+// Note, that these are added to the master output list 
+//
+// Corrections used: 
+//   - None
+//
+#include "AliFMDMCEventInspector.h"
+#include "AliLog.h"
+#include "AliAODForwardMult.h"
+#include "AliForwardUtil.h"
+#include "AliCentrality.h"
+#include "AliESDEvent.h"
+#include "AliMCEvent.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliGenDPMjetEventHeader.h"
+#include "AliGenHijingEventHeader.h"
+// #include "AliGenHydjetEventHeader.h"
+#include "AliGenEposEventHeader.h"
+#include "AliGenHerwigEventHeader.h"
+#include "AliGenGeVSimEventHeader.h"
+#include "AliHeader.h"
+#include <TList.h>
+//====================================================================
+AliFMDMCEventInspector::AliFMDMCEventInspector()
+  : AliFMDEventInspector(), 
+    fHVertex(0),
+    fHPhiR(0), 
+    fHB(0)
+{
+  // 
+  // Constructor 
+  //
+}
+
+//____________________________________________________________________
+AliFMDMCEventInspector::AliFMDMCEventInspector(const char* /* name */)
+  : AliFMDEventInspector("fmdEventInspector"), 
+    fHVertex(0),
+    fHPhiR(0), 
+    fHB(0)
+{
+  // 
+  // Constructor 
+  // 
+  // Parameters:
+  //   name Name of object
+  //
+}
+
+//____________________________________________________________________
+AliFMDMCEventInspector::AliFMDMCEventInspector(const AliFMDMCEventInspector& o)
+  : AliFMDEventInspector(o), 
+    fHVertex(0),
+    fHPhiR(0), 
+    fHB(0)
+{
+  // 
+  // Copy constructor 
+  // 
+  // Parameters:
+  //   o Object to copy from 
+  //
+}
+
+//____________________________________________________________________
+AliFMDMCEventInspector::~AliFMDMCEventInspector()
+{
+  // 
+  // Destructor 
+  //
+}
+//____________________________________________________________________
+AliFMDMCEventInspector&
+AliFMDMCEventInspector::operator=(const AliFMDMCEventInspector& o)
+{
+  // 
+  // Assignement operator
+  // 
+  // Parameters:
+  //   o Object to assign from 
+  // 
+  // Return:
+  //    Reference to this object
+  //
+  AliFMDEventInspector::operator=(o);
+  return *this;
+}
+
+//____________________________________________________________________
+void
+AliFMDMCEventInspector::Init(const TAxis& vtxAxis)
+{
+  // 
+  // Initialize the object 
+  // 
+  // Parameters:
+  //   vtxAxis Vertex axis in use 
+  //
+  AliFMDEventInspector::Init(vtxAxis);
+
+  fHVertex = new TH1F("vertex", "True vertex distribution", 
+                     vtxAxis.GetNbins(), 
+                     vtxAxis.GetXmin(), 
+                     vtxAxis.GetXmax());
+  fHVertex->SetXTitle("v_{z} [cm]");
+  fHVertex->SetYTitle("# of events");
+  fHVertex->SetFillColor(kGreen+1);
+  fHVertex->SetFillStyle(3001);
+  fHVertex->SetDirectory(0);
+  // fHVertex->Sumw2();
+  fList->Add(fHVertex);
+
+  fHPhiR = new TH1F("phiR", "Event plane", 120, 0, 2*TMath::Pi());
+  fHPhiR->SetXTitle("#Phi_{R} [radians]");
+  fHPhiR->SetYTitle("# of events");
+  fHPhiR->SetFillColor(kGreen+1);
+  fHPhiR->SetFillStyle(3001);
+  fHPhiR->SetDirectory(0);
+  fList->Add(fHPhiR);
+
+  fHB = new TH1F("b", "Impact parameter", 125, 0, 25);
+  fHB->SetXTitle("b [fm]");
+  fHB->SetYTitle("# of events");
+  fHB->SetFillColor(kGreen+1);
+  fHB->SetFillStyle(3001);
+  fHB->SetDirectory(0);
+  fList->Add(fHB);
+  
+}
+
+//____________________________________________________________________
+UInt_t
+AliFMDMCEventInspector::ProcessMC(AliMCEvent*       event, 
+                                 UInt_t&           triggers,
+                                 UShort_t&         ivz, 
+                                 Double_t&         vz,
+                                 Double_t&         b,
+                                 Double_t&         phiR)
+{
+  // 
+  // Process the event 
+  // 
+  // Parameters:
+  //   event     Input event 
+  //   triggers  On return, the triggers fired 
+  //   ivz       On return, the found vertex bin (1-based).  A zero
+  //                  means outside of the defined vertex range
+  //   vz        On return, the z position of the interaction
+  //   b         On return, the impact parameter - < 0 if not available
+  //   phiR      On return, the event plane angle - < 0 if not available 
+  // 
+  // Return:
+  //    0 (or kOk) on success, otherwise a bit mask of error codes 
+  //
+
+  // Check that we have an event 
+  if (!event) { 
+    AliWarning("No MC event found for input event");
+    return kNoEvent;
+  }
+
+  //Assign MC only triggers : True NSD etc.
+  AliHeader*               header          = event->Header();
+  AliGenEventHeader*       genHeader       = header->GenEventHeader();
+  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);
+  // AliGenHydjetEventHeader* hydjetHeader    = 
+  //   dynamic_cast<AliGenHydjetEventHeader*>(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;
+  if(pythiaHeader) {
+    Int_t pythiaType = pythiaHeader->ProcessType();
+    if (pythiaType==92 || pythiaType==93) sd = kTRUE;
+    b = pythiaHeader->GetImpactParameter();
+  }
+  if(dpmHeader) {
+    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;
+  }
+  // if (hydjetHeader) {
+  //   b    = hydjetHeader->ImpactParameter();
+  //   phi  = hydjetHeader->ReactionPlaneAngle();
+  // }
+  if (eposHeader) {
+    b    = eposHeader->ImpactParameter();
+    phi  = eposHeader->ReactionPlaneAngle();
+  }
+
+  // Normalize event plane angle to [0,2pi]
+  if (phi <= -1111) phiR = -1;
+  else { 
+    while (true) {
+      if      (phi < 0)             phi += 2*TMath::Pi();
+      else if (phi > 2*TMath::Pi()) phi -= 2*TMath::Pi();
+      else                          break;
+    }
+  }
+
+  // Set NSD flag
+  if(!sd) {
+    triggers |= AliAODForwardMult::kMCNSD;
+    fHTriggers->Fill(kMCNSD+0.5);
+  }
+  
+  // Get the primary vertex from EG 
+  TArrayF vtx;
+  genHeader->PrimaryVertex(vtx);
+  vz = vtx[2];
+
+  fHVertex->Fill(vz);
+  fHPhiR->Fill(phiR);
+  fHB->Fill(b);
+
+  // Check for the vertex bin 
+  ivz = fHEventsTr->GetXaxis()->FindBin(vz);
+  if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) { 
+    if (fDebug > 3) {
+      AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
+                     vz, fHEventsTr->GetXaxis()->GetXmin(), 
+                     fHEventsTr->GetXaxis()->GetXmax())); }
+    ivz = 0;
+    return kBadVertex;
+  }
+
+  
+  return kOk;
+}
+//____________________________________________________________________
+Bool_t
+AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd, Double_t& cent)
+{
+  // 
+  // Read centrality from event 
+  // 
+  // Parameters:
+  //    esd  Event 
+  //    cent On return, the centrality or negative if not found
+  // 
+  // Return:
+  //    False on error, true otherwise 
+  //
+  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);
+
+  return true;
+}
+
+  
+//
+// EOF
+//
+
diff --git a/PWG2/FORWARD/analysis2/AliFMDMCEventInspector.h b/PWG2/FORWARD/analysis2/AliFMDMCEventInspector.h
new file mode 100644 (file)
index 0000000..50812f3
--- /dev/null
@@ -0,0 +1,110 @@
+// 
+// This class inspects the event 
+//
+#ifndef ALIFMDMCEVENTINSPECTOR_H
+#define ALIFMDMCEVENTINSPECTOR_H
+#include "AliFMDEventInspector.h"
+class AliMCEvent;
+
+/** 
+ * This class inspects the event 
+ *
+ * @par Input:
+ *   - AliESDFMD object possibly corrected for sharing
+ *
+ * @par Output:
+ *   - A histogram of v_z of events with triggers. 
+ *   - A histogram of v_z of events with vertex and triggers 
+ *   - A histogram of trigger counters 
+ * 
+ * Note, that these are added to the master output list 
+ *
+ * @par Corrections used: 
+ *   - None
+ *
+ * @ingroup pwg2_forward_algo 
+ */
+class AliFMDMCEventInspector : public AliFMDEventInspector
+{
+public:
+  /** 
+   * Constructor 
+   */
+  AliFMDMCEventInspector();
+  /** 
+   * Constructor 
+   * 
+   * @param name Name of object
+   */
+  AliFMDMCEventInspector(const char* name);
+  /** 
+   * Copy constructor 
+   * 
+   * @param o Object to copy from 
+   */
+  AliFMDMCEventInspector(const AliFMDMCEventInspector& o);
+  /** 
+   * Destructor 
+   */
+  virtual ~AliFMDMCEventInspector();
+  /** 
+   * Assignement operator
+   * 
+   * @param o Object to assign from 
+   * 
+   * @return Reference to this object
+   */
+  AliFMDMCEventInspector& operator=(const AliFMDMCEventInspector&);
+
+  /** 
+   * Initialize the object 
+   * 
+   * @param vtxAxis Vertex axis in use 
+   */
+  void Init(const TAxis& vtxAxis);
+  /** 
+   * Process MC truth event.  Note, returned values are the MC truth
+   * values
+   * 
+   * @param event     Input event 
+   * @param triggers  On return, the triggers fired 
+   * @param lowFlux   On return, true if the event is considered a low-flux 
+   *                  event (according to the setting of fLowFluxCut) 
+   * @param ivz       On return, the found vertex bin (1-based).  A zero
+   *                  means outside of the defined vertex range
+   * @param vz        On return, the z position of the interaction
+   * @param cent      On return, the centrality (in percent) or < 0 
+   *                  if not found
+   * 
+   * @return 0 (or kOk) on success, otherwise a bit mask of error codes 
+   */
+  UInt_t ProcessMC(AliMCEvent*       event, 
+                  UInt_t&           triggers,
+                  UShort_t&         ivz, 
+                  Double_t&         vz,
+                  Double_t&         b,
+                  Double_t&         phiR);
+protected:
+  /** 
+   * Read centrality from event 
+   * 
+   * @param esd  Event 
+   * @param cent On return, the centrality or negative if not found
+   * 
+   * @return False on error, true otherwise 
+   */
+  virtual Bool_t ReadCentrality(const AliESDEvent* esd, Double_t& cent);
+
+  TH1F* fHVertex; // Histogram of vertex 
+  TH1F* fHPhiR;   // Histogram of event plane 
+  TH1F* fHB;      // Histogram of impact parameter 
+  ClassDef(AliFMDMCEventInspector,1); // Inspect the event 
+};
+
+#endif
+// Local Variables:
+//   mode: C++
+// End:
+
+
+
index 26b8f9a..ec5a90a 100644 (file)
@@ -199,6 +199,7 @@ AliForwardMCMultiplicityTask::UserCreateOutputObjects()
   // 
   //
   fList = new TList;
+  fList->SetOwner();
 
   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
   AliAODHandler*      ah = 
@@ -243,7 +244,8 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   //  
 
   // Get the input data 
-  AliESDEvent* esd = GetESDEvent();
+  AliESDEvent* esd     = GetESDEvent();
+  AliMCEvent*  mcEvent = MCEvent();
 
   // Clear stuff 
   fHistos.Clear();
@@ -261,6 +263,12 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   Double_t cent     = 0;
   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, 
                                              ivz, vz, cent);
+  UShort_t ivzMC    = 0;
+  Double_t vzMC     = 0;
+  Double_t phiR     = 0;
+  Double_t b        = 0;
+  // UInt_t   foundMC  = 
+  fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, phiR);
   
   
   //Store all events
@@ -278,7 +286,8 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   //All events should be stored - HHD
   //if (isAccepted) MarkEventForStore();
 
-  if (found & AliFMDEventInspector::kNoSPD)     isAccepted = false; // return;
+  // Disable this check on SPD - will bias data 
+  // if (found & AliFMDEventInspector::kNoSPD)  isAccepted = false; // return;
   if (found & AliFMDEventInspector::kNoFMD)     isAccepted = false; // return;
   if (found & AliFMDEventInspector::kNoVertex)  isAccepted = false; // return;
 
@@ -295,7 +304,7 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
 
   // Get FMD data 
   AliESDFMD*  esdFMD  = esd->GetFMDData();
-  AliMCEvent* mcEvent = MCEvent();
+
   // Apply the sharing filter (or hit merging or clustering if you like)
   if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
     AliWarning("Sharing filter failed!");
index f7eb539..e208497 100644 (file)
@@ -5,7 +5,7 @@
 #define ALIFORWARDMCMULTIPLICITYTASK_H
 #include "AliForwardMultiplicityBase.h"
 #include "AliForwardUtil.h"
-#include "AliFMDEventInspector.h"
+#include "AliFMDMCEventInspector.h"
 #include "AliFMDEnergyFitter.h"
 #include "AliFMDMCSharingFilter.h"
 #include "AliFMDMCDensityCalculator.h"
@@ -190,7 +190,7 @@ protected:
   AliAODForwardMult      fMCAODFMD;     // MC Output object
   TH2D*                  fPrimary;      // Per event primary particles 
 
-  AliFMDEventInspector      fEventInspector;    // Algorithm
+  AliFMDMCEventInspector    fEventInspector;    // Algorithm
   AliFMDEnergyFitter        fEnergyFitter;      // Algorithm
   AliFMDMCSharingFilter     fSharingFilter;     // Algorithm
   AliFMDMCDensityCalculator fDensityCalculator; // Algorithm
index b775e57..9c768ad 100644 (file)
 
 //____________________________________________________________________
 AliForwarddNdetaTask::AliForwarddNdetaTask()
-  : AliBasedNdetaTask(),
-    fSumPrimary(0),    //  Sum of primary histograms
-    fSNNString(0),
-    fSysString(0)
+  : AliBasedNdetaTask()
 {
   //
   // Constructor 
@@ -26,10 +23,7 @@ AliForwarddNdetaTask::AliForwarddNdetaTask()
 
 //____________________________________________________________________
 AliForwarddNdetaTask::AliForwarddNdetaTask(const char* /* name */)
-  : AliBasedNdetaTask("Forward"), 
-    fSumPrimary(0),    //  Sum of primary histograms
-    fSNNString(0),
-    fSysString(0)
+  : AliBasedNdetaTask("Forward")
 {
   // 
   // Constructor
@@ -40,10 +34,7 @@ AliForwarddNdetaTask::AliForwarddNdetaTask(const char* /* name */)
 
 //____________________________________________________________________
 AliForwarddNdetaTask::AliForwarddNdetaTask(const AliForwarddNdetaTask& o)
-  : AliBasedNdetaTask(o),
-    fSumPrimary(o.fSumPrimary),           // Sum of primary histograms
-    fSNNString(o.fSNNString),     //  
-    fSysString(o.fSysString)      //  
+  : AliBasedNdetaTask(o)
 {
   // 
   // Copy constructor
@@ -51,6 +42,65 @@ AliForwarddNdetaTask::AliForwarddNdetaTask(const AliForwarddNdetaTask& o)
 }
 
 //____________________________________________________________________
+AliBasedNdetaTask::CentralityBin*
+AliForwarddNdetaTask::MakeCentralityBin(const char* name, Short_t l, Short_t h) 
+  const 
+{
+  // 
+  // Make a new centrality bin
+  // 
+  // Parameters:
+  //    name   Histogram names
+  //    l      Lower cut
+  //    h      Upper cut
+  // 
+  // Return:
+  //    Newly allocated object (of our type)
+  //
+  return new AliForwarddNdetaTask::CentralityBin(name, l, h);
+}
+
+//____________________________________________________________________
+void
+AliForwarddNdetaTask::UserExec(Option_t* option)
+{
+  // 
+  // Called at each event 
+  //
+  // This is called once in the master 
+  // 
+  // Parameters:
+  //    option Not used 
+  //
+  AliBasedNdetaTask::UserExec(option);
+  
+  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
+  if (!aod) {
+    AliError("Cannot get the AOD event");
+    return;
+  } 
+
+  TObject* obj = aod->FindListObject("Forward");
+  if (!obj) { 
+    AliWarning("No forward object found");
+    return;
+  }
+  AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
+  TObject* oPrimary   = aod->FindListObject("primary");
+  if (!oPrimary) return;
+  
+  TH2D* primary   = static_cast<TH2D*>(oPrimary);
+
+  // Loop over centrality bins 
+  TIter next(fListOfCentralities);
+  CentralityBin* bin = 0;
+  while ((bin = static_cast<CentralityBin*>(next()))) 
+    bin->ProcessPrimary(forward, fTriggerMask, fVtxMin, fVtxMax, primary);
+}  
+  
+
+//____________________________________________________________________
 TH2D*
 AliForwarddNdetaTask::GetHistogram(const AliAODEvent* aod, Bool_t mc)
 {
@@ -74,110 +124,82 @@ AliForwarddNdetaTask::GetHistogram(const AliAODEvent* aod, Bool_t mc)
     return 0;
   }
   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
+  return &(forward->GetHistogram());
+}
 
-  // Here, we update the primary stuff 
-  if (mc) { 
-    TObject* oPrimary   = aod->FindListObject("primary");
-    if (oPrimary) {
-      TH2D* primary   = static_cast<TH2D*>(oPrimary);
-      // Add contribtion from MC 
-      
-      if(fTriggerMask == AliAODForwardMult::kInel) {
-       if (primary && !fSumPrimary) fSumPrimary = CloneHist(primary, "truth");
-       else if (primary) fSumPrimary->Add(primary);
-      }
-      
-      else if(fTriggerMask == AliAODForwardMult::kNSD && forward->IsTriggerBits(AliAODForwardMult::kMCNSD)) {
-       if (primary && !fSumPrimary) fSumPrimary = CloneHist(primary, "truth"); 
-       else if (primary) fSumPrimary->Add(primary);
-      }
-      
-      
-      
-      /*   if(fTriggerMask == AliAODForwardMult::kNSD && forward->IsTriggerBits(AliAODForwardMult::kMCNSD)) {
-          if (primary) fSumPrimary->Add(primary); 
-          
-          }
-          
-          else if (primary) fSumPrimary->Add(primary);*/
-      
-    }    
+//========================================================================
+void
+AliForwarddNdetaTask::CentralityBin::ProcessPrimary(const AliAODForwardMult* 
+                                                   forward, 
+                                                   Int_t triggerMask,
+                                                   Double_t vzMin, 
+                                                   Double_t vzMax, 
+                                                   const TH2D* primary)
+{ 
+  // Check the centrality class unless this is the 'all' bin 
+  if (!IsAllBin()) { 
+    Double_t centrality = forward->GetCentrality();
+    if (centrality < fLow || centrality >= fHigh) return;
   }
+  // Check if we have an event of interest. 
+  if (!forward->IsTriggerBits(triggerMask)) return;
   
-  // Here, we get the update 
-  if (!fSNNString) { 
-    UShort_t sNN = forward->GetSNN();
-    fSNNString = new TNamed("sNN", "");
-    fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
-    fSNNString->SetUniqueID(sNN);
-    fSums->Add(fSNNString);
+  //Check for pileup
+  if (forward->IsTriggerBits(AliAODForwardMult::kPileUp)) return;
   
-    UShort_t sys = forward->GetSystem();
-    fSysString = new TNamed("sys", "");
-    fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
-    fSysString->SetUniqueID(sys);
-    fSums->Add(fSysString);
+  // Check that we have a valid vertex
+  if (!forward->HasIpZ()) return;
+
+  // Check that vertex is within cuts 
+  if (!forward->InRange(vzMin, vzMax)) return;
+
+  if (!fSumPrimary) { 
+    fSumPrimary = static_cast<TH2D*>(primary->Clone("truth"));
+    fSumPrimary->SetDirectory(0);
+    fSumPrimary->Reset();
+    fSums->Add(fSumPrimary);
   }
 
-  return &(forward->GetHistogram());
+  if (triggerMask == AliAODForwardMult::kInel ||
+      (triggerMask == AliAODForwardMult::kNSD && 
+       forward->IsTriggerBits(AliAODForwardMult::kMCNSD)))
+    fSumPrimary->Add(primary);
 }
 
 //________________________________________________________________________
-void 
-AliForwarddNdetaTask::Terminate(Option_t *) 
+void
+AliForwarddNdetaTask::CentralityBin::End(TList*      sums, 
+                                        TList*      results,
+                                        const TH1*  shapeCorr, 
+                                        Double_t    trigEff,
+                                        Bool_t      symmetrice,
+                                        Int_t       rebin, 
+                                        Bool_t      corrEmpty, 
+                                        Bool_t      cutEdges, 
+                                        Double_t    vzMin, 
+                                        Double_t    vzMax, 
+                                        Int_t       triggerMask)
 {
-  // Draw result to screen, or perform fitting, normalizations Called
-  // once at the end of the query
-  AliBasedNdetaTask::Terminate("");
-  if(!fSums) {
-    AliError("Could not retrieve TList fSums"); 
-    return; 
-  }
-  if (!fOutput) { 
-    AliError("Could not retrieve TList fOutput"); 
-    return; 
-  }
+  AliBasedNdetaTask::CentralityBin::End(sums, results, shapeCorr, trigEff, 
+                                       symmetrice, rebin, corrEmpty, cutEdges,
+                                       vzMin, vzMax, triggerMask);
 
   fSumPrimary     = static_cast<TH2D*>(fSums->FindObject("truth"));
-  
-  fOutput->Add(fSNNString->Clone());
-  fOutput->Add(fSysString->Clone());
 
-  if (fSumPrimary) { 
-    Int_t nAll        = Int_t(fTriggers->GetBinContent(kAll));
-    Int_t nNSD        = Int_t(fTriggers->GetBinContent(kMCNSD));
-    
-    //  for(Int_t nn =1; nn<=fSumPrimary->GetNbinsX(); nn++) {
-    //  for(Int_t mm =1; mm<=fSumPrimary->GetNbinsY(); mm++) {
-    // if(fSumPrimary->GetBinContent(nn,mm) > 0) std::cout<<fSumPrimary->GetBinContent(nn,mm)<<" +/-  "<<fSumPrimary->GetBinError(nn,mm)<<std::endl;
-    //  }
-    // }
-    
-    //TH1D* dndetaTruth = ProjectX(fSumPrimary,"dndetaTruth",
-    //                  1,fSumPrimary->GetNbinsY(),false,true);
-    TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",1,fSumPrimary->GetNbinsY(),"e");
+  if (!fSumPrimary) return;
+  Int_t n           = (triggerMask == AliAODForwardMult::kNSD ? 
+                      Int_t(fTriggers->GetBinContent(kMCNSD)) : 
+                      Int_t(fTriggers->GetBinContent(kAll)));
 
     
-    if(fTriggerMask == AliAODForwardMult::kNSD)
-      dndetaTruth->Scale(1./nNSD, "width");
-    else
-      dndetaTruth->Scale(1./nAll, "width");
-    
-    //for(Int_t nn =1; nn<=dndetaTruth->GetNbinsX(); nn++) {
-    //  if(dndetaTruth->GetBinContent(nn) > 0) std::cout<<dndetaTruth->GetBinContent(nn)<<" +/-  "<<dndetaTruth->GetBinError(nn)<<std::endl;
-    // }
+  TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",1,
+                                              fSumPrimary->GetNbinsY(),"e");
+  dndetaTruth->Scale(1./n, "width");
 
+  SetHistogramAttributes(dndetaTruth, kBlue+3, 22, "Monte-Carlo truth");
     
-    SetHistogramAttributes(dndetaTruth, kBlue+3, 22, "Monte-Carlo truth");
-    
-    fOutput->Add(dndetaTruth);
-    fOutput->Add(Rebin(dndetaTruth));
-  }
-  // If only there was a away to get sqrt{s_NN} and beam type 
-  
-
-  PostData(2, fOutput);
+  fOutput->Add(dndetaTruth);
+  fOutput->Add(Rebin(dndetaTruth, rebin, cutEdges));
 }
 
 //________________________________________________________________________
index ff12e59..74fd328 100644 (file)
@@ -32,13 +32,13 @@ public:
    */
   virtual ~AliForwarddNdetaTask() {}
   /** 
-   * Called at end of event processing.. 
+   * Called at each event 
    *
    * This is called once in the master 
    * 
    * @param option Not used 
    */
-  virtual void Terminate(Option_t* option);
+  virtual void UserExec(Option_t* option);
 protected:
   /** 
    * Copy constructor 
@@ -60,9 +60,98 @@ protected:
    * @return Retrieved histogram or null
    */
   TH2D* GetHistogram(const AliAODEvent* aod, Bool_t mc);
-  TH2D*           fSumPrimary;    //  Sum of primary histograms
-  TNamed*         fSNNString;     //  sqrt(s_NN) string 
-  TNamed*         fSysString;     //  Collision system string 
+  /** 
+   * Make a new centrality bin
+   * 
+   * @param name   Histogram names
+   * @param l      Lower cut
+   * @param h      Upper cut
+   * 
+   * @return Newly allocated object (of our type)
+   */
+  AliBasedNdetaTask::CentralityBin* 
+  MakeCentralityBin(const char* name, Short_t l, Short_t h) const;
+
+  struct CentralityBin : public AliBasedNdetaTask::CentralityBin 
+  {
+    /** 
+     * Constructor 
+     */
+    CentralityBin() : AliBasedNdetaTask::CentralityBin(), fSumPrimary(0) {}
+    /** 
+     * Constructor 
+     * 
+     * @param name Name used for histograms (e.g., Forward)
+     * @param low  Lower centrality cut in percent 
+     * @param high Upper centrality cut in percent 
+     */
+    CentralityBin(const char* name, Short_t low, Short_t high)
+      : AliBasedNdetaTask::CentralityBin(name, low, high), 
+       fSumPrimary(0)
+    {}
+    /** 
+     * Copy constructor 
+     * 
+     * @param other Object to copy from 
+     */
+    CentralityBin(const CentralityBin& other) 
+      : AliBasedNdetaTask::CentralityBin(other), 
+       fSumPrimary(other.fSumPrimary) 
+    {}
+    /** 
+     * Destructor 
+     */
+    virtual ~CentralityBin() {}
+    /** 
+     * Assignement operator 
+     * 
+     * 
+     * @return 
+     */
+    CentralityBin& operator=(const CentralityBin&) { return *this; }
+    /** 
+     * Process an event
+     * 
+     * @param forward     Forward data (for trigger, vertex, & centrality)
+     * @param triggerMask Trigger mask 
+     * @param vzMin       Minimum IP z coordinate
+     * @param vzMax       Maximum IP z coordinate
+     * @param primary     MC truth histogram
+     */
+    virtual void ProcessPrimary(const AliAODForwardMult* forward, 
+                               Int_t triggerMask,
+                               Double_t vzMin, Double_t vzMax, 
+                               const TH2D* primary);
+    /** 
+     * End of processing 
+     * 
+     * @param sums        List of sums
+     * @param results     Output list of results
+     * @param shapeCorr   Shape correction or nil
+     * @param trigEff     Trigger efficiency 
+     * @param symmetrice  Whether to symmetrice the results
+     * @param rebin       Whether to rebin the results
+     * @param corrEmpty   Whether to correct for empty bins
+     * @param cutEdges    Whether to cut edges when rebinning
+     * @param vzMin       Minimum IP z coordinate
+     * @param vzMax      Maximum IP z coordinate
+     * @param triggerMask Trigger mask 
+     */
+    virtual void End(TList*      sums, 
+                    TList*      results,
+                    const TH1*  shapeCorr, 
+                    Double_t    trigEff,
+                    Bool_t      symmetrice,
+                    Int_t       rebin, 
+                    Bool_t      corrEmpty, 
+                    Bool_t      cutEdges, 
+                    Double_t    vzMin, 
+                    Double_t    vzMax, 
+                    Int_t       triggerMask);
+  protected: 
+    TH2D*           fSumPrimary;    //  Sum of primary histograms
+    ClassDef(CentralityBin,1); // A centrality bin     
+  };
 
   ClassDef(AliForwarddNdetaTask,1); // Determine multiplicity in forward region
 };
index 91cedf7..937a566 100644 (file)
@@ -1,3 +1,8 @@
+/**
+ * Script to visualise the dN/deta 
+ *
+ * This script is independent of any AliROOT code - and should stay that way. 
+ */
 #include <TH1.h>
 #include <THStack.h>
 #include <TGraphErrors.h>
@@ -48,19 +53,18 @@ struct dNdetaDrawer
       fRebin(5),               // UShort_t
       fCutEdges(false),                // Bool_t
       fTitle(""),              // TString
-      fHHDFile(""),            // TString
       fTrigString(0),          // TNamed*
       fSNNString(0),           // TNamed*
       fSysString(0),           // TNamed*
       fVtxAxis(0),             // TAxis*
-      fForward(0),             // TH1*
-      fForwardMC(0),           // TH1*
-      fForwardHHD(0),          // TH1*
-      fTruth(0),               // TH1*
-      fCentral(0),             // TH1*
-      fForwardSym(0),          // TH1*
-      fForwardMCSym(0),                // TH1*
-      fForwardHHDSym(0),       // TH1*
+      // fForward(0),          // TH1*
+      // fForwardMC(0),                // TH1*
+      // fForwardHHD(0),       // TH1*
+      // fTruth(0),            // TH1*
+      // fCentral(0),          // TH1*
+      // fForwardSym(0),       // TH1*
+      // fForwardMCSym(0),     // TH1*
+      // fForwardHHDSym(0),    // TH1*
       fTriggers(0),            // TH1*
       fRangeParam(0)
 
@@ -127,13 +131,6 @@ struct dNdetaDrawer
    * @param x Title
    */
   void SetTitle(TString x)        { fTitle = x; }
-  //__________________________________________________________________
-  /** 
-   * Set the file name of the file containing the HHD results
-   * 
-   * @param fn File name 
-   */
-  void SetHHDFile(const char* fn) { fHHDFile = fn; }
   /* @} */
   //==================================================================  
   /** 
@@ -200,12 +197,13 @@ struct dNdetaDrawer
    */
   void Run(const char* filename="forward_dndeta.C") 
   {
-    if (!Open(filename)) return;
-
     Double_t max = 0;
 
+    if (!Open(filename, max)) return;
+    Info("Run", "Got data from %s with maximum %f", filename, max);
+
     // Create our stack of results
-    THStack* results = StackResults(max);
+    THStack* results = fResults; // StackResults(max);
 
     // Create our stack of other results 
     TMultiGraph* other = 0;
@@ -221,41 +219,15 @@ struct dNdetaDrawer
 
     Plot(results, other, max, ratios, smax, leftright, amax);
   }
-    
+
   //__________________________________________________________________
   /** 
-   * Open input file, and find data 
+   * Fetch the information on the run from the results list
    * 
-   * @param filename File name
-   * 
-   * @return true on success 
+   * @param results  Results list
    */
-  Bool_t Open(const char* filename)
+  void FetchInformation(const TList* results)
   {
-    TFile* file = TFile::Open(filename, "READ");
-    if (!file) { 
-      Error("Open", "Cannot open %s", filename);
-      return false;
-    }
-    
-    TList* results = static_cast<TList*>(file->Get("ForwardResults"));
-    if (!results) { 
-      Error("Open", "Couldn't find list ForwardResults");
-      return false;
-    }
-
-    fForward   = GetResult(results, "dndetaForward");
-    fForwardMC = GetResult(results, "dndetaForwardMC");
-    fTruth     = GetResult(results, "dndetaTruth");
-    if (!fTruth) results->ls();
-
-    TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
-    if (!clusters) 
-      Warning("Open", "Couldn't find list CentralResults");
-    else {
-      fCentral   = GetResult(clusters, "dndetaCentral");
-      if (fCentral) fCentral->SetMarkerColor(kGreen+1);
-    }
     if (!fTrigString) 
       fTrigString = static_cast<TNamed*>(results->FindObject("trigString"));
     if (!fSNNString) 
@@ -264,7 +236,6 @@ struct dNdetaDrawer
       fSysString  = static_cast<TNamed*>(results->FindObject("sys"));
     if (!fVtxAxis)
       fVtxAxis    = static_cast<TAxis*>(results->FindObject("vtxAxis"));
-    
     if (!fTrigString) fTrigString = new TNamed("trigString", "unknown");
     if (!fSNNString)  fSNNString  = new TNamed("sNN", "unknown");
     if (!fSysString)  fSysString  = new TNamed("sys", "unknown");
@@ -274,7 +245,7 @@ struct dNdetaDrawer
       fVtxAxis->SetTitle("v_{z} range unspecified");
     }
 
-    Info("Open", 
+    Info("FetchInformation", 
         "Initialized for\n"
         "   Trigger:    %s  (%d)\n"
         "   sqrt(sNN):  %s  (%dGeV)\n"
@@ -284,39 +255,69 @@ struct dNdetaDrawer
         fSNNString->GetTitle(),  fSNNString->GetUniqueID(), 
         fSysString->GetTitle(),  fSysString->GetUniqueID(), 
         fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
-
-    TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
-    if (sums) 
-      fTriggers = GetResult(sums, "triggers");
-
-    if (!fForward) { 
-      Error("Open", "Couldn't find the result of the forward analysis");
+  }
+    
+  //__________________________________________________________________
+  /** 
+   * Open input file, and find data 
+   * 
+   * @param filename File name
+   * 
+   * @return true on success 
+   */
+  Bool_t Open(const char* filename, Double_t& max)
+  {
+    TFile* file = TFile::Open(filename, "READ");
+    if (!file) { 
+      Error("Open", "Cannot open %s", filename);
       return false;
     }
-    file->Close();
+    
+    TList* results = static_cast<TList*>(file->Get("ForwardResults"));
+    if (!results) { 
+      Error("Open", "Couldn't find list ForwardResults");
+      return false;
+    }
+
+    // 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");
     
-    fForwardHHD = GetHHD();
+    file->Close();
 
     return true;
   }
   //__________________________________________________________________
-  /** 
-   * Make a histogram stack of results 
-   * 
-   * @param max On return, the maximum value in the stack 
-   * 
-   * @return Newly allocated stack
-   */
-  THStack* StackResults(Double_t& max)
+  void FetchResults(THStack* stack, const TList* list, const char* name, 
+                   Double_t& max)
   {
-    THStack* stack = new THStack("results", "Stack of Results");
-    max = TMath::Max(max, AddHistogram(stack, fTruth,      "e5 p"));
-    max = TMath::Max(max, AddHistogram(stack, fForwardHHD, "", fForwardHHDSym));
-    max = TMath::Max(max, AddHistogram(stack, fForwardMC,  "", fForwardMCSym));
-    max = TMath::Max(max, AddHistogram(stack, fCentral,    ""));
-    max = TMath::Max(max, AddHistogram(stack, fForward,    "", fForwardSym));
-    return stack;
+    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;
+
+    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));
+    
+    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);
   }
   //__________________________________________________________________
   /** 
@@ -326,7 +327,7 @@ struct dNdetaDrawer
    * 
    * @return Newly allocated stack
    */
-  TMultiGraph* StackOther(Double_t& max) const
+  TMultiGraph* StackOther(Double_t& max)
   {
     gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/OtherData.C");
     Int_t    error = 0;
@@ -349,8 +350,10 @@ struct dNdetaDrawer
 
     TGraphAsymmErrors* o      = 0;
     TIter              next(other->GetListOfGraphs());
-    while ((o = static_cast<TGraphAsymmErrors*>(next()))) 
+    while ((o = static_cast<TGraphAsymmErrors*>(next()))) {
       max = TMath::Max(max, TMath::MaxElement(o->GetN(), o->GetY()));
+      fDenominators.Add(o);
+    }
 
     return other;
   }
@@ -367,49 +370,31 @@ struct dNdetaDrawer
     THStack* ratios = new THStack("ratios", "Ratios");
 
     if (others) {
-      TGraphAsymmErrors* ua5_1  = 0;
-      TGraphAsymmErrors* ua5_2  = 0;
-      TGraphAsymmErrors* alice  = 0;
-      TGraphAsymmErrors* cms    = 0;
-      TGraphAsymmErrors* o      = 0;
-      TIter              nextG(others->GetListOfGraphs());
-      while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
-       ratios->Add(Ratio(fForward,          o, max));
-       ratios->Add(Ratio(fForwardSym,       o, max));
-       ratios->Add(Ratio(fForwardHHD,       o, max));
-       ratios->Add(Ratio(fForwardHHDSym,    o, max));
-       ratios->Add(Ratio(fCentral,          o, max));
-       TString oName(o->GetName());
+      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 = o; else ua5_1 = o; }
-       if (oName.Contains("alice")) alice = o;
-       if (oName.Contains("cms"))   cms = o;
+       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 (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));
     }
-
-    // Check if we have a primaries from MC 
-    if (fTruth) {
-      ratios->Add(Ratio(fForward,    fTruth, max));
-      ratios->Add(Ratio(fForwardSym, fTruth, max));
-      ratios->Add(Ratio(fCentral,    fTruth, max));
-    }
-
-    // If we have data from HHD's analysis, then do the ratio of 
-    // our result to that data. 
-    if (fForwardHHD) { 
-      ratios->Add(Ratio(fForward,    fForwardHHD,    max));
-      ratios->Add(Ratio(fForwardSym, fForwardHHDSym, max));
-    }
-
-    // Do comparison to MC 
-    if (fForwardMC) { 
-      ratios->Add(Ratio(fForward,    fForwardMC,    max));
-      ratios->Add(Ratio(fForwardSym, fForwardMCSym, max));
-    }
-
     // Check if we have ratios 
     if (!ratios->GetHists() || 
        (ratios->GetHists()->GetEntries() <= 0)) { 
@@ -429,9 +414,8 @@ struct dNdetaDrawer
   THStack* StackLeftRight(Double_t& max)
   {
     THStack* ret = new THStack("leftright", "Left-right asymmetry");
-    ret->Add(Asymmetry(fForward,    max));
-    ret->Add(Asymmetry(fForwardHHD, max));
-    ret->Add(Asymmetry(fForwardMC,  max));
+    // ret->Add(Asymmetry(fForward,    max));
+    // ret->Add(Asymmetry(fForwardMC,  max));
 
     if (!ret->GetHists() || 
        (ret->GetHists()->GetEntries() <= 0)) { 
@@ -538,6 +522,8 @@ 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);
 
@@ -674,8 +660,8 @@ struct dNdetaDrawer
 
     // Make a nice band from 0.9 to 1.1
     TGraphErrors* band = new TGraphErrors(2);
-    band->SetPoint(0, fForwardSym->GetXaxis()->GetXmin(), 1);
-    band->SetPoint(1, fForward->GetXaxis()->GetXmax(), 1);
+    band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
+    band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
     band->SetPointError(0, 0, .1);
     band->SetPointError(1, 0, .1);
     band->SetFillColor(kYellow+2);
@@ -765,8 +751,8 @@ struct dNdetaDrawer
 #endif
     // Make a nice band from 0.9 to 1.1
     TGraphErrors* band = new TGraphErrors(2);
-    band->SetPoint(0, fForwardSym->GetXaxis()->GetXmin(), 1);
-    band->SetPoint(1, fForward->GetXaxis()->GetXmax(), 1);
+    band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
+    band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
     band->SetPointError(0, 0, .05);
     band->SetPointError(1, 0, .05);
     band->SetFillColor(kYellow+2);
@@ -801,60 +787,24 @@ struct dNdetaDrawer
    * 
    * @return 
    */
-  TH1* GetResult(TList* list, const char* name) const 
+  TH1* FetchResult(const TList* list, const char* name) const 
   {
     if (!list) return 0;
     
-    TH1* ret = static_cast<TH1*>(list->FindObject(name));
-    if (!ret) 
-      Warning("GetResult", "Histogram %s not found", name);
-    
-    return ret;
-  }
-  //__________________________________________________________________
-  /** 
-   * Get the result from previous analysis code 
-   * 
-   * @param fn  File to open 
-   * @param nsd Whether this is NSD
-   * 
-   * @return null or result of previous analysis code 
-   */
-  TH1* GetHHD() 
-  {
-    if (fHHDFile.IsNull()) return 0;
-    const char* fn = fHHDFile.Data();
-    Bool_t nsd = (fTrigString ? fTrigString->GetUniqueID() & 0x4 : false);
-    TDirectory* savdir = gDirectory;
-    if (gSystem->AccessPathName(fn)) { 
-      Warning("GetHHD", "Output of HHD analysis (%s) not available", fn);
-      return 0;
-    }
-    TFile* file = TFile::Open(fn, "READ");
-    if (!file) { 
-      Warning("GetHHD", "couldn't open HHD file %s", fn);
+    TList* all = static_cast<TList*>(list->FindObject("all"));
+    if (!all) { 
+      Warning("GetResult", "No 'all' list find in %s", list->GetName());
+      // list->ls();
       return 0;
     }
-    TString hist(Form("dNdeta_dNdeta%s", (nsd ? "NSD" : "")));
-    TH1* h = static_cast<TH1*>(file->Get(hist.Data()));
-    if (!h) { 
-      Warning("GetHHD", "Couldn't find HHD histogram %s in %s", 
-             hist.Data(), fn);
-      file->Close();
-      savdir->cd();
-      return 0;
+
+    TH1* ret = static_cast<TH1*>(all->FindObject(name));
+    if (!ret) {
+      // all->ls();
+      Warning("GetResult", "Histogram %s not found", name);
     }
-    TH1* r = static_cast<TH1*>(h->Clone("dndeta_hhd"));
-    r->SetTitle("ALICE Forward (HHD)");
-    r->SetFillStyle(0);
-    r->SetFillColor(0);
-    r->SetMarkerStyle(21);
-    r->SetMarkerColor(kPink+1);
-    r->SetDirectory(0);
 
-    file->Close();
-    savdir->cd();
-    return r;
+    return ret;
   }
   //__________________________________________________________________
   /** 
@@ -1135,6 +1085,36 @@ struct dNdetaDrawer
   }
   //__________________________________________________________________
   /** 
+   * Make the ratio of h1 to h2 
+   * 
+   * @param h1 First histogram (numerator) 
+   * @param h2 Second histogram (denominator)
+   * 
+   * @return h1 / h2
+   */
+  TH1* Ratio(const TObject* o1, const TObject* o2, Double_t& max) const
+  {
+    const TH1* h1 = dynamic_cast<const TH1*>(o1); 
+    if (h1) { 
+      const TH1* h2 = dynamic_cast<const TH1*>(o2); 
+      if (h2) return Ratio(h1,h2,max);
+      
+      const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
+      if (g2) return Ratio(h1,g2,max);      
+    }
+    
+    const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
+    if (g1) { 
+      const TGraphAsymmErrors* g2 = dynamic_cast<const TGraphAsymmErrors*>(o2);
+      if (g2) return Ratio(g1, g2, max);
+    }
+
+    Warning("Ratio", "Don't know how to divide a %s (%s) with a %s (%s)", 
+           o1->ClassName(), o1->GetName(), o2->ClassName(), o2->GetName());
+    return 0;
+  }
+  //__________________________________________________________________
+  /** 
    * Compute the ratio of @a h to @a g.  @a g is evaluated at the bin
    * centers of @a h 
    * 
@@ -1340,19 +1320,22 @@ struct dNdetaDrawer
   UShort_t    fRebin;        // Rebinning factor 
   Bool_t      fCutEdges;     // Whether to cut edges
   TString     fTitle;        // Title on plot
-  TString     fHHDFile;      // File name of old results
   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)
-  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
+  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 
   
@@ -1431,7 +1414,6 @@ DrawdNdeta(const char* filename="forward_dndeta.root",
   d.SetRebin(rebin);
   d.SetCutEdges(cutEdges);
   d.SetTitle(title);
-  d.SetHHDFile("");
   d.SetShowOthers(flags & 0x1);
   d.SetShowAlice(flags & 0x2);
   d.SetShowRatios(flags & 0x4);
index 0c6585e..686323f 100644 (file)
@@ -2,7 +2,7 @@
  * Configuration script for forward multiplicity task.  
  *
  * You can copy this to your working directory or to some other
- * directory in up-front the ROOT macro path, and edit it to suit your
+ * directory up-front in your ROOT macro path, and edit it to suit your
  * needs.
  * 
  */
@@ -12,13 +12,19 @@ ForwardAODConfig(AliForwardMultiplicityBase* task)
   if (!task) return;
 
   Info("ForwardAODConfig", "Setting up task %s (%p)", task->GetName(), task);
+
+  // --- General parameters ------------------------------------------
   // Whether to enable low flux specific code 
   task->SetEnableLowFlux(kFALSE);
+
+  // --- Event inspector ---------------------------------------------
   // Set the number of SPD tracklets for which we consider the event a
   // low flux event
   task->GetEventInspector().SetLowFluxCut(1000); 
   // Set the maximum error on v_z [cm]
   task->GetEventInspector().SetMaxVzErr(0.2);
+  
+  // --- Energy Loss Fitter ------------------------------------------
   // Set the eta axis to use - note, this overrides whatever is used
   // by the rest of the algorithms - but only for the energy fitter
   // algorithm. 
@@ -43,16 +49,22 @@ ForwardAODConfig(AliForwardMultiplicityBase* task)
   // Set the minimum number of entries in the distribution before
   // trying to fit to the data
   task->GetEnergyFitter().SetMinEntries(1000);
+
+  // --- Sharing filter ----------------------------------------------
   // Set the low cut used for sharing - overrides settings in eloss fits
   task->GetSharingFilter().SetLowCut(0.3);
   // Set the number of xi's (width of landau peak) to stop at 
   task->GetSharingFilter().SetNXi(1);
+
+  // --- Density calculator 
   // Set the maximum number of particle to try to reconstruct 
   task->GetDensityCalculator().SetMaxParticles(3);
   // Wet whether to use poisson statistics to estimate N_ch
   task->GetDensityCalculator().SetUsePoisson(false);
   // Set the lower multiplicity cut.  Overrides setting in energy loss fits.
   task->GetDensityCalculator().SetMultCut(0.3); //was 0.3
+
+  // --- Corrector ---------------------------------------------------
   // Whether to use the secondary map correction
   task->GetCorrections().SetUseSecondaryMap(true);
   // Whether to use the vertex bias correction
@@ -61,15 +73,32 @@ ForwardAODConfig(AliForwardMultiplicityBase* task)
   task->GetCorrections().SetUseAcceptance(true);
   // Whether to use the merging efficiency correction 
   task->GetCorrections().SetUseMergingEfficiency(false);
+
+  // --- Histogram Collector -----------------------------------------
   // Set the number of extra bins (beyond the secondary map border) 
   task->GetHistCollector().SetNCutBins(2);
   // Set the correction cut, that is, when bins in the secondary map 
   // is smaller than this, they are considered empty 
   task->GetHistCollector().SetCorrectionCut(0.5);
+  // How to calculate the value of overlapping bins. 
+  // Possible values are 
+  //    kStraightMean 
+  //    kStraightMeanNoZero 
+  //    kWeightedMean 
+  //    kLeastError 
+  task->GetHistCollector().SetMergeMethod(AliFMDHistCollector::kStraightMean);
+  // How to find the fiducial area of the secondary maps 
+  // Possible values are 
+  //   kByCut    Only bins larger that cut are trusted 
+  //   kDistance Only bins that are more than half the size of it neighbors
+  task->GetHistCollector().SetFiducialMethod(AliFMDHistCollector::kByCut);
+
+  // --- Debug -------------------------------------------------------
   // Set the overall debug level (1: some output, 3: a lot of output)
   task->SetDebug(0);
   // Set the debug level of a single algorithm 
   // task->GetEventInspector().SetDebug(4);
+
   // --- Set limits on fits the energy -------------------------------
   // Maximum relative error on parameters 
   AliFMDCorrELossFit::ELossFit::fgMaxRelError = .12;
index 3c82c57..f251c08 100644 (file)
@@ -43,6 +43,11 @@ void MakeAOD(const char* esddir,
   TChain* chain = MakeChain("ESD", esddir,true);
   // If 0 or less events is select, choose all 
   if (nEvents <= 0) nEvents = chain->GetEntries();
+  
+  // --- Set the macro path ------------------------------------------
+  gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:"
+                          "$ALICE_ROOT/ANALYSIS/macros",
+                          gROOT->GetMacroPath()));
 
   // --- Creating the manager and handlers ---------------------------
   AliAnalysisManager *mgr  = new AliAnalysisManager("Forward Train", 
@@ -94,7 +99,7 @@ void MakeAOD(const char* esddir,
 
   // --- Add tasks ---------------------------------------------------
   // Physics selection 
-  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  gROOT->LoadMacro("AddTaskPhysicsSelection.C");
   // test HHD
   AddTaskPhysicsSelection(mc, kTRUE, kFALSE);
 
@@ -107,15 +112,15 @@ void MakeAOD(const char* esddir,
   }
 #endif
   if(centrality) {
-    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
+    gROOT->LoadMacro("AddTaskCentrality.C");
     AddTaskCentrality();
   }
   // FMD 
-  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskFMD.C");
+  gROOT->LoadMacro("AddTaskFMD.C");
   AddTaskFMD(mc);
 
   // Central 
-  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskCentral.C");
+  gROOT->LoadMacro("AddTaskCentral.C");
   AddTaskCentral();
   
   // --- Run the analysis --------------------------------------------
index 3d1c45c..88004b2 100644 (file)
@@ -31,6 +31,11 @@ void MakeELossFits(const char* esddir,
   if (nEvents <= 0) nEvents = chain->GetEntries();
   Info("MakeELossFits", "Will analyse %d events", nEvents);
 
+  // --- Set the macro path ------------------------------------------
+  gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:"
+                          "$ALICE_ROOT/ANALYSIS/macros",
+                          gROOT->GetMacroPath()));
+
   // --- Creating the manager and handlers ---------------------------
   AliAnalysisManager *mgr  = new AliAnalysisManager("Forward ELoss Train", 
                                                    "Forward energy loss");
@@ -57,18 +62,15 @@ void MakeELossFits(const char* esddir,
   aodHandler->SetOutputFileName("AliAODs.root");
 
   // --- Add tasks ---------------------------------------------------
-  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  gROOT->LoadMacro("AddTaskPhysicsSelection.C");
   AddTaskPhysicsSelection(mc, kTRUE, kFALSE);
   
-  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
+  gROOT->LoadMacro("AddTaskCentrality.C");
   AddTaskCentrality();
   
   AliFMDEnergyFitterTask* task = new AliFMDEnergyFitterTask("fmdEnergyFitter");
   mgr->AddTask(task);
   
-  
-  
   // --- Make the output container and connect it --------------------
   TString outputfile = AliAnalysisManager::GetCommonFileName();
   AliAnalysisDataContainer* histOut = 
index e4639db..9b7536d 100644 (file)
@@ -44,6 +44,11 @@ void MakedNdeta(const char* aoddir=".",
   // If 0 or less events is select, choose all 
   if (nEvents <= 0) nEvents = chain->GetEntries();
 
+  // --- Set the macro path ------------------------------------------
+  gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:"
+                          "$ALICE_ROOT/ANALYSIS/macros",
+                          gROOT->GetMacroPath()));
+
   // --- Creating the manager and handlers ---------------------------
   AliAnalysisManager *mgr  = new AliAnalysisManager("Forward Train", 
                                                    "Forward dN/deta");
@@ -55,10 +60,10 @@ void MakedNdeta(const char* aoddir=".",
        
   // --- Add tasks ---------------------------------------------------
   // Forward 
-  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskForwarddNdeta.C");
+  gROOT->LoadMacro("AddTaskForwarddNdeta.C");
   AddTaskForwarddNdeta(trig, vzMin, vzMax, centLow, centHigh, true);
   // Central
-  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskCentraldNdeta.C");
+  gROOT->LoadMacro("AddTaskCentraldNdeta.C");
   AddTaskCentraldNdeta(trig, vzMin, vzMax, centLow, centHigh);
 
   
index 684be24..97ca8e8 100755 (executable)
@@ -17,6 +17,10 @@ hhd=1
 comp=1
 cent=0
 tit=
+others=0
+published=1
+ratios=1
+asymm=1 
 dopass1=0
 dopass2=0
 dopass3=0
@@ -50,9 +54,11 @@ Options:
        -g,--gdb                Run in GDB mode             ($gdb)
        -E,--eloss              Run energy loss script      
         -r,--rebin              Rebin factor                ($rebin)
-        -A,--cent               Run centrality task         ($cent)
-       -c,--cent-low           Lower centrality cut        ($centlow)
-       -C,--cent-high          Upper centrality cut        ($centhigh)
+        -C,--use-centrality     Run centrality task         ($cent)
+       -O,--show-older         Show older data             ($others)
+       -J,--show-published     Show ALICE published data   ($published)
+       -R,--show-ratios        Show ratios to other data   ($ratios)
+       -Z,--show-asymmetry     Show asymmetry              ($asymm)
 
 TYPE is a comma or space separated list of 
  
@@ -100,12 +106,10 @@ while test $# -gt 0 ; do
        -1|--pass1|-A|--aod)  dopass1=`toggle $dopass1`   ;; 
        -b|--batch)           batch=`toggle $batch`       ;; 
        -P|--proof)           proof=$2            ; shift ;; 
-       -B|--use-cent)        cent=`toggle $cent` ;;
+       -C|--use-centrality)  cent=`toggle $cent` ;;
        -M|--mc)              mc=`toggle $mc`     ;; 
        -g|--gdb)             gdb=`toggle $gdb`   ;;
        -r|--rebin)           rebin=$2            ; shift ;;
-       -c|--cent-low)        centlow=$2          ; shift ;;
-       -C|--cent-high)       centhigh=$2         ; shift ;;
        -v|--vz-min)          vzmin=$2            ; shift ;; 
        -V|--vz-max)          vzmax=$2            ; shift ;; 
        -E|--eloss)           pass1=MakeELossFits.C 
@@ -113,7 +117,11 @@ while test $# -gt 0 ; do
                              pass3=scripts/DrawAnaELoss.C 
                              output1=forward_eloss.root 
                              dopass2=1 
-                            ;;
+                             ;;
+       -O|--show-older)      others=`toggle $others`   ;;
+       -J|--show-published)  published=`toggle $published`     ;;
+       -R|--show-ratios)     ratios=`toggle $ratios`   ;;
+       -Z|--show-asymmetry)  asymm=`toggle $asymm`     ;;
        -t|--type)           
            #if test "x$type" = "x" ; then type=$2 ; else type="$type|$2"; fi
            type=$2
@@ -188,7 +196,13 @@ fi
 
 if test $dopass3 -gt 0 ; then
     tit=`echo $tit | tr ' ' '@'` 
-    args=(\(\"${output2}\"\,0xf,\"\",$rebin \))
+    flags=0
+    if test $others    -gt 0 ; then let flags=$(($flags|0x1)); fi
+    if test $published -gt 0 ; then let flags=$(($flags|0x2)); fi
+    if test $ratios    -gt 0 ; then let flags=$(($flags|0x4)); fi
+    if test $asymm     -gt 0 ; then let flags=$(($flags|0x8)); fi
+
+    args=(\(\"${output2}\"\,${flags},\"\",$rebin \))
     if test "x$pass1" = "xMakeELossFits.C" ; then 
        args=(\(\"${output1}\"\))
     fi