A better way to specify the Nch axis for the MultDists task.
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Sep 2013 12:55:55 +0000 (12:55 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Sep 2013 12:55:55 +0000 (12:55 +0000)
We can now easily make the same binning as used by CMS
for example.  This is done by introducing the nested class
AliForwardMultDists::BinSpec which provides a simple
interface to define Nch axis ranges.

MC Nch axis is computed so we have integer bins.

PWGLF/FORWARD/analysis2/AddTaskForwardMultDists.C
PWGLF/FORWARD/analysis2/AliForwardMultDists.cxx
PWGLF/FORWARD/analysis2/AliForwardMultDists.h
PWGLF/FORWARD/analysis2/scripts/UnfoldMultDists.C
PWGLF/PWGLFforward2LinkDef.h

index 44de299e91944e2811006dbadcfaa288274c3d2b..d1398298d05579ec49a5f29072a574698ee26356 100644 (file)
  *
  * @ingroup pwglf_forward_dndeta
  */
+
 AliAnalysisTask*
 AddTaskForwardMultDists(const char* trig      = "V0AND", 
                        Double_t    vzMin     = -4, 
                        Double_t    vzMax     = +4,
-                       UShort_t    maxN      = 150,
-                       UShort_t    nDiv      = 1,
                        Bool_t      usePhiAcc = true,
                        Bool_t      useAsymm  = false)
 {
@@ -54,26 +53,76 @@ AddTaskForwardMultDists(const char* trig      = "V0AND",
   task->SetIpZRange(vzMin, vzMax);
   // Set the trigger mask to use (INEL,INEL>0,NSD)
   task->SetTriggerMask(trig);
-  // Set the maximum number of charged particles
-  task->SetMaxN(maxN);
-  // Set number of divisions per particle number 
-  task->SetNDivisions(nDiv);
   // Set whether to use stored phi acceptance 
   task->SetUsePhiAcc(usePhiAcc);
   // add the task 
   mgr->AddTask(task);
 
+  // Variable size axis objects
+  // for |eta|<0.5 from CMD
+  AliForwardMultDists::BinSpec b05(-0.5, 0.5, -0.5);
+  b05.Push(21, 1);
+  b05.Push(1, 3); 
+  b05.Push(1, 5);
+  b05.Push(3, 5); // <-- Extra 
+
+  // for |eta|<0.5 from ALICE
+  AliForwardMultDists::BinSpec a05(-0.5, 0.5, -0.5);
+  a05.Push(21, 1);
+  a05.Push(3, 2);
+
+  // For |eta|<1 from CMS 
+  AliForwardMultDists::BinSpec b10(-1.0, 1.0, -0.5);
+  b10.Push(35, 1);
+  b10.Push(1, 3);
+  b10.Push(4, 5);
+  b10.Push(2, 5); // <-- Extra 
+
+  // For |eta|<1 from ALICE 
+  AliForwardMultDists::BinSpec a10(-1.0, 1.0, -0.5);
+  a10.Push(41, 1);
+  a10.Push(1,  2);
+
+  // For |eta|<1.3 from ALICE 
+  AliForwardMultDists::BinSpec a13(-1.3, 1.3, -0.5);
+  a13.Push(41, 1);
+  a13.Push(7,  2);
+
+  // For |eta|<1.5 from CMS 
+  AliForwardMultDists::BinSpec b15(-1.5, 1.5, -0.5);
+  b15.Push(46, 1);
+  b15.Push(1,  2);
+  b15.Push(1,  3); 
+  b15.Push(4,  5);
+  b15.Push(4,  5); // <-- Extra
+
+  // for |eta|<2.0 from CMS 
+  AliForwardMultDists::BinSpec b20(-2.0, 2.0, -0.5);
+  b20.Push(62, 1);
+  b20.Push(2,  2);
+  b20.Push(1,  4);
+  b20.Push(3,  10);
+  b20.Push(2,  10); // <-- Extra
+
+  // for |eta|<2 from CMS 
+  AliForwardMultDists::BinSpec b24(-2.4, 2.4, -0.5);
+  b24.Push(57, 1);
+  b24.Push(1,  2);
+  b24.Push(1,  3); 
+  b24.Push(3,  10);
+  b24.Push(4,  10); // <-- Extra
+  
   // Add bins 
-  Double_t limits[] = { 1, 1.5, 2, 3, 0 };
-  Double_t* plimit  = limits;
-  while (*plimit) {
-    Double_t eta = *plimit;
-    task->AddBin(-eta, +eta);
+  AliForwardMultDists::BinSpec*  bs[] = { &b05, &b10, &b15, &b20, &b24, 0 };
+  AliForwardMultDists::BinSpec** pb   = bs;
+  while (*pb) {
+    AliForwardMultDists::BinSpec* b = *pb;
+    task->AddBin(*b);
     if (useAsymm) { 
-      task->AddBin(-eta, 0);
-      task->AddBin(0,    +eta);
+      task->AddBin(b->fEtaMin, 0,          b->Axis());
+      task->AddBin(0,          b->fEtaMax, b->Axis());
     }
-    plimit++;
+    pb++;
   }
 
   // --- create containers for input/output --------------------------
index 30865266cbac217001b061458ed4a6a3a7f1d947..29dd01dcff6f26ede91c69aa4a2641db609f80c5 100644 (file)
@@ -9,9 +9,60 @@
 #include <TVector2.h>
 #include <THStack.h>
 #include <TParameter.h>
+#include <TError.h>
 #include <iostream>
 #include <iomanip> 
 
+//____________________________________________________________________
+AliForwardMultDists::BinSpec::BinSpec(Double_t etaMin, 
+                                     Double_t etaMax, 
+                                     Double_t nchLow)
+  : fEtaMin(etaMin), fEtaMax(etaMax), fLow(nchLow), fN(0), fD(0), fAxis(1,0,1)
+{
+}
+//____________________________________________________________________
+void
+AliForwardMultDists::BinSpec::Push(UShort_t n, Double_t d)
+{
+  Int_t l = fN.GetSize();
+  fN.Set(l+1); // Not terribly efficient, but that's ok because we do
+  fD.Set(l+1); // this infrequently 
+  fN[l] = n;
+  fD[l] = d;
+}
+//____________________________________________________________________
+const TAxis&
+AliForwardMultDists::BinSpec::Axis() const
+{
+  if (fAxis.GetNbins() > 1) return fAxis;
+  if (fN.GetSize() <= 0) return fAxis;
+  if (fN.GetSize() == 1) { 
+    // Exactly one spec, 
+    fAxis.Set(fN[0], fLow, fN[0]*fD[0]+fLow);
+  }
+  else { 
+    Int_t n = fN.GetSum();
+    Int_t l = fN.GetSize();
+    Int_t i = 0;
+    TArrayD b(n+1);
+    b[0]    = fLow;
+    for (Int_t j = 0; j < l; j++) { 
+      for (Int_t k = 0; k < fN[j]; k++) {
+       b[i+1] = b[i] + fD[j];
+       i++;
+      }
+    }
+    if (i != n) { 
+      ::Warning("Axis", "Assigned bins, and number of bins not equal");
+      n = TMath::Min(i, n);
+    }
+    fAxis.Set(n, b.GetArray());
+  }
+  return fAxis;
+}
+      
+
+
 //____________________________________________________________________
 AliForwardMultDists::AliForwardMultDists()
   : AliAnalysisTaskSE(), 
@@ -32,8 +83,6 @@ AliForwardMultDists::AliForwardMultDists()
     fForwardCache(0),
     fCentralCache(0),
     fMCCache(0),
-    fMaxN(200),
-    fNDivisions(1),
     fUsePhiAcc(true)
 {}
 
@@ -57,8 +106,6 @@ AliForwardMultDists::AliForwardMultDists(const char* name)
     fForwardCache(0),
     fCentralCache(0),
     fMCCache(0),
-    fMaxN(200),
-    fNDivisions(1),
     fUsePhiAcc(true)
 {
   /** 
@@ -90,8 +137,6 @@ AliForwardMultDists::AliForwardMultDists(const AliForwardMultDists& o)
     fForwardCache(o.fForwardCache),
     fCentralCache(o.fCentralCache),
     fMCCache(o.fMCCache),
-    fMaxN(o.fMaxN),
-    fNDivisions(o.fNDivisions),
     fUsePhiAcc(o.fUsePhiAcc)
 {}
 //____________________________________________________________________
@@ -116,23 +161,11 @@ AliForwardMultDists::operator=(const AliForwardMultDists& o)
   fForwardCache                = o.fForwardCache;
   fCentralCache                = o.fCentralCache;
   fMCCache             = o.fMCCache;
-  fMaxN                        = o.fMaxN;
-  fNDivisions           = o.fNDivisions;
   fUsePhiAcc           = o.fUsePhiAcc;
   
   return *this;
 }
 
-//____________________________________________________________________
-void
-AliForwardMultDists::SetMaxN(UShort_t maxN, 
-                            UShort_t nDivisions)
-{
-  const UShort_t one = 1;
-  fMaxN              = maxN;
-  fNDivisions        = TMath::Max(nDivisions, one); 
-}
-
 //____________________________________________________________________
 void
 AliForwardMultDists::UserCreateOutputObjects()
@@ -402,7 +435,7 @@ AliForwardMultDists::Terminate(Option_t* /*option=""*/)
   EtaBin*  bin    = 0;
   TIter    next(&fBins);
   while ((bin = static_cast<EtaBin*>(next()))) {
-    bin->Terminate(in, out, fMaxN);
+    bin->Terminate(in, out);
 
     sta = oth;
     if      (bin->IsSymmetric()) sta = sym;
@@ -451,7 +484,6 @@ AliForwardMultDists::StoreInformation(const AliAODForwardMult* aod)
   fList->Add(AliForwardUtil::MakeParameter("trigger", ULong_t(fTriggerMask)));
   fList->Add(AliForwardUtil::MakeParameter("minIpZ", fMinIpZ));
   fList->Add(AliForwardUtil::MakeParameter("maxIpZ", fMaxIpZ));
-  fList->Add(AliForwardUtil::MakeParameter("maxN", UShort_t(fMaxN)));
   fList->Add(AliForwardUtil::MakeParameter("count", UShort_t(1)));
 }
 
@@ -539,7 +571,7 @@ AliForwardMultDists::SetupForData(const TH2& hist, Bool_t useMC)
   TIter   next(&fBins);
   EtaBin* bin = 0;
   while ((bin = static_cast<EtaBin*>(next()))) {
-    bin->SetupForData(fList, hist, fMaxN, fNDivisions, useMC);
+    bin->SetupForData(fList, hist, useMC);
   }
     
 }
@@ -615,7 +647,15 @@ AliForwardMultDists::ProjectX(const TH2* input, TH1* cache)
 }
 //____________________________________________________________________
 void
-AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax) 
+AliForwardMultDists::AddBin(const BinSpec& spec)
+{
+  AddBin(spec.fEtaMin, spec.fEtaMax, spec.Axis());
+}
+
+//____________________________________________________________________
+void
+AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax, 
+                           const TAxis& mAxis) 
 {
   /** 
    * Add an @f$\eta@f$ bin
@@ -628,7 +668,28 @@ AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax)
     etaLow = -TMath::Abs(etaLow);
     etaMax = +TMath::Abs(etaLow);
   }
-  EtaBin* bin = new EtaBin(etaLow, etaMax);
+  EtaBin* bin = new EtaBin(etaLow, etaMax, mAxis);
+  // AliInfoF("Adding bin %f, %f: %s", etaLow, etaMax, bin->GetName());
+  fBins.Add(bin);
+}
+//____________________________________________________________________
+void
+AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax, 
+                           UShort_t nMax, UShort_t nDiv)
+{
+  /** 
+   * Add an @f$\eta@f$ bin
+   * 
+   * @param etaLow Low cut on @f$\eta@f$
+   * @param etaMax High cut on @f$\eta@f$
+   */
+  // Symmetric bin
+  if (etaMax >= kInvalidEta) { 
+    etaLow = -TMath::Abs(etaLow);
+    etaMax = +TMath::Abs(etaLow);
+  }
+  TAxis mAxis((nMax+1)/nDiv, -.5, nMax+.5);
+  EtaBin* bin = new EtaBin(etaLow, etaMax, mAxis);
   // AliInfoF("Adding bin %f, %f: %s", etaLow, etaMax, bin->GetName());
   fBins.Add(bin);
 }
@@ -656,7 +717,6 @@ AliForwardMultDists::Print(Option_t* /*option=""*/) const
            << "  Trigger mask:        " 
            << AliAODForwardMult::GetTriggerString(fTriggerMask) << "\n"
            << "  IpZ range:           " << fMinIpZ <<" - "<< fMaxIpZ << "\n"
-           << "  Max Nch:             " << fMaxN << "\n"
            << "  Use phi acceptance:  " << fUsePhiAcc << "\n"
            << "  Bins:" << std::endl;
   TIter   next(&fBins);
@@ -669,6 +729,8 @@ AliForwardMultDists::Print(Option_t* /*option=""*/) const
 //====================================================================
 AliForwardMultDists::EtaBin::EtaBin() 
   : fName(""),
+    fMAxis(1,0,0),
+    fTAxis(1,0,0),
     fMinEta(0),
     fMaxEta(0),
     fMinBin(0), 
@@ -686,8 +748,11 @@ AliForwardMultDists::EtaBin::EtaBin()
 }
 
 //____________________________________________________________________
-AliForwardMultDists::EtaBin::EtaBin(Double_t minEta, Double_t maxEta) 
+AliForwardMultDists::EtaBin::EtaBin(Double_t minEta, Double_t maxEta,
+                                   const TAxis& mAxis) 
   : fName(""),
+    fMAxis(1,0,0),
+    fTAxis(1,0,0),
     fMinEta(minEta), 
     fMaxEta(maxEta),
     fMinBin(0), 
@@ -709,11 +774,40 @@ AliForwardMultDists::EtaBin::EtaBin(Double_t minEta, Double_t maxEta)
   fName.ReplaceAll("-", "m");
   fName.ReplaceAll("+", "p");
   fName.ReplaceAll(".", "d");
+
+  // Copy to other our object
+  mAxis.Copy(fMAxis);
+  
+  if (mAxis.GetXbins() && mAxis.GetXbins()->GetArray()) {
+    const TArrayD& mA  = *(mAxis.GetXbins());
+    TArrayD        tA(mA.GetSize());
+    Int_t          j   = 0;
+    Double_t       min = mA[0];
+    tA[0]              = min;
+    for (Int_t i = 1; i < mA.GetSize(); i++) { 
+      Double_t d = mA[i] - min;
+      if (d < 1) 
+       // Not full integer bin
+       continue;  
+      tA[j+1] = tA[j] + Int_t(d);
+      min     = tA[j+1];
+      j++;
+    }
+    fTAxis.Set(j, tA.GetArray());
+  }
+  else {
+    // Rounded down maximum and minimum
+    Int_t max = mAxis.GetXmax(); 
+    Int_t min = mAxis.GetXmin();
+    fTAxis.Set((max-min)+1, min-.5, max+.5);
+  }
 }
 //____________________________________________________________________
 AliForwardMultDists::EtaBin::EtaBin(const EtaBin& o) 
   : TObject(o),
     fName(o.fName),
+    fMAxis(o.fMAxis), 
+    fTAxis(o.fTAxis),
     fMinEta(o.fMinEta),
     fMaxEta(o.fMaxEta),
     fMinBin(o.fMinBin), 
@@ -732,6 +826,8 @@ AliForwardMultDists::EtaBin::operator=(const EtaBin& o)
   if (&o == this) return *this;
   
   fName                 = o.fName;
+  fMAxis         = o.fMAxis;
+  fTAxis         = o.fTAxis;
   fMinEta       = o.fMinEta;
   fMaxEta       = o.fMaxEta;
   fMinBin       = o.fMinBin; 
@@ -802,10 +898,59 @@ AliForwardMultDists::EtaBin::FindParent(TList* l, Bool_t create) const
   return p;
 }
 
+//____________________________________________________________________
+TH1*
+AliForwardMultDists::EtaBin::CreateH1(const char*  name, 
+                                     const char*  title, 
+                                     const TAxis& xAxis)
+{
+  TH1* ret = 0;
+  
+  if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray()) 
+    ret = new TH1D(name,title,xAxis.GetNbins(), xAxis.GetXbins()->GetArray());
+  else 
+    ret = new TH1D(name,title,xAxis.GetNbins(),xAxis.GetXmin(),xAxis.GetXmax());
+  return ret;
+}
+//____________________________________________________________________
+TH2*
+AliForwardMultDists::EtaBin::CreateH2(const char*  name, 
+                                     const char*  title, 
+                                     const TAxis& xAxis,
+                                     const TAxis& yAxis)
+{
+  TH2* ret = 0;
+  
+  if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray()) {
+    if (yAxis.GetXbins() && yAxis.GetXbins()->GetArray()) 
+      // Variable variable
+      ret = new TH2D(name,title, 
+                    xAxis.GetNbins(), xAxis.GetXbins()->GetArray(),
+                    yAxis.GetNbins(), yAxis.GetXbins()->GetArray());
+    else 
+      // variable fixed
+      ret = new TH2D(name,title, 
+                    xAxis.GetNbins(), xAxis.GetXbins()->GetArray(),
+                    yAxis.GetNbins(),yAxis.GetXmin(),yAxis.GetXmax());
+  }
+  else {
+    if (yAxis.GetXbins() && yAxis.GetXbins()->GetArray()) 
+      // fixed variable 
+      ret = new TH2D(name,title, 
+                    xAxis.GetNbins(), xAxis.GetXmin(), xAxis.GetXmax(),
+                    yAxis.GetNbins(), yAxis.GetXbins()->GetArray());
+    else 
+      // fixed fixed
+      ret = new TH2D(name,title, 
+                    xAxis.GetNbins(), xAxis.GetXmin(), xAxis.GetXmax(),
+                    yAxis.GetNbins(), yAxis.GetXmin(), yAxis.GetXmax());
+  }
+  return ret;
+}
+
 //____________________________________________________________________
 void
 AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist, 
-                                         UShort_t max, UShort_t nDiv, 
                                          Bool_t useMC)
 {
   /** 
@@ -842,9 +987,8 @@ AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
   fMinBin = hist.GetXaxis()->FindBin(fMinEta);
   fMaxBin = hist.GetXaxis()->FindBin(fMaxEta-.00001);
 
-  TString eta(Form("%+5.2f<#eta<%+5.2f", fMinEta, fMaxEta));
-  fSum = new TH1D("rawDist", Form("Raw P(#it{N}_{ch}) in %s", eta.Data()), 
-                 (max+1)*nDiv, -.5, max+.5);
+  TString t(Form("%+5.2f<#eta<%+5.2f", fMinEta, fMaxEta));
+  fSum = CreateH1("rawDist",Form("Raw P(#it{N}_{ch}) in %s",t.Data()),fMAxis);
   fSum->SetDirectory(0);
   fSum->GetXaxis()->SetTitle("#it{N}_{ch}");
   fSum->GetYaxis()->SetTitle("Raw P(#it{N}_{ch})");
@@ -857,8 +1001,7 @@ AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
   fSum->Sumw2();
   l->Add(fSum);
   
-  fCorr = new TH2D("corr", Form("C_{SPD,FMD} in %s", eta.Data()),
-                  max+2, -1.5, max+.5, max+2, -1.5, max+.5);
+  fCorr = CreateH2("corr",Form("C_{SPD,FMD} in %s", t.Data()),fTAxis,fTAxis);
   fCorr->SetDirectory(0);
   fCorr->GetXaxis()->SetTitle("Forward");
   fCorr->GetYaxis()->SetTitle("Central");
@@ -874,16 +1017,15 @@ AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
   l->Add(fCoverage);
 
   if (!useMC) return;
-  fResponse = new TH2D("response", Form("Reponse matrix in %s", eta.Data()),
-                      nDiv*(max+1), -0.5, max+0.5, max+1, -.5, max+.5);
+  fResponse = CreateH2("response", Form("Reponse matrix in %s", t.Data()),
+                      fMAxis, fTAxis);
   fResponse->SetDirectory(0);
   fResponse->SetYTitle("MC");
   fResponse->SetXTitle("Signal");
   fResponse->SetOption("colz");
   l->Add(fResponse);
   
-  fTruth = new TH1D("truth", Form("True P(#it{N}_{ch}) in %s", eta.Data()),
-                   (max+1), -0.5, max+0.5);
+  fTruth = CreateH1("truth",Form("True P(#it{N}_{ch}) in %s",t.Data()),fTAxis);
   fTruth->SetXTitle(fSum->GetXaxis()->GetTitle());
   fTruth->SetYTitle(fSum->GetYaxis()->GetTitle());
   fTruth->SetLineColor(kBlack);
@@ -898,7 +1040,7 @@ AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
 
   fTruthAccepted = static_cast<TH1D*>(fTruth->Clone("truthAccepted"));
   fTruthAccepted->SetTitle(Form("True (accepted) P(#it{N}_{ch}) in %s", 
-                               eta.Data()));
+                               t.Data()));
   fTruthAccepted->SetLineColor(kGray+2);
   fTruthAccepted->SetFillColor(kOrange+2);
   fTruthAccepted->SetDirectory(0);
@@ -996,7 +1138,7 @@ AliForwardMultDists::EtaBin::Process(const TH1& sumForward,
     // Fill with weight 
     Double_t w = 1; // e2sum > 0 ? 1/TMath::Sqrt(e2sum) : 1
     fSum->Fill(sum, w);
-    fCorr->Fill(fsum, csum);
+    fCorr->Fill((fsum<0?0:fsum), (csum<0?0:csum));
     
     Int_t nTotal = fMaxBin - fMinBin + 1;
     fCoverage->Fill(100*float(covered) / nTotal);
@@ -1018,7 +1160,7 @@ AliForwardMultDists::EtaBin::Process(const TH1& sumForward,
 }
 //____________________________________________________________________
 void
-AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out, UShort_t maxN)
+AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out)
 {
   /** 
    * Called at the end of the final processing of the job on the
@@ -1054,6 +1196,7 @@ AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out, UShort_t maxN)
   while (*phist) { 
     TH1* h = *phist;
     if (h) { 
+      Int_t    maxN = h->GetNbinsX();
       Double_t intg = h->Integral(1, maxN);
       h->Scale(1. / intg, "width");
     }
index f5e6517b2c04bbfba8e11e491a8c2cb015f83513..8e44fc84a8f26fb7bb90543aa161aa5366f3fcc6 100644 (file)
@@ -24,6 +24,52 @@ public:
     kVertex   = 4, 
     kTriggerVertex = 5
   };
+  /**
+   * Structure to define @f$\eta@f$ bins with an @f$ N_{ch}@f$ axis 
+   * 
+   */
+  struct BinSpec 
+  {
+    /** 
+     * Constructor
+     * 
+     * @param etaMin Low cut on @f$\eta@f$
+     * @param etaMax High cut on @f$\eta@f$
+     * @param nchLow Lowest @f$ N_{ch}@f$ (e.g., -0.5);
+     */
+    BinSpec(Double_t etaMin, Double_t etaMax, Double_t nchLow);
+    /** 
+     * Push @a n bins of with @a d in @f$ N_{ch}@f$ onto the axis 
+     * 
+     * If only a single push is done, then we will get an axis of
+     * equally sized bins (@a d) from @f$ l@f$ edge to @f$ nd+low@f$ -
+     * e.g., if one only does 
+     *
+     * @code
+     * BinSpec b(e1, e2, -.5);
+     * b.Push(10, 1);
+     * @endcode 
+     *
+     * One ends up with 10 bins from -0.5 to 9.5. 
+     *
+     * @param n Number of bins to push
+     * @param d Bin width of each of the bins
+     */
+    void Push(UShort_t n, Double_t d);
+    /** 
+     * Get the axis computed from the setup using Push
+     * 
+     * @return Reference to the axis 
+     */
+    const TAxis& Axis() const;
+    Double_t fEtaMin;
+    Double_t fEtaMax;
+    Double_t fLow;
+    TArrayI  fN;
+    TArrayD  fD;
+    mutable TAxis    fAxis;
+  };
+
   /** 
    * Default constructor
    */
@@ -98,22 +144,30 @@ public:
    * @param cache 
    */
   static void ProjectX(const TH2* input, TH1* cache);
+  /** 
+   * Add an @f$\eta@f$ bin
+   * 
+   * @param spec Bin specification 
+   */
+  void AddBin(const BinSpec& spec);
   /** 
    * Add an @f$\eta@f$ bin
    * 
    * @param etaLow Low cut on @f$\eta@f$
    * @param etaMax High cut on @f$\eta@f$
+   * @param nAxis  Axis to use for measured @f$ N_{ch}@f$ 
    */
-  void AddBin(Double_t etaLow, Double_t etaMax=kInvalidEta); 
+  void AddBin(Double_t etaLow, Double_t etaMax, const TAxis& nAxis); 
   /** 
-   * Set the maximum @f$N_{ch}@f$ to consider 
+   * Add an @f$\eta@f$ bin
    * 
-   * @param maxN Maximum
-   * @param nBins (Optionally) the number of bins per particle number 
+   * @param etaLow Low cut on @f$\eta@f$
+   * @param etaMax High cut on @f$\eta@f$
+   * @param nMax   Maximum @f$ N_{ch}@f$ 
+   * @param nDiv   Number of subdivisions per @f$ N_{ch}@f$
    */
-  void SetMaxN(UShort_t maxN, UShort_t nBins=-1);
-  void SetNDivisions(UInt_t nDiv) { fNDivisions = nDiv; }
-  /** 
+  void AddBin(Double_t etaLow, Double_t etaMax, UShort_t nMax, UShort_t nDiv); 
+ /** 
    * Set the range of valid interaction points 
    * 
    * @param z1 Least Z coordinate 
@@ -158,8 +212,10 @@ public:
      * 
      * @param minEta Least @f$\eta@f$ to consider 
      * @param maxEta Largest @f$\eta@f$ to consider 
+     * @param mAxis  The @f$ N_{ch}@f$ axis to use for measured data
+     * @param tAxis  The @f$ N_{ch}@f$ axis to use for truth data
      */
-    EtaBin(Double_t minEta, Double_t maxEta); 
+    EtaBin(Double_t minEta, Double_t maxEta, const TAxis& mAxis); 
     /** 
      * Copy constructor
      *
@@ -211,17 +267,37 @@ public:
      * @return Container, or null
      */
     TList* FindParent(TList* l, Bool_t create=true) const;
+    /** 
+     * Create a 1D histogram with specified axis 
+     * 
+     * @param name  Name of histogram 
+     * @param title Title of histogram 
+     * @param xAxis X-axis to use 
+     * 
+     * @return Created histogram
+     */
+    static TH1* CreateH1(const char* name, const char* title, 
+                        const TAxis& xAxis);
+    /** 
+     * Create a 2D histogram with specified axis 
+     * 
+     * @param name  Name of histogram 
+     * @param title Title of histogram 
+     * @param xAxis X-axis to use 
+     * @param yAxis Y-axis to use 
+     * 
+     * @return Created histogram
+     */
+    static TH2* CreateH2(const char* name, const char* title, 
+                        const TAxis& xAxis, const TAxis& yAxis);
     /** 
      * Set-up internal structures on first event. 
      * 
      * @param list  List to add information to
      * @param hist  Template histogram 
-     * @param max   Maximum number of particles 
-     * @param nDiv  Number of divisions per charged particle bin
      * @param useMC Whether to set-up for MC input 
      */
-    void SetupForData(TList* list, const TH2& hist, UShort_t max, 
-                     UShort_t nDiv, Bool_t useMC);
+    void SetupForData(TList* list, const TH2& hist, Bool_t useMC);
     /** 
      * Process a single event 
      * 
@@ -241,11 +317,12 @@ public:
      * 
      * @param in    Input list
      * @param out   Output list 
-     * @param maxN  Maximum number of @f$N_{ch}@f$ to consider
      */
-    void Terminate(TList* in, TList* out, UShort_t maxN);
+    void Terminate(TList* in, TList* out);
       
     TString  fName;          // Name of this bin
+    TAxis    fMAxis;         // Axis used for measured Nch
+    TAxis    fTAxis;         // Axis used for true Nch
     Double_t fMinEta;        // Least @f$\eta@f$ to consider
     Double_t fMaxEta;        // Largest @f$\eta@f$ to consider
     Int_t    fMinBin;        // Least @f$\eta@f$ bin to consider
@@ -257,7 +334,7 @@ public:
     TH1*     fTruthAccepted; // `true' distribution for accepted events
     TH1*     fCoverage;      // How much was covered
 
-    ClassDef(EtaBin,1);
+    ClassDef(EtaBin,2);
   };
   TList    fBins;         // List of bins 
   TList*   fSymmetric;    // Bins symmetric around 0
@@ -276,8 +353,6 @@ public:
   TH1*     fForwardCache; // Projection cache 
   TH1*     fCentralCache; // Projection cache 
   TH1*     fMCCache;      // Projection cache 
-  UShort_t fMaxN;         // Maximum of @f$N_{ch}@f$ 
-  UShort_t fNDivisions;   // Number of particle number sub-divions
   Bool_t   fUsePhiAcc;    // If true, scale by phi acceptance 
 
   ClassDef(AliForwardMultDists,1);
index 4ceb8853bf218d85f035fea8d08d6c842e3dc0ed..168284661e5711800f5afc6f2fab861c8fce5a75 100644 (file)
@@ -280,13 +280,11 @@ struct Unfolder
     // Get some info from the input collection 
     UShort_t sys;
     UShort_t sNN; 
-    UShort_t maxN; 
     ULong_t  trig; 
     Double_t minZ; 
     Double_t maxZ;
     GetParameter(mTop, "sys",     sys);          
     GetParameter(mTop, "snn",     sNN);          
-    GetParameter(mTop, "maxN",    maxN);         
     GetParameter(mTop, "trigger", trig);
     GetParameter(mTop, "minIpZ",  minZ); 
     GetParameter(mTop, "maxIpZ",  maxZ); 
@@ -307,7 +305,8 @@ struct Unfolder
     if (mId == 0xDeadBeef) return;
 
     // Store information 
-    SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,maxN);
+    SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,
+                   corrFile.IsNull());
 
     // Load other data 
     TString savPath(gROOT->GetMacroPath());
@@ -336,6 +335,12 @@ struct Unfolder
 
     SaveSummarize();
   }
+  /** 
+   * Append an & to a string and the next term.
+   * 
+   * @param trg  Output string
+   * @param what Term
+   */
   static void AppendAnd(TString& trg, const TString& what)
   {
     if (!trg.IsNull()) trg.Append(" & ");
@@ -353,7 +358,6 @@ struct Unfolder
    * @param trigger  Trigger mask 
    * @param minIpZ   Least z coordinate of interaction point
    * @param maxIpZ   Largest z coordinate of interaction point
-   * @param maxN     Largest @f$N_{ch}@f$ to consider 
    */
   void SaveInformation(TDirectory* dir, 
                       const TString& method,
@@ -364,15 +368,20 @@ struct Unfolder
                       UInt_t         trigger,
                       Double_t       minIpZ, 
                       Double_t       maxIpZ, 
-                      UShort_t       maxN) const
+                      Bool_t         self) const
   {
     dir->cd();
 
+    TParameter<bool>* pM = new TParameter<bool>("self", self);
+    pM->SetBit(BIT(19));
+    pM->Write();
+
     TNamed* outMeth = new TNamed("method", method.Data());
     outMeth->SetUniqueID(mId);
     outMeth->Write();
 
     TParameter<double>* pR = new TParameter<double>("regParam", regParam);
+    pR->SetBit(BIT(19));
     pR->Write();
     
     TString tS = (sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "?");
@@ -435,13 +444,12 @@ struct Unfolder
     pT->Write();
     
     TParameter<double>* pY = new TParameter<double>("minIpZ", minIpZ);
+    pY->SetBit(BIT(19));
     pY->Write();
 
     TParameter<double>* pZ = new TParameter<double>("maxIpZ", maxIpZ);
+    pZ->SetBit(BIT(19));
     pZ->Write();
-
-    TParameter<int>* pN = new TParameter<int>("maxN", maxN);
-    pN->Write();
   }
   /** 
    * Save a script to do a summary of this step 
@@ -493,6 +501,8 @@ struct Unfolder
                                         "Unfolded P(#it{N}_{ch})");
     THStack*  allCorrected = new THStack("corrected",
                                         "Corrected P(#it{N}_{ch})");
+    THStack*  allRatio     = (sys != 1 ? 0 : 
+                             new THStack("ratios", "Ratios to other"));
     TMultiGraph* allALICE  = (sys != 1 ? 0 : 
                              new TMultiGraph("alice", "ALICE Published"));
     TMultiGraph* allCMS    = (sys != 1 ? 0 : 
@@ -525,9 +535,14 @@ struct Unfolder
       THStack* binS = ProcessBin(mBin, cBin, method, r, dir);
       if (!binS) continue;
 
+      TH1* result = 0;
       Bin2Stack(binS, i, allMeasured, allTruth, allTruthA, 
-               allUnfolded, allCorrected);
-      Other2Stack(o->GetName(), i, sNN, allALICE, allCMS);
+               allUnfolded, allCorrected, result);
+
+      TGraph* alice = 0;
+      TGraph* cms   = 0;
+      Other2Stack(o->GetName(), i, sNN, allALICE, allCMS, alice, cms);
+      Ratio2Stack(i, result, alice, cms, allRatio);
       i++;
     }
     dir->Add(allMeasured);
@@ -547,6 +562,12 @@ struct Unfolder
       else 
        delete allCMS;
     }
+    if (allRatio && allRatio->GetHists()) { 
+      if (allRatio->GetHists()->GetEntries() > 0) 
+       dir->Add(allRatio);
+      else 
+       delete allRatio;
+    }
   }
   /** 
    * Process a single eta bin 
@@ -749,7 +770,8 @@ struct Unfolder
                 THStack* truth, 
                 THStack* accepted, 
                 THStack* unfolded,
-                THStack* corrected)
+                THStack* corrected,
+                TH1*&    result)
   {
     Int_t open, closed;
     Double_t factor; 
@@ -776,6 +798,7 @@ struct Unfolder
       cln->SetMarkerStyle(sty);
       cln->SetMarkerSize(size);
       cln->Scale(factor); // Scale by 10^i
+      if (col == kColorCorrected) result = cln;
 
       // Make sure we do not get the old legend 
       TObject* tst = cln->FindObject("legend");
@@ -816,7 +839,8 @@ struct Unfolder
    * @param allCMS   Stack of CMS data 
    */
   void Other2Stack(const TString& name, Int_t i,
-                  UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS) 
+                  UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS,
+                  TGraph*& alice, TGraph*& cms) 
   {
     if (!allALICE && !allCMS) return;
 
@@ -824,6 +848,7 @@ struct Unfolder
     tmp.ReplaceAll("p", "+");
     tmp.ReplaceAll("m", "-");
     tmp.ReplaceAll("_", " ");
+    tmp.ReplaceAll("d", ".");
     TObjArray* tokens = tmp.Tokenize(" ");
     if (!tokens || tokens->GetEntriesFast() < 2) { 
       Error("Other2Stack", "Failed to decode eta range from %s", name.Data());
@@ -851,6 +876,7 @@ struct Unfolder
        g->SetMarkerColor(kColorALICE);
        g->SetMarkerSize(size);
        allALICE->Add(g, "p same");
+       alice = g;
       }
     }
     if (allCMS) {
@@ -860,11 +886,74 @@ struct Unfolder
        g->SetMarkerColor(kColorCMS);
        g->SetMarkerSize(size);
        allCMS->Add(g, "p same");
+       cms = g;
+      }
+    }
+  }
+  /** 
+   * Create ratios to other data 
+   * 
+   * @param i 
+   * @param res 
+   * @param alice 
+   * @param cms 
+   * @param all 
+   */
+  void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
+  {
+    if (!all || !res || !(alice || cms)) return;
+
+    Int_t        off  = 5*ib;
+    TGraph*      gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
+    TGraph**     pg   = gs;
+    while (*pg) { 
+      TGraph*     g = *pg;
+      const char* n = (g == alice ? "ALICE" : "CMS");
+
+      TH1*    r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
+      TString tit(r->GetTitle());
+      tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
+      r->SetTitle(tit);
+      r->SetMarkerColor(g->GetMarkerColor());
+      r->SetLineColor(g->GetLineColor());
+
+      TObject* tst = r->FindObject("legend");
+      if (tst) r->GetListOfFunctions()->Remove(tst);
+
+      for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
+       Double_t c = r->GetBinContent(i);
+       Double_t e = r->GetBinError(i);
+       Double_t o = g->Eval(r->GetBinCenter(i));
+       if (o < 1e-12) { 
+         r->SetBinContent(i, 0);
+         r->SetBinError(i, 0);
+         continue;
+       }
+       r->SetBinContent(i, (c - o) / o + off);
+       r->SetBinError(i, e / o);
       }
+      all->Add(r);
+      pg++;
     }
+    TLegend* leg = StackLegend(all);
+    if (!leg) return;
+      
+    TString   txt      = res->GetTitle();
+    txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
+    if      (ib == 0) txt.Append(" "); // (#times1)");
+    // else if (ib == 1) txt.Append(" (#times10)");
+    else              txt.Append(Form(" (+%d)", off));
 
-    
+    TObject* dummy = 0;
+    TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
+    e->SetMarkerStyle(res->GetMarkerStyle());
+    e->SetMarkerSize(res->GetMarkerSize());
+    e->SetMarkerColor(kBlack);
+    e->SetFillColor(0);
+    e->SetFillStyle(0);
+    e->SetLineColor(kBlack);
   }
+
   /** 
    * Get or create a stack legend.  This is done by adding a TLegend
    * object to the list of functions for the first histogram in the
@@ -905,16 +994,16 @@ struct Unfolder
   TGraphAsymmErrors* GetOther(UShort_t type, Double_t eta, UShort_t sNN,
                              Int_t factor)
   {
-    TString oC = Form("OtherPNch::GetData(%hu,%f,%hu)", 
+    TString oC = Form("OtherPNch::GetData(%hu,%f,%hu);", 
                       type, eta, sNN);
     TGraphAsymmErrors* g = 
       reinterpret_cast<TGraphAsymmErrors*>(gROOT->ProcessLine(oC));
     if (!g) { 
-      Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
-             type, eta, sNN);
+      // Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
+      //         type, eta, sNN);
       return 0;
     }
-
+  
 
     for (Int_t j = 0; j < g->GetN(); j++) { 
       g->SetPoint(j, g->GetX()[j], g->GetY()[j]*factor);
index d90096f991c4628a3dfcd325172f06d6b75f750d..fd26f79ba7abba8d2e8bef0c8d419a857b564da4 100644 (file)
 #pragma link C++ class AliForwardMultiplicityDistribution::Bin+;
 #pragma link C++ class AliForwardMultDists+;
 #pragma link C++ class AliForwardMultDists::EtaBin+;
+#pragma link C++ class AliForwardMultDists::BinSpec+;
 
 #else
 # error Not for compilation