]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Updates to the P(Nch) analysis.
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 Sep 2013 09:12:35 +0000 (09:12 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 Sep 2013 09:12:35 +0000 (09:12 +0000)
- Better summary.
- Separate script with external data
- Summary after unfoldomg
- Create trigger/vertex bias correction

PWGLF/FORWARD/analysis2/AddTaskForwardMultDists.C
PWGLF/FORWARD/analysis2/AliForwardMultDists.cxx
PWGLF/FORWARD/analysis2/AliForwardMultDists.h
PWGLF/FORWARD/analysis2/DrawUnfoldedSummary.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/OtherPNchData.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/SummaryMultDistsDrawer.C
PWGLF/FORWARD/analysis2/scripts/SummaryUnfoldedDrawer.C [new file with mode: 0644]
PWGLF/FORWARD/analysis2/scripts/UnfoldMultDists.C

index 9c851a8708314297235d7de4480a5d56246507e4..44de299e91944e2811006dbadcfaa288274c3d2b 100644 (file)
@@ -28,6 +28,7 @@ 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)
 {
@@ -55,6 +56,8 @@ AddTaskForwardMultDists(const char* trig      = "V0AND",
   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 
index 0395aed693c5a3fdcac858989a147c2b60a496eb..30865266cbac217001b061458ed4a6a3a7f1d947 100644 (file)
@@ -22,6 +22,9 @@ AliForwardMultDists::AliForwardMultDists()
     fList(0),
     fTriggers(0),
     fStatus(0),
+    fVertex(0),
+    fMCVertex(0),
+    fDiag(0),
     fTriggerMask(0),
     fMinIpZ(0),
     fMaxIpZ(0),
@@ -30,6 +33,7 @@ AliForwardMultDists::AliForwardMultDists()
     fCentralCache(0),
     fMCCache(0),
     fMaxN(200),
+    fNDivisions(1),
     fUsePhiAcc(true)
 {}
 
@@ -43,6 +47,9 @@ AliForwardMultDists::AliForwardMultDists(const char* name)
     fList(0),
     fTriggers(0),
     fStatus(0),
+    fVertex(0),
+    fMCVertex(0),
+    fDiag(0),
     fTriggerMask(AliAODForwardMult::kV0AND),
     fMinIpZ(-4),
     fMaxIpZ(4),
@@ -51,6 +58,7 @@ AliForwardMultDists::AliForwardMultDists(const char* name)
     fCentralCache(0),
     fMCCache(0),
     fMaxN(200),
+    fNDivisions(1),
     fUsePhiAcc(true)
 {
   /** 
@@ -72,6 +80,9 @@ AliForwardMultDists::AliForwardMultDists(const AliForwardMultDists& o)
     fList(o.fList),
     fTriggers(o.fTriggers),
     fStatus(o.fStatus),
+    fVertex(o.fVertex),
+    fMCVertex(o.fMCVertex),
+    fDiag(o.fDiag),
     fTriggerMask(o.fTriggerMask),
     fMinIpZ(o.fMinIpZ),
     fMaxIpZ(o.fMaxIpZ),
@@ -80,6 +91,7 @@ AliForwardMultDists::AliForwardMultDists(const AliForwardMultDists& o)
     fCentralCache(o.fCentralCache),
     fMCCache(o.fMCCache),
     fMaxN(o.fMaxN),
+    fNDivisions(o.fNDivisions),
     fUsePhiAcc(o.fUsePhiAcc)
 {}
 //____________________________________________________________________
@@ -94,6 +106,9 @@ AliForwardMultDists::operator=(const AliForwardMultDists& o)
   fList                        = o.fList;
   fTriggers            = o.fTriggers;
   fStatus              = o.fStatus;
+  fVertex               = o.fVertex;
+  fMCVertex             = o.fMCVertex;
+  fDiag                 = o.fDiag;
   fTriggerMask         = o.fTriggerMask;
   fMinIpZ              = o.fMinIpZ;
   fMaxIpZ              = o.fMaxIpZ;
@@ -102,11 +117,22 @@ AliForwardMultDists::operator=(const AliForwardMultDists& o)
   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()
@@ -171,9 +197,71 @@ AliForwardMultDists::UserExec(Option_t* /*option=""*/)
     SetupForData(forwardData, primary != 0);
     fFirstEvent = false;
   }
-    
-  if (!forward->CheckEvent(fTriggerMask, fMinIpZ, fMaxIpZ, 0, 0, 
-                          fTriggers, fStatus)) return;
+
+  Double_t vz  = forward->GetIpZ();
+  Bool_t acc   = forward->CheckEvent(fTriggerMask, fMinIpZ, fMaxIpZ, 0, 0, 
+                                    fTriggers, fStatus);
+  Bool_t trg   = forward->IsTriggerBits(fTriggerMask);
+  Bool_t vtx   = (vz <= fMaxIpZ && vz >= fMinIpZ);
+  Bool_t ok    = true;  // Should bins process this event?
+  Bool_t mcTrg = false;
+  Bool_t mcVtx = false;
+  fVertex->Fill(vz);
+  // If the event was not accepted for analysis 
+  if (primary) {
+    // For MC, we need to check if we should process the event for
+    // MC information or not.
+    if ((fTriggerMask & (AliAODForwardMult::kV0AND|AliAODForwardMult::kNSD))
+       && !forward->IsTriggerBits(AliAODForwardMult::kMCNSD)) 
+      // Bail out if this event is not MC NSD event
+      ok = false;
+      else 
+       mcTrg = true;
+    Double_t mcVz = primary->GetBinContent(0,0);
+    fMCVertex->Fill(mcVz);
+    if (mcVz > fMaxIpZ || mcVz < fMinIpZ) 
+      // Bail out if this event was not in the valid range 
+      ok = false;
+    else 
+      mcVtx = true;
+  }
+#if 0
+  Printf("accepted=%3s triggered=%3s vertex=%3s "
+        "mcTriggered=%3s mcVertex=%3s ok=%3s vz=%f", 
+        (acc ? "yes" : "no"), (trg ? "yes" : "no"), (vtx ? "yes" : "no"), 
+        (mcTrg ? "yes" : "no"), (mcVtx ? "yes" : "no"), (ok ? "yes" : "no"),
+        vz);
+  if (mcTrg && trg) 
+    Printf("Both MC and analysis trigger present");
+#endif
+  if (trg) { 
+    fDiag->Fill(kTrigger, kAnalysis);
+    if (mcTrg)          fDiag->Fill(kTrigger,       kTrigger);
+    if (mcVtx)          fDiag->Fill(kTrigger,       kVertex);
+    if (mcTrg && mcVtx) fDiag->Fill(kTrigger,       kTriggerVertex);
+  }
+  if (vtx) { 
+    fDiag->Fill(kVertex, kAnalysis);
+    if (mcTrg)          fDiag->Fill(kVertex,        kTrigger);
+    if (mcVtx)          fDiag->Fill(kVertex,        kVertex);
+    if (mcTrg && mcVtx) fDiag->Fill(kVertex,        kTriggerVertex);
+  }
+  if (trg && vtx) {
+    fDiag->Fill(kTriggerVertex, kAnalysis);
+    if (mcTrg)          fDiag->Fill(kTriggerVertex, kTrigger);
+    if (mcVtx)          fDiag->Fill(kTriggerVertex, kVertex);
+    if (mcTrg && mcVtx) fDiag->Fill(kTriggerVertex, kTriggerVertex);
+  }
+  if (primary) {
+    if (mcTrg)          fDiag->Fill(kTrigger,       kMC);
+    if (mcVtx)          fDiag->Fill(kVertex,        kMC);
+    if (mcTrg && mcVtx) fDiag->Fill(kTriggerVertex, kMC);
+  }
+  if (!acc && !ok) {
+    // Printf("Event is neither accepted by analysis nor by MC");
+    return; 
+  }
+
   // forward->Print();
   // Info("", "Event vz=%f", forward->GetIpZ());
   
@@ -185,8 +273,8 @@ AliForwardMultDists::UserExec(Option_t* /*option=""*/)
   EtaBin* bin = 0;
   while ((bin = static_cast<EtaBin*>(next()))) {
     bin->Process(*fForwardCache, *fCentralCache,
-                forwardData,    centralData,
-                fMCCache);
+                forwardData,    centralData, 
+                acc,    fMCCache);
   }
 
   PostData(1, fList);
@@ -304,6 +392,7 @@ AliForwardMultDists::Terminate(Option_t* /*option=""*/)
   }
   out->Add(in->FindObject("triggers")->Clone());
   out->Add(in->FindObject("status")->Clone());
+  out->Add(in->FindObject("diagnostics")->Clone());
 
   THStack* sym    = new THStack("all", "All distributions");
   THStack* pos    = new THStack("all", "All distributions");
@@ -320,7 +409,7 @@ AliForwardMultDists::Terminate(Option_t* /*option=""*/)
     else if (bin->IsNegative())  sta = neg;
     else if (bin->IsPositive())  sta = pos;
 
-    TH1*  hh[] = { bin->fSum, bin->fTruth, 0 };
+    TH1*  hh[] = { bin->fSum, bin->fTruth, bin->fTruthAccepted, 0 };
     TH1** ph   = hh;
 
     Int_t idx     = sta->GetHists() ? sta->GetHists()->GetEntries() : 0;
@@ -363,6 +452,7 @@ AliForwardMultDists::StoreInformation(const AliAODForwardMult* aod)
   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)));
 }
 
 //____________________________________________________________________
@@ -404,17 +494,52 @@ AliForwardMultDists::SetupForData(const TH2& hist, Bool_t useMC)
   fCentralCache->GetXaxis()->SetTitle("#eta");
   fCentralCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
   fCentralCache->Sumw2();
+
+  fVertex = new TH1D("vertex", "Vertex distribution",
+                    Int_t(fMaxIpZ-fMinIpZ+.5), 2*fMinIpZ, 2*fMaxIpZ);
+  fVertex->SetDirectory(0);
+  fVertex->SetXTitle("IP_{z} [cm]");
+  fVertex->SetFillColor(kRed+2);
+  fVertex->SetFillStyle(3002);
+  fList->Add(fVertex);
+
   if (useMC) {
     fMCCache->SetDirectory(0);
     fMCCache->GetXaxis()->SetTitle("#eta");
     fMCCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
     fMCCache->Sumw2();
+
+    fMCVertex = static_cast<TH1*>(fVertex->Clone("mcVertex"));
+    fMCVertex->SetTitle("Vertex distribution from MC");
+    fMCVertex->SetDirectory(0);
+    fMCVertex->SetFillColor(kBlue+2);
+    fList->Add(fMCVertex);
   }
 
+  UShort_t xBase = kTrigger-1;
+  UShort_t yBase = kAnalysis-1;
+  fDiag = new TH2I("diagnostics", "Event selection", 
+                  kTriggerVertex-xBase, kTrigger-.5,  kTriggerVertex+.5, 
+                  kTriggerVertex-yBase, kAnalysis-.5, kTriggerVertex+.5);
+  fDiag->SetDirectory(0);
+  fDiag->GetXaxis()->SetTitle("Selection");
+  fDiag->GetXaxis()->SetBinLabel(kTrigger      -xBase, "Trigger");
+  fDiag->GetXaxis()->SetBinLabel(kVertex       -xBase, "Vertex");
+  fDiag->GetXaxis()->SetBinLabel(kTriggerVertex-xBase, "Trigger+Vertex");
+  fDiag->GetYaxis()->SetTitle("Type/MC selection");
+  fDiag->GetYaxis()->SetBinLabel(kAnalysis     -yBase, "Analysis");
+  fDiag->GetYaxis()->SetBinLabel(kMC           -yBase, "MC");
+  fDiag->GetYaxis()->SetBinLabel(kTrigger      -yBase, "MC Trigger");
+  fDiag->GetYaxis()->SetBinLabel(kVertex       -yBase, "MC Vertex");
+  fDiag->GetYaxis()->SetBinLabel(kTriggerVertex-yBase, "MC Trigger+Vertex");
+  fDiag->SetMarkerSize(3);
+  fDiag->SetMarkerColor(kWhite);
+  fList->Add(fDiag);
+
   TIter   next(&fBins);
   EtaBin* bin = 0;
   while ((bin = static_cast<EtaBin*>(next()))) {
-    bin->SetupForData(fList, hist, fMaxN, useMC);
+    bin->SetupForData(fList, hist, fMaxN, fNDivisions, useMC);
   }
     
 }
@@ -552,6 +677,7 @@ AliForwardMultDists::EtaBin::EtaBin()
     fCorr(0),
     fResponse(0), 
     fTruth(0),
+    fTruthAccepted(0),
     fCoverage(0)
 {
   /** 
@@ -570,6 +696,7 @@ AliForwardMultDists::EtaBin::EtaBin(Double_t minEta, Double_t maxEta)
     fCorr(0),
     fResponse(0), 
     fTruth(0),
+    fTruthAccepted(0),
     fCoverage(0)
 {
   /** 
@@ -595,6 +722,7 @@ AliForwardMultDists::EtaBin::EtaBin(const EtaBin& o)
     fCorr(o.fCorr),
     fResponse(o.fResponse), 
     fTruth(o.fTruth),
+    fTruthAccepted(o.fTruthAccepted),
     fCoverage(o.fCoverage)
 {}
 //____________________________________________________________________
@@ -603,16 +731,17 @@ AliForwardMultDists::EtaBin::operator=(const EtaBin& o)
 {
   if (&o == this) return *this;
   
-  fName                = o.fName;
-  fMinEta      = o.fMinEta;
-  fMaxEta      = o.fMaxEta;
-  fMinBin      = o.fMinBin; 
-  fMaxBin      = o.fMaxBin; 
-  fSum         = o.fSum;
-  fCorr                = o.fCorr;
-  fResponse    = o.fResponse; 
-  fTruth       = o.fTruth;
-  fCoverage     = o.fCoverage;
+  fName                 = o.fName;
+  fMinEta       = o.fMinEta;
+  fMaxEta       = o.fMaxEta;
+  fMinBin       = o.fMinBin; 
+  fMaxBin       = o.fMaxBin; 
+  fSum          = o.fSum;
+  fCorr                 = o.fCorr;
+  fResponse     = o.fResponse; 
+  fTruth        = o.fTruth;
+  fTruthAccepted = o.fTruthAccepted;
+  fCoverage      = o.fCoverage;
 
   return *this;
 }
@@ -676,7 +805,8 @@ AliForwardMultDists::EtaBin::FindParent(TList* l, Bool_t create) const
 //____________________________________________________________________
 void
 AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist, 
-                                         UShort_t max, Bool_t useMC)
+                                         UShort_t max, UShort_t nDiv, 
+                                         Bool_t useMC)
 {
   /** 
    * Set-up internal structures on first event. 
@@ -711,11 +841,10 @@ AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
   
   fMinBin = hist.GetXaxis()->FindBin(fMinEta);
   fMaxBin = hist.GetXaxis()->FindBin(fMaxEta-.00001);
-  
-  fSum = new TH1D("rawDist", 
-                 Form("Raw P(#it{N}_{ch}) in %+5.2f<#eta<%+5.2f", 
-                      fMinEta, fMaxEta),
-                 max+1, -.5, max+.5);
+
+  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);
   fSum->SetDirectory(0);
   fSum->GetXaxis()->SetTitle("#it{N}_{ch}");
   fSum->GetYaxis()->SetTitle("Raw P(#it{N}_{ch})");
@@ -728,7 +857,7 @@ AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
   fSum->Sumw2();
   l->Add(fSum);
   
-  fCorr = new TH2D("corr", "Correlation of SPD and FMD signals",
+  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->SetDirectory(0);
   fCorr->GetXaxis()->SetTitle("Forward");
@@ -745,17 +874,18 @@ AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
   l->Add(fCoverage);
 
   if (!useMC) return;
-  fResponse = new TH2D("response", "Reponse matrix", 
-                      max+1, -.5, max+.5, max+1, -.5, max+.5);
+  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->SetDirectory(0);
-  fResponse->SetXTitle("MC");
-  fResponse->SetYTitle("Signal");
+  fResponse->SetYTitle("MC");
+  fResponse->SetXTitle("Signal");
   fResponse->SetOption("colz");
   l->Add(fResponse);
   
-  fTruth = static_cast<TH1D*>(fSum->Clone("truth"));
-  fTruth->SetTitle(Form("True P(#it{N}_{ch}) in %+5.2f<#eta<%+5.2f", 
-                       fMinEta, fMaxEta));
+  fTruth = new TH1D("truth", Form("True P(#it{N}_{ch}) in %s", eta.Data()),
+                   (max+1), -0.5, max+0.5);
+  fTruth->SetXTitle(fSum->GetXaxis()->GetTitle());
+  fTruth->SetYTitle(fSum->GetYaxis()->GetTitle());
   fTruth->SetLineColor(kBlack);
   fTruth->SetFillColor(kBlue+1);
   fTruth->SetFillStyle(0);
@@ -763,8 +893,18 @@ AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
   /// fTruth->SetOption("");
   fTruth->SetMarkerColor(kBlue+1);
   fTruth->SetMarkerStyle(24);
-  // fTruth->Sumw2();
+  fTruth->Sumw2();
   l->Add(fTruth);
+
+  fTruthAccepted = static_cast<TH1D*>(fTruth->Clone("truthAccepted"));
+  fTruthAccepted->SetTitle(Form("True (accepted) P(#it{N}_{ch}) in %s", 
+                               eta.Data()));
+  fTruthAccepted->SetLineColor(kGray+2);
+  fTruthAccepted->SetFillColor(kOrange+2);
+  fTruthAccepted->SetDirectory(0);
+  /// fTruth->SetOption("");
+  fTruthAccepted->SetMarkerColor(kOrange+2);
+  l->Add(fTruthAccepted);
 }
 //____________________________________________________________________
 void
@@ -772,6 +912,7 @@ AliForwardMultDists::EtaBin::Process(const TH1& sumForward,
                                     const TH1& sumCentral,
                                     const TH2& forward,   
                                     const TH2& central,
+                                    Bool_t     accepted, 
                                     const TH1* mc)
 {
   /** 
@@ -782,6 +923,12 @@ AliForwardMultDists::EtaBin::Process(const TH1& sumForward,
    * @param forward     The original forward data 
    * @param central     The original central data
    */
+  if (!mc && !accepted) 
+    // If we're not looking at MC data, and the event wasn't selected,
+    // we bail out - this check is already done, but we do it again
+    // for safety.
+    return;
+
   Double_t sum        = 0;
   Double_t e2sum      = 0;
   Int_t    covered    = 0;
@@ -790,29 +937,26 @@ AliForwardMultDists::EtaBin::Process(const TH1& sumForward,
   Double_t mcsum      = 0;
   Double_t mce2sum    = 0;
   for (Int_t iEta = fMinBin; iEta <= fMaxBin; iEta++) { 
+    if (mc) {
+      Double_t cM = mc->GetBinContent(iEta);
+      Double_t eM = mc->GetBinError(iEta);
+      mcsum   += cM;
+      mce2sum += eM * eM;
+    }
+    if (!accepted) 
+      // Event wasn't selected, but we still need to get the rest of
+      // the MC data.
+      continue;
+
     Double_t cF = sumForward.GetBinContent(iEta);
     Double_t eF = sumForward.GetBinError(iEta);
     Bool_t   bF = forward.GetBinContent(iEta, 0) > 0;
     Double_t cC = sumCentral.GetBinContent(iEta);
     Double_t eC = sumCentral.GetBinError(iEta);
     Bool_t   bC = central.GetBinContent(iEta, 0) > 0;
-    Double_t cM = (mc ? mc->GetBinContent(iEta) : 0);
-    Double_t eM = (mc ? mc->GetBinError(iEta)   : 0);
-
-    /*
-      Info("", 
-      "bF=%d cF=%7.4f+/-%8.5f "
-      "bC=%d cC=%7.4f+/-%8.5f "
-      "cM=%7.4f+/-%48.5f", 
-      bF, cF, eF, bC, cC, eC, cM, eM); 
-    */
-    if (mc) { 
-      mcsum   += cM;
-      mce2sum += eM * eM;
-    }
-
     Double_t c  = 0;
     Double_t e  = 0;
+
     // If we have an overlap - as given by the eta-coverage, 
     // calculate the mean 
     if (bF & bC) { 
@@ -847,23 +991,29 @@ AliForwardMultDists::EtaBin::Process(const TH1& sumForward,
     e2sum += e*e;
   }
       
-  // Fill with weight 
-  Double_t w = 1; // e2sum > 0 ? 1/TMath::Sqrt(e2sum) : 1
-  fSum->Fill(sum, w);
-  fCorr->Fill(fsum, csum);
-
-  Int_t nTotal = fMaxBin - fMinBin + 1;
-  fCoverage->Fill(100*float(covered) / nTotal);
-  // Info(GetName(), "covered: %3d, nTotal: %3d -> %3d%% mc: %3d", 
-  //      covered, nTotal, int(100*float(covered)/nTotal), int(mcsum));
+  if (accepted) {
+    // Only update the histograms if the event was triggered. 
+    // Fill with weight 
+    Double_t w = 1; // e2sum > 0 ? 1/TMath::Sqrt(e2sum) : 1
+    fSum->Fill(sum, w);
+    fCorr->Fill(fsum, csum);
+    
+    Int_t nTotal = fMaxBin - fMinBin + 1;
+    fCoverage->Fill(100*float(covered) / nTotal);
+  }
 
   if (mc) {
+    Double_t w = 1; // mce2sum > 0 ? 1/TMath::Sqrt(mce2sum) : 1
     if (fTruth) {
-      w = 1; // mce2sum > 0 ? 1/TMath::Sqrt(mce2sum) : 1
       fTruth->Fill(mcsum, w);
     }
-    if (fResponse) 
-      fResponse->Fill(mcsum, sum);
+    if (accepted) {
+      if (fResponse) 
+       // Only fill response matrix for accepted events
+       fResponse->Fill(sum, mcsum);
+      if (fTruthAccepted) 
+       fTruthAccepted->Fill(mcsum, w);
+    }
   }
 }
 //____________________________________________________________________
@@ -895,18 +1045,32 @@ AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out, UShort_t maxN)
   o->SetOwner();
   op->Add(o);
 
-  fSum   = static_cast<TH1*>(o->FindObject("rawDist"));
-  fTruth = static_cast<TH1*>(o->FindObject("truth"));
-  TH1*  hists[] = { fSum, fTruth, 0 };
+  fSum           = static_cast<TH1*>(o->FindObject("rawDist"));
+  fTruth         = static_cast<TH1*>(o->FindObject("truth"));
+  fTruthAccepted = static_cast<TH1*>(o->FindObject("truthAccepted"));
+
+  TH1*  hists[] = { fSum, fTruth, fTruthAccepted, 0 };
   TH1** phist   = hists;
   while (*phist) { 
     TH1* h = *phist;
     if (h) { 
       Double_t intg = h->Integral(1, maxN);
-      h->Scale(1. / intg);
+      h->Scale(1. / intg, "width");
     }
     phist++;
   }
+
+  if (fTruth && fTruthAccepted) {
+    TH1*  trgVtx  = static_cast<TH1*>(fTruthAccepted->Clone("triggerVertex"));
+    TString tit(fTruth->GetTitle());
+    tit.ReplaceAll("True P(#it{N}_{ch})", "C_{trigger,vertex}");
+    trgVtx->SetTitle(tit);
+    trgVtx->SetYTitle("P_{MC}(#it{N}_{ch})/P_{MC,accepted}(#it{N}_{ch})");
+    trgVtx->Divide(fTruth);
+    trgVtx->SetDirectory(0);
+    o->Add(trgVtx);
+  }
+  
 }
 //====================================================================
 // 
index 7a806ed10e3f9116c14330b5bf0d52fc860288eb..f5e6517b2c04bbfba8e11e491a8c2cb015f83513 100644 (file)
@@ -17,6 +17,13 @@ public:
   enum { 
     kInvalidEta = 999
   };
+  enum { 
+    kAnalysis = 1, 
+    kMC       = 2, 
+    kTrigger  = 3, 
+    kVertex   = 4, 
+    kTriggerVertex = 5
+  };
   /** 
    * Default constructor
    */
@@ -102,8 +109,10 @@ public:
    * Set the maximum @f$N_{ch}@f$ to consider 
    * 
    * @param maxN Maximum
+   * @param nBins (Optionally) the number of bins per particle number 
    */
-  void SetMaxN(UInt_t maxN) { fMaxN = maxN; }
+  void SetMaxN(UShort_t maxN, UShort_t nBins=-1);
+  void SetNDivisions(UInt_t nDiv) { fNDivisions = nDiv; }
   /** 
    * Set the range of valid interaction points 
    * 
@@ -208,9 +217,11 @@ public:
      * @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, Bool_t useMC);
+    void SetupForData(TList* list, const TH2& hist, UShort_t max, 
+                     UShort_t nDiv, Bool_t useMC);
     /** 
      * Process a single event 
      * 
@@ -218,10 +229,12 @@ public:
      * @param sumCentral  Projection of the central data
      * @param forward     The original forward data 
      * @param central     The original central data
+     * @param accepted    True if event is accepted for analysis
+     * @param mc          Distribution of primary particles from MC
      */
     void Process(const TH1& sumForward, const TH1& sumCentral,
                 const TH2& forward,    const TH2& central,
-                const TH1* mc);
+                Bool_t     accepted,   const TH1* mc);
     /** 
      * Called at the end of the final processing of the job on the
      * full data set (merged data)
@@ -232,16 +245,17 @@ public:
      */
     void Terminate(TList* in, TList* out, UShort_t maxN);
       
-    TString  fName;     // Name of this bin
-    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
-    Int_t    fMaxBin;   // Largest @f$\eta@f$ bin to consider
-    TH1*     fSum;      // Distribution 
-    TH2*     fCorr;     // Correlation between forward and central
-    TH2*     fResponse; // Response matrix (for MC)
-    TH1*     fTruth;    // `true' distribution 
-    TH1*     fCoverage; // How much was covered
+    TString  fName;          // Name of this bin
+    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
+    Int_t    fMaxBin;        // Largest @f$\eta@f$ bin to consider
+    TH1*     fSum;           // Distribution 
+    TH2*     fCorr;          // Correlation between forward and central
+    TH2*     fResponse;      // Response matrix (for MC)
+    TH1*     fTruth;         // `true' distribution 
+    TH1*     fTruthAccepted; // `true' distribution for accepted events
+    TH1*     fCoverage;      // How much was covered
 
     ClassDef(EtaBin,1);
   };
@@ -252,6 +266,9 @@ public:
   TList*   fList;         // Output 
   TH1*     fTriggers;     // Histogram of triggers
   TH1*     fStatus;       // Histogram of event selection status 
+  TH1*     fVertex;       // Histogram of IpZ
+  TH1*     fMCVertex;     // Histogram of MC IpZ
+  TH2*     fDiag;         // Diagnostics
   UInt_t   fTriggerMask;  // Trigger mask
   Double_t fMinIpZ;       // Least @f$IP_{z}@f$ to consider 
   Double_t fMaxIpZ;       // Largest @f$IP_{z}@f$ to consider 
@@ -259,7 +276,8 @@ public:
   TH1*     fForwardCache; // Projection cache 
   TH1*     fCentralCache; // Projection cache 
   TH1*     fMCCache;      // Projection cache 
-  UInt_t   fMaxN;         // Maximum of @f$N_{ch}@f$ 
+  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);
diff --git a/PWGLF/FORWARD/analysis2/DrawUnfoldedSummary.C b/PWGLF/FORWARD/analysis2/DrawUnfoldedSummary.C
new file mode 100644 (file)
index 0000000..a5ecfe4
--- /dev/null
@@ -0,0 +1,23 @@
+/**
+ * @file   DrawAODSummary.C
+ * @author Christian Holm Christensen <cholm@master.hehi.nbi.dk>
+ * @date   Tue Oct 30 09:47:30 2012
+ * 
+ * @brief  Script to draw summary of AOD pass into a PDF 
+ * 
+ * 
+ */
+//____________________________________________________________________
+void DrawUnfoldedSummary(const char* fname="forward_unfolded.root")
+{
+  gROOT->SetMacroPath(Form("%s:$ALICE_ROOT/../trunk/PWGLF/FORWARD/analysis2/scripts",
+                          gROOT->GetMacroPath()));
+  gROOT->LoadMacro("SummaryUnfoldedDrawer.C++g");
+  
+  SummaryUnfoldedDrawer d;
+  d.Run(fname);
+}
+
+//
+// EOF
+//
diff --git a/PWGLF/FORWARD/analysis2/scripts/OtherPNchData.C b/PWGLF/FORWARD/analysis2/scripts/OtherPNchData.C
new file mode 100644 (file)
index 0000000..51a1f3d
--- /dev/null
@@ -0,0 +1,1326 @@
+#include <TGraphAsymmErrors.h>
+
+struct OtherPNch 
+{
+  /* ===================================================================
+   * 
+   * CMS data 
+   */
+  /**
+   * CMS @f$ P(N_{ch})@f$ at @f$\sqrt{s}=900GeV@f$ for @f$|\eta|\lt0.5@f$
+   * http://hepdata.cedar.ac.uk/view/p8068/d2
+   *
+   * @return p8068_d2x1y1
+   */
+  TGraphAsymmErrors* CMSsqrts900eta05()
+  {
+    double x[] = { 0.0,        1.0, 
+                  2.0, 3.0,
+                  4.0, 5.0,
+                  6.0, 7.0,
+                  8.0, 9.0,
+                  10.0,        11.0,
+                  12.0,        13.0,
+                  14.0,        15.0,
+                  16.0,        17.0,
+                  18.0,        19.0,
+                  20.0,        22.0,
+                  26.0 };
+    double xem[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.5,
+                    2.5 };
+    double xep[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.5,
+                    2.5 };
+    double y[] = { 0.193,      0.1504,
+                  0.13507,     0.11647,
+                  0.09566,     0.07567,
+                  0.05817,     0.04427,
+                  0.0337,      0.02572,
+                  0.01948,     0.01457,
+                  0.01069,     0.00769,
+                  0.0055,      0.00396,
+                  0.00289,     0.002112,
+                  0.001544,    0.001118,
+                  7.52E-4,     3.35E-4,
+                  9.8E-5 };
+    double yem[] = { 0.02189155088156159,      0.012180722474467595,
+                    0.007909551188278638,      0.006072832946821442,
+                    0.004838605170914444,      0.003828276374558138,
+                    0.003016836091006603,      0.002365797962633327,
+                    0.0018847015678881366,     0.0015610893632332517,
+                    0.001308166656049603,      0.0011140017953306899,
+                    9.219544457292888E-4,      7.203471385380802E-4,
+                    5.508175741568165E-4,      4.3011626335213136E-4,
+                    3.3837848631377264E-4,     2.746961958236772E-4,
+                    2.2377220560203628E-4,     1.9070920271449933E-4,
+                    1.6011246047700347E-4,     8.134494452638099E-5,
+                    2.5059928172283337E-5 };
+    double yep[] = { 0.031663227883461285,     0.015861904047118684,
+                    0.010948538715280683,      0.008490583018850945,
+                    0.007046453008429135,      0.005675746646917919,
+                    0.004395201929377079,      0.003365070578754627,
+                    0.0026049952015310893,     0.0020535335400231475,
+                    0.0016522106403240478,     0.0013497407158413798,
+                    0.001118033988749895,      9.060905032059435E-4,
+                    7.061161377563892E-4,      5.457105459856901E-4,
+                    4.148493702538308E-4,      3.257759966602819E-4,
+                    2.6689698387205504E-4,     2.292705825002414E-4,
+                    1.904205871222962E-4,      9.646242791885347E-5,
+                    3.417601498127012E-5 };
+    int np = 23;
+    TGraphAsymmErrors* g = new TGraphAsymmErrors(np, 
+                                                x, y, 
+                                                xem, 
+                                                xep, 
+                                                yem, 
+                                                yep);
+    g->SetName("/HepData/8068/d2x1y1");
+    g->SetTitle("/HepData/8068/d2x1y1");
+
+    return g;
+  }
+  /**
+   * CMS @f$ P(N_{ch})@f$ at @f$\sqrt{s}=900GeV@f$ for @f$|\eta|\lt1.0@f$
+   * http://hepdata.cedar.ac.uk/view/p8068/d3
+   *
+   * @return p8068_d3x1y1
+   */
+  TGraphAsymmErrors* CMSsqrts900eta10()
+  {
+
+    double x[] = { 0.0,        1.0,
+                  2.0, 3.0,
+                  4.0, 5.0,
+                  6.0, 7.0,
+                  8.0, 9.0,
+                  10.0,        11.0,
+                  12.0,        13.0,
+                  14.0,        15.0,
+                  16.0,        17.0,
+                  18.0,        19.0,
+                  20.0,        21.0,
+                  22.0,        23.0,
+                  24.0,        25.0,
+                  26.0,        27.0,
+                  28.0,        29.0,
+                  30.0,        31.0,
+                  32.0,        33.0,
+                  34.0,        36.0,
+                  40.0,        45.0,
+                  50.0,        55.0 };
+    double xem[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.5,
+                    2.5,       2.5,
+                    2.5,       2.5 };
+    double xep[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.5,
+                    2.5,       2.5,
+                    2.5,       2.5 };
+    double y[] = { 0.1076,     0.0655,
+                  0.07574,     0.08168,
+                  0.0813,      0.07622,
+                  0.06843,     0.06034,
+                  0.05259,     0.0458,
+                  0.03987,     0.03472,
+                  0.03035,     0.02661,
+                  0.02328,     0.02018,
+                  0.01728,     0.0147,
+                  0.0125,      0.01067,
+                  0.00916,     0.00786,
+                  0.00671,     0.0057,
+                  0.00478,     0.00394,
+                  0.00319,     0.00251,
+                  0.00196,     0.00156,
+                  0.001301,    0.001121,
+                  9.84E-4,     8.25E-4,
+                  6.39E-4,     3.78E-4,
+                  1.96E-4,     6.0E-5,
+                  1.25E-5,     1.6E-6 };
+    double yem[] = { 0.015802847844613326,     0.009464142856064674,
+                    0.008214195030555337,      0.00676426640516176,
+                    0.0058916636020736966,     0.005200240379059414,
+                    0.004441733895676326,      0.003671511950137164,
+                    0.0029149442533262966,     0.0024150569351466646,
+                    0.0020554318281081475,     0.001767059704707229,
+                    0.0015500322577288513,     0.0013800362314084368,
+                    0.0012684636376341263,     0.0011763077828527702,
+                    0.001084158659975559,      9.82344135219425E-4,
+                    8.805679985100526E-4,      7.692203845452875E-4,
+                    6.705221845696084E-4,      5.913543776789007E-4,
+                    5.375872022286245E-4,      5.154609587543949E-4,
+                    5.124451190127582E-4,      4.904079934095691E-4,
+                    4.4922154890432404E-4,     3.7E-4,
+                    3.0083217912982643E-4,     2.4166091947189146E-4,
+                    2.1108529081866412E-4,     1.8252944967867517E-4,
+                    1.662077013859466E-4,      1.4663219291819924E-4,
+                    1.2229881438509532E-4,     7.826237921249263E-5,
+                    4.74236228055175E-5,       2.1095023109728988E-5,
+                    6.161980201201558E-6,      1.3601470508735443E-6 };
+    double yep[] = { 0.02725949375905576,      0.011154371340420759,
+                    0.009740395269186974,      0.008248078564126314,
+                    0.007285581651453781,      0.006454339625399333,
+                    0.005595158621522718,      0.004804039966528172,
+                    0.00405504623894722,       0.0034544174617437305,
+                    0.0029644561052577585,     0.002545584412271571,
+                    0.0021980445855350615,     0.0019087430418995638,
+                    0.0016690416411821486,     0.001468911161370898,
+                    0.0013179529581893276,     0.0011866338946785568,
+                    0.001065129100156408,      9.339164844888434E-4,
+                    8.052328855678958E-4,      6.868041933477111E-4,
+                    5.94810894318522E-4,       5.440588203494177E-4,
+                    5.220153254455275E-4,      5.192301994298867E-4,
+                    4.972926703662542E-4,      4.368065933568311E-4,
+                    3.5735136770411276E-4,     2.8792360097775935E-4,
+                    2.4200206610688264E-4,     2.0310588371585893E-4,
+                    1.670748335327616E-4,      1.364587849865299E-4,
+                    1.1970797801316335E-4,     9.470480452437458E-5,
+                    5.758472019555187E-5,      2.195449840010015E-5,
+                    7.481310045707236E-6,      1.7029386365926402E-6 };
+    int np = 40;
+    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
+                                                x,
+                                                y,
+                                                xem,
+                                                xep,
+                                                yem,
+                                                yep);
+    g->SetName("/HepData/8068/d3x1y1");
+    g->SetTitle("/HepData/8068/d3x1y1");
+    return g;
+  }
+  /**
+   * CMS @f$ P(N_{ch})@f$ at @f$\sqrt{s}=900GeV@f$ for @f$|\eta|\lt1.5@f$
+   * http://hepdata.cedar.ac.uk/view/p8068/d4
+   *
+   * @return p8068_d4x1y1
+   */
+  TGraphAsymmErrors* CMSsqrts900eta15()
+  {
+    double x[] = { 0.0,        1.0,
+                  2.0, 3.0,
+                  4.0, 5.0,
+                  6.0, 7.0,
+                  8.0, 9.0,
+                  10.0,        11.0,
+                  12.0,        13.0,
+                  14.0,        15.0,
+                  16.0,        17.0,
+                  18.0,        19.0,
+                  20.0,        21.0,
+                  22.0,        23.0,
+                  24.0,        25.0,
+                  26.0,        27.0,
+                  28.0,        29.0,
+                  30.0,        31.0,
+                  32.0,        33.0,
+                  34.0,        35.0,
+                  36.0,        37.0,
+                  38.0,        39.0,
+                  40.0,        41.0,
+                  42.0,        43.0,
+                  44.0,        45.0,
+                  46.5,        49.0,
+                  53.0,        58.0,
+                  63.0,        68.0 };
+    double xem[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    1.0,       1.5,
+                    2.5,       2.5,
+                    2.5,       2.5 };
+    double xep[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    1.0,       1.5,
+                    2.5,       2.5,
+                    2.5,       2.5 };
+    double y[] = { 0.0753,     0.03557,
+                  0.04347,     0.05068,
+                  0.05566,     0.05786,
+                  0.05835,     0.05688,
+                  0.05348,     0.0489,
+                  0.04422,     0.04002,
+                  0.03647,     0.03337,
+                  0.03055,     0.02797,
+                  0.02556,     0.02329,
+                  0.02122,     0.01936,
+                  0.01765,     0.01604,
+                  0.01451,     0.01304,
+                  0.01168,     0.01045,
+                  0.00937,     0.0084,
+                  0.00751,     0.00667,
+                  0.00592,     0.00526,
+                  0.00466,     0.00413,
+                  0.00362,     0.00316,
+                  0.00272,     0.00233,
+                  0.00203,     0.00178,
+                  0.00156,     0.00136,
+                  0.0012,      0.001053,
+                  9.13E-4,     7.84E-4,
+                  5.85E-4,     3.73E-4,
+                  2.13E-4,     9.6E-5,
+                  3.6E-5,      7.4E-6 };
+    double yem[] = { 0.014524806367039803,     0.0061507154055443014,
+                    0.006900384047283165,      0.0068944397886992964,
+                    0.00601252858621063,       0.004686747272896204,
+                    0.0042880881520789655,     0.0039679339712248235,
+                    0.0035288099977187778,     0.0031197756329582422,
+                    0.0027127108213003464,     0.0023229722340140013,
+                    0.002033937068839643,      0.0017866449003649271,
+                    0.001591257364476281,      0.0014504137340772805,
+                    0.0013412307780542466,     0.0012490796611905903,
+                    0.0011594826432508594,     0.0010893117092916976,
+                    0.0010384603988597735,     9.972462083156796E-4,
+                    9.588013350011566E-4,      8.984430978086481E-4,
+                    8.160882305241266E-4,      7.178439941937245E-4,
+                    6.390618123468182E-4,      5.88727441181401E-4,
+                    5.571355310873647E-4,      5.255473337388365E-4,
+                    4.846648326421054E-4,      4.6615448083226656E-4,
+                    4.4384682042344296E-4,     4.3081318457076036E-4,
+                    4.0853396431630995E-4,     3.992492955535426E-4,
+                    3.676955262170047E-4,      3.2695565448543625E-4,
+                    2.9068883707497264E-4,     2.505992817228334E-4,
+                    2.1954498400100151E-4,     1.8601075237738274E-4,
+                    1.8601075237738274E-4,     1.7042300314218148E-4,
+                    1.7998055450520203E-4,     1.6058953888718904E-4,
+                    1.4985659811966905E-4,     8.561541917201597E-5,
+                    6.324555320336759E-5,      4.657252408878007E-5,
+                    2.4166091947189146E-5,     5.590169943749474E-6 };
+    double yep[] = { 0.03405304685340212,      0.00724989655098609,
+                    0.007817985674072318,      0.007751728839426725,
+                    0.006919566460407762,      0.00571192611997039,
+                    0.005232991496266739,      0.004813003220443552,
+                    0.00435332057170156,       0.0038838769290491172,
+                    0.003445692963686695,      0.0030549959083442323,
+                    0.0027252339349127445,     0.0024367396249907374,
+                    0.002189794510907359,      0.001969568480657629,
+                    0.0017906702655709678,     0.0016395731151735808,
+                    0.0015002999700059986,     0.0013710215169719256,
+                    0.0012712198865656563,     0.0011812705024675763,
+                    0.0011035397591387453,     0.0010139033484509261,
+                    9.121403400793104E-4,      8.132035415564789E-4,
+                    7.242237223399962E-4,      6.545991139621257E-4,
+                    6.135144660071188E-4,      5.724508712544685E-4,
+                    5.692099788303082E-4,      5.126402247190518E-4,
+                    4.531004303683677E-4,      4.4011362169330773E-4,
+                    4.1785164831552356E-4,     4.0853396431630995E-4,
+                    3.8626415831655934E-4,     3.54682957019364E-4,
+                    3.1780497164141406E-4,     2.7730849247724094E-4,
+                    2.459674775249769E-4,      2.3706539182259396E-4,
+                    2.1954498400100151E-4,     2.0068881383873889E-4,
+                    1.5864425612041553E-4,     1.6058953888718904E-4,
+                    1.4338758663147937E-4,     1.0245974819410792E-4,
+                    6.609841147864296E-5,      4.657252408878007E-5,
+                    2.5079872407968906E-5,     5.941380311005179E-6 };
+    int np = 52;
+    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
+                                                x,
+                                                y,
+                                                xem,
+                                                xep,
+                                                yem,
+                                                yep);
+    g->SetName("/HepData/8068/d4x1y1");
+    g->SetTitle("/HepData/8068/d4x1y1");
+    return g;
+  }
+  /**
+   * CMS @f$ P(N_{ch})@f$ at @f$\sqrt{s}=900GeV@f$ for @f$|\eta|\lt2.4@f$
+   * http://hepdata.cedar.ac.uk/view/p8068/d5
+   *
+   * @return p8068_d5x1y1
+   */
+  TGraphAsymmErrors* CMSsqrts900eta24()
+  {
+    double x[] = { 0.0,        1.0,
+                  2.0, 3.0,
+                  4.0, 5.0,
+                  6.0, 7.0,
+                  8.0, 9.0,
+                  10.0,        11.0,
+                  12.0,        13.0,
+                  14.0,        15.0,
+                  16.0,        17.0,
+                  18.0,        19.0,
+                  20.0,        21.0,
+                  22.0,        23.0,
+                  24.0,        25.0,
+                  26.0,        27.0,
+                  28.0,        29.0,
+                  30.0,        31.0,
+                  32.0,        33.0,
+                  34.0,        35.0,
+                  36.0,        37.0,
+                  38.0,        39.0,
+                  40.0,        41.0,
+                  42.0,        43.0,
+                  44.0,        45.0,
+                  46.0,        47.0,
+                  48.0,        49.0,
+                  50.0,        51.0,
+                  52.0,        53.0,
+                  54.0,        55.0,
+                  56.0,        57.5,
+                  60.0,        66.5,
+                  76.5,        86.5 };
+    double xem[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.0,
+                    1.5,       5.0,
+                    5.0,       5.0 };
+    double xep[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.0,
+                    1.5,       5.0,
+                    5.0,       5.0 };
+    double y[] = { 0.0561,     0.02321,
+                  0.02645,     0.031,
+                  0.0362,      0.04081,
+                  0.04402,     0.04548,
+                  0.04592,     0.04515,
+                  0.04342,     0.04101,
+                  0.03836,     0.03573,
+                  0.03316,     0.03079,
+                  0.02868,     0.02679,
+                  0.02506,     0.02343,
+                  0.02187,     0.02039,
+                  0.01901,     0.0178,
+                  0.01668,     0.01556,
+                  0.01445,     0.01334,
+                  0.01229,     0.01132,
+                  0.01045,     0.00967,
+                  0.00896,     0.00829,
+                  0.00763,     0.00699,
+                  0.00637,     0.00578,
+                  0.00527,     0.00482,
+                  0.00441,     0.00412,
+                  0.00363,     0.00328,
+                  0.00292,     0.00261,
+                  0.00233,     0.00209,
+                  0.00185,     0.00164,
+                  0.00146,     0.00131,
+                  0.00117,     0.00106,
+                  9.7E-4,      8.9E-4,
+                  8.2E-4,      6.63E-4,
+                  4.53E-4,     2.12E-4,
+                  7.7E-5,      1.12E-5 };
+    double yem[] = { 0.01767766952966369,      0.003660737630587584,
+                    0.005079094407470687,      0.006069810540700591,
+                    0.006147235150862541,      0.005458580401532985,
+                    0.004473712105176193,      0.0035970126494078384,
+                    0.003456848275524976,      0.003228761372415125,
+                    0.0029797483115189447,     0.002739361239413305,
+                    0.0025401181074902798,     0.002313114783144148,
+                    0.0020573040611440983,     0.001820109886792553,
+                    0.0016324827717314506,     0.001484318025222358,
+                    0.0013727709204379296,     0.001287866452703851,
+                    0.001217538500417954,      0.0011279184367674819,
+                    0.0010412012293500234,     9.741663102366044E-4,
+                    9.360555539069249E-4,      9.139474820797965E-4,
+                    9.013878188659973E-4,      8.697700845625814E-4,
+                    8.190848551890091E-4,      7.433034373659253E-4,
+                    6.741661516273269E-4,      6.239390995922599E-4,
+                    5.961543424315552E-4,      5.8309518948453E-4,
+                    5.738466694161429E-4,      5.515432893255072E-4,
+                    5.292447448959697E-4,      4.976946855251721E-4,
+                    4.609772228646444E-4,      4.295346318982906E-4,
+                    4.0249223594996216E-4,     4.386342439892262E-4,
+                    3.712142238654117E-4,      3.6674241641784496E-4,
+                    3.6249137920783716E-4,     3.5341194094144583E-4,
+                    3.3105890714493693E-4,     3.130495168499705E-4,
+                    2.8178005607210744E-4,     2.5553864678361276E-4,
+                    2.3323807579381201E-4,     2.1633307652783938E-4,
+                    1.8601075237738274E-4,     1.6401219466856724E-4,
+                    1.562049935181331E-4,      1.562049935181331E-4,
+                    1.5480633061990717E-4,     1.395886814895821E-4,
+                    1.0756393447619885E-4,     6.676076692189808E-5,
+                    2.9966648127543395E-5,     6.551335741663679E-6 };
+    double yep[] = { 0.09323352401362935,      0.004710169848317574,
+                    0.005736898116578331,      0.006618164398078972,
+                    0.006715772479767314,      0.006096630216767292,
+                    0.00519042387479096,       0.004362247127341595,
+                    0.004172217635742412,      0.003913693907295255,
+                    0.003644283194264683,      0.0033638668225719043,
+                    0.0031145144083789367,     0.0028765604460883488,
+                    0.002629087294100369,      0.0023903974564912843,
+                    0.002181604913819182,      0.0020026232796010335,
+                    0.0018514858897652987,     0.0017280335644888382,
+                    0.00160822883943797,       0.0014787156589419076,
+                    0.0013612494260788505,     0.0012539936203984452,
+                    0.0011763077828527702,     0.0011157060544785082,
+                    0.001084158659975559,      0.001033247308247159,
+                    9.630160954002794E-4,      9.244457799135652E-4,
+                    8.35224520712844E-4,       6.989277502002622E-4,
+                    6.519202405202648E-4,      6.203224967708329E-4,
+                    6.10982814815605E-4,       5.88727441181401E-4,
+                    5.664803615307418E-4,      5.348831648126533E-4,
+                    4.976946855251721E-4,      4.477722635447622E-4,
+                    4.204759208325728E-4,      3.9357337308308855E-4,
+                    3.981205847478877E-4,      3.8483762809787716E-4,
+                    3.807886552931954E-4,      3.716180835212409E-4,
+                    3.584689665786984E-4,      3.4928498393145964E-4,
+                    3.1780497164141406E-4,     2.9068883707497264E-4,
+                    2.505992817228334E-4,      2.2472205054244234E-4,
+                    2.0248456731316588E-4,     1.8027756377319944E-4,
+                    1.6401219466856724E-4,     1.6401219466856724E-4,
+                    1.5640012787718557E-4,     1.4223923509355636E-4,
+                    1.0288342918079667E-4,     6.86804193347711E-5,
+                    3.361547262794322E-5,      7.071067811865475E-6 };
+    int np = 62;
+    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
+                                                x,
+                                                y,
+                                                xem,
+                                                xep,
+                                                yem,
+                                                yep);
+    g->SetName("/HepData/8068/d5x1y1");
+    g->SetTitle("/HepData/8068/d5x1y1");
+    return g;
+  }
+  /**
+   * CMS @f$ P(N_{ch})@f$ at @f$\sqrt{s}=900GeV@f$ for @f$|\eta|\lt2.0@f$
+   * http://hepdata.cedar.ac.uk/view/p8068/d6
+   *
+   * @return p8068_d6x1y1
+   */
+  TGraphAsymmErrors* CMSsqrts900eta20()
+  {
+    double x[] = { 0.0,        1.0,
+                  2.0, 3.0,
+                  4.0, 5.0,
+                  6.0, 7.0,
+                  8.0, 9.0,
+                  10.0,        11.0,
+                  12.0,        13.0,
+                  14.0,        15.0,
+                  16.0,        17.0,
+                  18.0,        19.0,
+                  20.0,        21.0,
+                  22.0,        23.0,
+                  24.0,        25.0,
+                  26.0,        27.0,
+                  28.0,        29.0,
+                  30.0,        31.0,
+                  32.0,        33.0,
+                  34.0,        35.0,
+                  36.0,        37.0,
+                  38.0,        39.0,
+                  40.0,        41.0,
+                  42.0,        43.0,
+                  44.0,        45.0,
+                  46.0,        47.0,
+                  48.0,        49.0,
+                  50.0,        51.0,
+                  52.0,        53.0,
+                  54.0,        55.0,
+                  56.0,        57.0,
+                  58.0,        59.0,
+                  60.0,        61.0,
+                  62.5,        64.5,
+                  67.5,        74.5,
+                  84.5,        94.5 };
+    double xem[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    1.0,       1.0,
+                    2.0,       5.0,
+                    5.0,       5.0 };
+    double xep[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    1.0,       1.0,
+                    2.0,       5.0,
+                    5.0,       5.0 };
+    double y[] = { 0.0494,     0.01789,
+                  0.01895,     0.02174,
+                  0.0255,      0.02961,
+                  0.0335,      0.03646,
+                  0.03826,     0.03906,
+                  0.0392,      0.03868,
+                  0.03753,     0.03596,
+                  0.0341,      0.03206,
+                  0.03,        0.02802,
+                  0.02616,     0.02448,
+                  0.023,       0.0217,
+                  0.02058,     0.01959,
+                  0.01866,     0.01776,
+                  0.01685,     0.01595,
+                  0.01504,     0.01415,
+                  0.01326,     0.01241,
+                  0.01161,     0.01087,
+                  0.01015,     0.00948,
+                  0.00887,     0.00829,
+                  0.00777,     0.00734,
+                  0.00681,     0.00624,
+                  0.00578,     0.00535,
+                  0.00494,     0.00457,
+                  0.00422,     0.00388,
+                  0.00356,     0.00326,
+                  0.00298,     0.00274,
+                  0.00253,     0.00234,
+                  0.00212,     0.00189,
+                  0.00167,     0.0015,
+                  0.00136,     0.00123,
+                  0.00112,     0.001041,
+                  9.17E-4,     7.64E-4,
+                  4.92E-4,     2.45E-4,
+                  9.9E-5,      1.64E-5 };
+    double yem[] = { 0.031045933711196384,     0.0025678006153126453,
+                    0.003808017857100988,      0.004928792955683978,
+                    0.005475372133471842,      0.005444713031923721,
+                    0.004928792955683978,      0.0041622710147226115,
+                    0.003405994715204356,      0.0030985964564621835,
+                    0.0029811742652854096,     0.002811547616527239,
+                    0.0026404734423962684,     0.002492308167141455,
+                    0.0023525518060183073,     0.0021654098919142305,
+                    0.001967053634245899,      0.0017886866690396059,
+                    0.0016107451691686056,     0.0014552319402761885,
+                    0.001326989073052224,      0.001225275479229059,
+                    0.001164817582284883,      0.0011168258592994702,
+                    0.0010756393447619887,     0.001027861858422619,
+                    9.7082439194738E-4,        9.202716990106781E-4,
+                    8.825531145489205E-4,      8.509406559801923E-4,
+                    8.228000972289684E-4,      8.099382692526635E-4,
+                    7.912016177940992E-4,      7.725283166331187E-4,
+                    7.316419889536139E-4,      6.815423684555495E-4,
+                    6.407807737440318E-4,      5.950630218724736E-4,
+                    5.636488268416782E-4,      5.636488268416782E-4,
+                    5.322593352868506E-4,      4.919349550499538E-4,
+                    4.7853944456021594E-4,     4.7423622805517506E-4,
+                    4.651881339845203E-4,      4.428317965096906E-4,
+                    4.248529157249601E-4,      3.9824615503479755E-4,
+                    3.671511950137164E-4,      3.3615472627943223E-4,
+                    3.138470965295043E-4,      3.138470965295043E-4,
+                    3.176476034853718E-4,      3.3541019662496844E-4,
+                    3.3105890714493693E-4,     2.996664812754339E-4,
+                    2.641968962724581E-4,      2.3853720883753126E-4,
+                    2.1633307652783938E-4,     1.9416487838947602E-4,
+                    1.6401219466856724E-4,     1.5181897114655994E-4,
+                    1.5966527487215246E-4,     1.3961733416735901E-4,
+                    9.740636529508736E-5,      6.800735254367721E-5,
+                    3.733630940518893E-5,      9.437160589923222E-6 };
+    double yep[] = { 0.027627884464793896,     0.0038148263394288343,
+                    0.004395600072800072,      0.005347317084295637,
+                    0.005874325493194942,      0.0058935897380119695,
+                    0.005446999173857106,      0.004729587719875803,
+                    0.004031935019317648,      0.0037237615390892046,
+                    0.0035759474269066098,     0.00339607125955861,
+                    0.0032050585018061684,     0.003026549190084311,
+                    0.0028468403537957655,     0.00266865134478073,
+                    0.002469412885687608,      0.0022703523955544874,
+                    0.0020813697413001857,     0.0019043371550227129,
+                    0.00175524927004685,       0.001643471934655411,
+                    0.0015533190271158081,     0.001475127113166862,
+                    0.0014046351839534703,     0.0013267252918370102,
+                    0.0012490796611905903,     0.0011691449867317568,
+                    0.0011208925015361642,     0.0011472575996697516,
+                    9.932774033471214E-4,      9.52102935611481E-4,
+                    9.139474820797965E-4,      8.75956619930462E-4,
+                    8.34865258589672E-4,       7.749193506423748E-4,
+                    7.244998274671983E-4,      6.685057965343307E-4,
+                    6.092618484691126E-4,      5.727128425310541E-4,
+                    5.412947441089743E-4,      5.280151512977634E-4,
+                    5.147815070493501E-4,      5.015974481593781E-4,
+                    4.924428900898052E-4,      4.792702786528704E-4,
+                    4.7010637094172637E-4,     4.428317965096906E-4,
+                    4.204759208325728E-4,      3.8910152916687437E-4,
+                    3.6674241641784496E-4,     3.488552708502481E-4,
+                    3.2649655434629013E-4,     3.0886890422961E-4,
+                    3.04138126514911E-4,       2.996664812754339E-4,
+                    2.9068883707497264E-4,     2.641968962724581E-4,
+                    2.418677324489565E-4,      2.1095023109728988E-4,
+                    1.972308292331602E-4,      1.739482681718907E-4,
+                    1.4475496537252186E-4,     1.2083045973594572E-4,
+                    1.097679370308106E-4,      6.992138442565336E-5,
+                    3.54682957019364E-5,       9.078546139112804E-6 };
+    int np = 68;
+    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
+                                                x,
+                                                y,
+                                                xem,
+                                                xep,
+                                                yem,
+                                                yep);
+    g->SetName("/HepData/8068/d6x1y1");
+    g->SetTitle("/HepData/8068/d6x1y1");
+    return   g;
+  }
+  /** 
+   * Get CMS data 
+   * 
+   * @param eta Limits on eta 
+   * @param sNN Center of mass energy 
+   * 
+   * @return Graph or null
+   */
+  TGraphAsymmErrors* GetCMS(Double_t eta, UShort_t sNN)
+  {
+    Double_t aEta = TMath::Abs(eta);
+    TGraphAsymmErrors* g = 0;
+    switch (sNN) { 
+    case 900: 
+      if      (aEta <= 0.51)  g = CMSsqrts900eta05();
+      else if (aEta <= 1.01)  g = CMSsqrts900eta10();
+      else if (aEta <= 1.51)  g = CMSsqrts900eta15();
+      else if (aEta <= 2.01)  g = CMSsqrts900eta20();
+      else if (aEta <= 2.41)  g = CMSsqrts900eta24();
+      break;
+    default: 
+      return g;
+    }
+    return g;
+  }
+
+  /* ===================================================================
+   * 
+   * ALICE data 
+   */
+  /**
+   * ALICE @f$ P(N_{ch})@f$ at @f$\sqrt{s}=900GeV@f$ for @f$|\eta|\lt0.5@f$
+   * http://hepdata.cedar.ac.uk/view/p7742/d8
+   *
+   * @return p7742_d8x1y1
+   */
+  TGraphAsymmErrors* ALICEsqrts900eta05()
+  {
+    double x[] = { 0.0,        1.0,
+                  2.0, 3.0,
+                  4.0, 5.0,
+                  6.0, 7.0,
+                  8.0, 9.0,
+                  10.0,        11.0,
+                  12.0,        13.0,
+                  14.0,        15.0,
+                  16.0,        17.0,
+                  18.0,        19.0,
+                  20.0,        21.5,
+                  23.5,        25.5 };
+    double xem[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.0,
+                    1.0,       1.0 };
+    double xep[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.0,
+                    1.0,       1.0 };
+    double y[] = { 0.179767,   0.155985,
+                  0.140806,    0.120933,
+                  0.097273,    0.075258,
+                  0.057395,    0.043794,
+                  0.033212,    0.024939,
+                  0.018642,    0.013927,
+                  0.010337,    0.007693,
+                  0.005665,    0.004134,
+                  0.002994,    0.002144,
+                  0.001527,    0.001093,
+                  8.03E-4,     4.82E-4,
+                  2.56E-4,     1.02E-4 };
+    double yem[] = { 0.028627623478032542,     0.006826681770816624,
+                    0.004736912179891031,      0.005286866936097409,
+                    0.004605827830043151,      0.003753897441326814,
+                    0.003172958871463669,      0.0023132784095305087,
+                    0.0017387768689512751,     0.0012885437516824954,
+                    9.891637882575362E-4,      7.64476945368531E-4,
+                    5.769662035162892E-4,      4.500277769204919E-4,
+                    3.7107950630558943E-4,     3.1458703088334716E-4,
+                    2.7004073766748605E-4,     2.210701246211256E-4,
+                    1.7819090885900997E-4,     1.3576450198781714E-4,
+                    1.2447088012864694E-4,     9.762171889492625E-5,
+                    6.095900261651268E-5,      3.4828149534536E-5 };
+    double yep[] = { 0.028627623478032542,     0.006826681770816624,
+                    0.004736912179891031,      0.005286866936097409,
+                    0.004605827830043151,      0.003753897441326814,
+                    0.003172958871463669,      0.0023132784095305087,
+                    0.0017387768689512751,     0.0012885437516824954,
+                    9.891637882575362E-4,      7.64476945368531E-4,
+                    5.769662035162892E-4,      4.500277769204919E-4,
+                    3.7107950630558943E-4,     3.1458703088334716E-4,
+                    2.7004073766748605E-4,     2.210701246211256E-4,
+                    1.7819090885900997E-4,     1.3576450198781714E-4,
+                    1.2447088012864694E-4,     9.762171889492625E-5,
+                    6.095900261651268E-5,      3.4828149534536E-5 };
+    int np = 24;
+    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
+                                                x,
+                                                y,
+                                                xem,
+                                                xep,
+                                                yem,
+                                                yep);
+    g->SetName("/HepData/7742/d8x1y1");
+    g->SetTitle("/HepData/7742/d8x1y1");
+    return g;
+  }
+  /**
+   * ALICE @f$ P(N_{ch})@f$ at @f$\sqrt{s}=900GeV@f$ for @f$|\eta|\lt1.0@f$
+   * http://hepdata.cedar.ac.uk/view/p7742/d9
+   *
+   * @return p7742_d9x1y1
+   */
+  TGraphAsymmErrors* ALICEsqrts900eta10()
+  {
+    double x[] = { 0.0,        1.0,
+                  2.0, 3.0,
+                  4.0, 5.0,
+                  6.0, 7.0,
+                  8.0, 9.0,
+                  10.0,        11.0,
+                  12.0,        13.0,
+                  14.0,        15.0,
+                  16.0,        17.0,
+                  18.0,        19.0,
+                  20.0,        21.0,
+                  22.0,        23.0,
+                  24.0,        25.0,
+                  26.0,        27.0,
+                  28.0,        29.0,
+                  30.0,        31.0,
+                  32.0,        33.0,
+                  34.0,        35.0,
+                  36.0,        37.0,
+                  38.0,        39.0,
+                  40.0,        41.5 };
+    double xem[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.0 };
+    double xep[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.0 };
+    double y[] = { 0.092396,   0.076087,
+                  0.073233,    0.078506,
+                  0.080885,    0.077084,
+                  0.070058,    0.062258,
+                  0.054437,    0.047195,
+                  0.040978,    0.035352,
+                  0.03048,     0.026336,
+                  0.022788,    0.019725,
+                  0.017071,    0.014674,
+                  0.012509,    0.010593,
+                  0.008935,    0.007527,
+                  0.006361,    0.005416,
+                  0.004665,    0.004041,
+                  0.003519,    0.00306,
+                  0.002637,    0.002234,
+                  0.001856,    0.001512,
+                  0.001213,    9.62E-4,
+                  7.61E-4,     6.02E-4,
+                  4.84E-4,     3.93E-4,
+                  3.27E-4,     2.76E-4,
+                  2.37E-4,     1.7E-4 };
+    double yem[] = { 0.021888152137629163,     0.0066076761421849355,
+                    0.004451948000594796,      0.0034909631908686747,
+                    0.002952765483407038,      0.0028889466938661224,
+                    0.0029957812002881653,     0.0032637882897026274,
+                    0.002922475833946279,      0.0026227742945209753,
+                    0.0023487871338203465,     0.002091141554271255,
+                    0.0017347841364273539,     0.0014762689456870655,
+                    0.0012594490859101847,     0.0010506721658062519,
+                    9.019456746389995E-4,      7.836938177630343E-4,
+                    6.432239112470866E-4,      5.679656679765071E-4,
+                    4.967293025381128E-4,      4.475187146924696E-4,
+                    4.040519768544636E-4,      3.3646693745448454E-4,
+                    2.9206163733020465E-4,     2.820017730440715E-4,
+                    2.6982401672201087E-4,     2.5470178640912594E-4,
+                    2.3519353732617736E-4,     2.1914607000811127E-4,
+                    1.9727392123643713E-4,     1.7402586014727813E-4,
+                    1.5060212481900778E-4,     1.3009611831257687E-4,
+                    1.1269871339105873E-4,     9.693812459502196E-5,
+                    8.48528137423857E-5,       7.427651041883969E-5,
+                    6.726812023536856E-5,      6.030754513325841E-5,
+                    5.544366510251645E-5,      4.601086828130936E-5 };
+    double yep[] = { 0.021888152137629163,     0.0066076761421849355,
+                    0.004451948000594796,      0.0034909631908686747,
+                    0.002952765483407038,      0.0028889466938661224,
+                    0.0029957812002881653,     0.0032637882897026274,
+                    0.002922475833946279,      0.0026227742945209753,
+                    0.0023487871338203465,     0.002091141554271255,
+                    0.0017347841364273539,     0.0014762689456870655,
+                    0.0012594490859101847,     0.0010506721658062519,
+                    9.019456746389995E-4,      7.836938177630343E-4,
+                    6.432239112470866E-4,      5.679656679765071E-4,
+                    4.967293025381128E-4,      4.475187146924696E-4,
+                    4.040519768544636E-4,      3.3646693745448454E-4,
+                    2.9206163733020465E-4,     2.820017730440715E-4,
+                    2.6982401672201087E-4,     2.5470178640912594E-4,
+                    2.3519353732617736E-4,     2.1914607000811127E-4,
+                    1.9727392123643713E-4,     1.7402586014727813E-4,
+                    1.5060212481900778E-4,     1.3009611831257687E-4,
+                    1.1269871339105873E-4,     9.693812459502196E-5,
+                    8.48528137423857E-5,       7.427651041883969E-5,
+                    6.726812023536856E-5,      6.030754513325841E-5,
+                    5.544366510251645E-5,      4.601086828130936E-5 };
+    int np = 42;
+    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
+                                                x,
+                                                y,
+                                                xem,
+                                                xep,
+                                                yem,
+                                                yep);
+    g->SetName("/HepData/7742/d9x1y1");
+    g->SetTitle("/HepData/7742/d9x1y1");
+    return g;
+  }
+
+  /**
+   * ALICE @f$ P(N_{ch})@f$ at @f$\sqrt{s}=900GeV@f$ for @f$|\eta|\lt1.3@f$
+   * http://hepdata.cedar.ac.uk/view/p7742/d10
+   *
+   * @return p7742_d10x1y1
+   */
+  TGraphAsymmErrors* ALICEsqrts900eta13()
+  {
+    double x[] = { 0.0,        1.0,
+                  2.0, 3.0,
+                  4.0, 5.0,
+                  6.0, 7.0,
+                  8.0, 9.0,
+                  10.0,        11.0,
+                  12.0,        13.0,
+                  14.0,        15.0,
+                  16.0,        17.0,
+                  18.0,        19.0,
+                  20.0,        21.0,
+                  22.0,        23.0,
+                  24.0,        25.0,
+                  26.0,        27.0,
+                  28.0,        29.0,
+                  30.0,        31.0,
+                  32.0,        33.0,
+                  34.0,        35.0,
+                  36.0,        37.0,
+                  38.0,        39.0,
+                  40.0,        41.5,
+                  43.5,        45.5,
+                  47.5,        49.5,
+                  51.5,        53.5 };
+    double xem[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.0,
+                    1.0,       1.0,
+                    1.0,       1.0,
+                    1.0,       1.0 };
+    double xep[] = { 0.5,      0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       0.5,
+                    0.5,       1.0,
+                    1.0,       1.0,
+                    1.0,       1.0,
+                    1.0,       1.0 };
+    double y[] = { 0.068695,   0.056382,
+                  0.049828,    0.056116,
+                  0.063345,    0.066198,
+                  0.065598,    0.061117,
+                  0.055042,    0.04927,
+                  0.044132,    0.039615,
+                  0.035906,    0.03262,
+                  0.029489,    0.026657,
+                  0.023839,    0.021302,
+                  0.018884,    0.016732,
+                  0.014707,    0.012944,
+                  0.01137,     0.009982,
+                  0.008753,    0.007661,
+                  0.006711,    0.005891,
+                  0.005179,    0.004572,
+                  0.004064,    0.003636,
+                  0.003278,    0.002968,
+                  0.002677,    0.002392,
+                  0.002111,    0.001829,
+                  0.001561,    0.001319,
+                  0.001115,    7.93E-4,
+                  5.24E-4,     3.51E-4,
+                  2.2E-4,      1.56E-4,
+                  1.27E-4,     8.5E-5 };
+    double yem[] = { 0.021322725646595934,     0.007630831737104416,
+                    0.004038016592338372,      0.0036354708085748677,
+                    0.0035003511252444373,     0.0030936373413831173,
+                    0.0026374775828431228,     0.002920906195001818,
+                    0.00276560481631053,       0.002864668392676542,
+                    0.0026843833556330957,     0.0022856071403458645,
+                    0.002325611317481922,      0.002006865217198205,
+                    0.0017203360136903488,     0.0014989629748596194,
+                    0.0012394615766533467,     0.0011637718848640398,
+                    9.151568171630477E-4,      9.078199160626517E-4,
+                    7.98723982361867E-4,       6.967302203866285E-4,
+                    6.308502199413107E-4,      5.846990678973245E-4,
+                    5.4516236113657E-4,        5.460009157501477E-4,
+                    4.930608481719068E-4,      4.3796118549478783E-4,
+                    3.4394330928221294E-4,     3.2323366161339073E-4,
+                    3.0869078379504626E-4,     2.9443165590676553E-4,
+                    2.756320010448714E-4,      2.630304164920856E-4,
+                    2.5628889948649746E-4,     2.4142286552851616E-4,
+                    2.245551157288562E-4,      2.0678491240900531E-4,
+                    1.8828170383762728E-4,     1.7255723688098393E-4,
+                    1.5911316727411343E-4,     1.3023056476879765E-4,
+                    1.0761505470890214E-4,     8.62670273047588E-5,
+                    6.129437168288782E-5,      6.037383539249432E-5,
+                    3.535533905932738E-5,      3.8078865529319545E-5 };
+    double yep[] = { 0.021322725646595934,     0.007630831737104416,
+                    0.004038016592338372,      0.0036354708085748677,
+                    0.0035003511252444373,     0.0030936373413831173,
+                    0.0026374775828431228,     0.002920906195001818,
+                    0.00276560481631053,       0.002864668392676542,
+                    0.0026843833556330957,     0.0022856071403458645,
+                    0.002325611317481922,      0.002006865217198205,
+                    0.0017203360136903488,     0.0014989629748596194,
+                    0.0012394615766533467,     0.0011637718848640398,
+                    9.151568171630477E-4,      9.078199160626517E-4,
+                    7.98723982361867E-4,       6.967302203866285E-4,
+                    6.308502199413107E-4,      5.846990678973245E-4,
+                    5.4516236113657E-4,        5.460009157501477E-4,
+                    4.930608481719068E-4,      4.3796118549478783E-4,
+                    3.4394330928221294E-4,     3.2323366161339073E-4,
+                    3.0869078379504626E-4,     2.9443165590676553E-4,
+                    2.756320010448714E-4,      2.630304164920856E-4,
+                    2.5628889948649746E-4,     2.4142286552851616E-4,
+                    2.245551157288562E-4,      2.0678491240900531E-4,
+                    1.8828170383762728E-4,     1.7255723688098393E-4,
+                    1.5911316727411343E-4,     1.3023056476879765E-4,
+                    1.0761505470890214E-4,     8.62670273047588E-5,
+                    6.129437168288782E-5,      6.037383539249432E-5,
+                    3.535533905932738E-5,      3.8078865529319545E-5 };
+    int np = 48;
+    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
+                                                x,
+                                                y,
+                                                xem,
+                                                xep,
+                                                yem,
+                                                yep);
+    g->SetName("/HepData/7742/d10x1y1");
+    g->SetTitle("/HepData/7742/d10x1y1");
+    return g;
+  }
+  /** 
+   * Get ALICE data 
+   * 
+   * @param eta Eta limit
+   * @param sNN Center of mass energy 
+   * 
+   * @return 
+   */
+  TGraphAsymmErrors* GetALICE(Double_t eta, UShort_t sNN)
+  {
+    Double_t aEta = TMath::Abs(eta);
+    TGraphAsymmErrors* g = 0;
+    switch (sNN) { 
+    case 900:
+      if      (aEta <= 0.51) g = ALICEsqrts900eta05();
+      else if (aEta <= 1.01) g = ALICEsqrts900eta10();
+      else if (aEta <= 1.31) g = ALICEsqrts900eta13();
+      break;
+    }
+    return g;
+  }
+
+  /* =================================================================
+   *
+   * Measurements from other sources, such as published ALICE, CMS, ...
+   */
+  static TGraphAsymmErrors* GetData(UShort_t type, Double_t eta, UShort_t sNN)
+  {
+    OtherPNch d;
+    TGraphAsymmErrors* g = 0;
+    switch (type) { 
+    case 0: g = d.GetCMS(eta, sNN);   break;
+    case 1: g = d.GetALICE(eta, sNN); break;
+    }
+    if (!g) {
+      // Warning("GetOther", "No other data found for type=%d, eta=%f, sNN=%d", 
+      //         type, eta, sNN);
+      return g;
+    }
+    return g;
+  }    
+};
+
+//
+// EOF
+//
index 58fe7c0f17226609ee9b2db768e4fb97e21571b1..fa3f8507c6a37de97c701b273974740ea10a8a9a 100644 (file)
@@ -36,9 +36,10 @@ public:
     // --- Force MB for pp ---------------------------------------------
     TCollection* c   = GetCollection(file, "ForwardMultSums");
     
-    fBody->Divide(1,2);
-    DrawInPad(fBody, 1, GetH1(c, "triggers"));
-    DrawInPad(fBody, 2, GetH1(c, "status"));
+    fBody->Divide(1,3);
+    DrawInPad(fBody, 1, GetH1(c, "triggers"),   "hist text30");
+    DrawInPad(fBody, 2, GetH1(c, "status"),     "hist text30");
+    DrawInPad(fBody, 3, GetH1(c, "diagnostics"),"colz text");
     PrintCanvas("Overview");
 
     DrawSumCollection(c, "symmetric");
@@ -108,8 +109,9 @@ protected:
       if (!bin) continue;
 
       fBody->Divide(2,2);
-      DrawInPad(fBody, 1, GetH1(bin, "rawDist"), "",     kLogy);
-      DrawInPad(fBody, 1, GetH1(bin, "truth"),   "same", kLogy|kLegend);
+      DrawInPad(fBody, 1, GetH1(bin, "rawDist"), "",          kLogy);
+      DrawInPad(fBody, 1, GetH1(bin, "truthAccepted"),"same", kLogy);
+      DrawInPad(fBody, 1, GetH1(bin, "truth"),        "same", kLogy|kLegend);
       DrawInPad(fBody, 2, GetH1(bin, "coverage"));
       DrawInPad(fBody, 3, GetH2(bin, "corr"),     "colz");
       DrawInPad(fBody, 4, GetH2(bin, "response"), "colz", kLogz);
@@ -148,13 +150,16 @@ protected:
       e = l->AddEntry("dummy", "MC truth", "p");
       e->SetMarkerStyle(24);
       e->SetMarkerColor(kBlue+1);
+      e = l->AddEntry("dummy", "MC truth selected", "p");
+      e->SetMarkerStyle(24);
+      e->SetMarkerColor(kOrange+1);
     }
     fBody->cd();
     l->Draw();
     
     PrintCanvas(Form("%s results", name.Data()));
 
-    return;
+    // return;
 
     TIter       nextO(c);
     TObject*    o = 0;
@@ -164,12 +169,14 @@ protected:
       TCollection* bin = GetEtaBin(o, etaMin, etaMax);
       if (!bin) continue;
 
-      fBody->Divide(2,2);
-      DrawInPad(fBody, 1, GetH1(bin, "rawDist"), "",     kLogy);
-      DrawInPad(fBody, 1, GetH1(bin, "truth"),   "same", kLogy);
+      fBody->Divide(2,3);
+      DrawInPad(fBody, 1, GetH1(bin, "rawDist"),      "",     kLogy);
+      DrawInPad(fBody, 1, GetH1(bin, "truthAccepted"),"same", kLogy);
+      DrawInPad(fBody, 1, GetH1(bin, "truth"),        "same", kLogy|kLegend);
       DrawInPad(fBody, 2, GetH1(bin, "coverage"));
       DrawInPad(fBody, 3, GetH2(bin, "corr"),     "colz");
       DrawInPad(fBody, 4, GetH2(bin, "response"), "colz", kLogz);
+      DrawInPad(fBody, 5, GetH1(bin, "triggerVertex"));
       
       PrintCanvas(Form("%+5.1f < #eta < %+5.1f", etaMin, etaMax));
     }
diff --git a/PWGLF/FORWARD/analysis2/scripts/SummaryUnfoldedDrawer.C b/PWGLF/FORWARD/analysis2/scripts/SummaryUnfoldedDrawer.C
new file mode 100644 (file)
index 0000000..f39eb95
--- /dev/null
@@ -0,0 +1,243 @@
+#include "SummaryDrawer.C"
+#include <TMultiGraph.h>
+#include <TKey.h>
+#include <TList.h>
+
+struct SummaryUnfoldedDrawer : public SummaryDrawer
+{
+
+  SummaryUnfoldedDrawer() : SummaryDrawer() {}
+  void Run(const char* fname)
+  {
+    TString filename(fname);
+    TFile* file = TFile::Open(filename,"READ");
+    if (!file) { 
+      Error("Run", "Failed to open \"%s\"", filename.Data());
+      return;
+    }
+    
+    TString pdfName(filename);
+    pdfName.ReplaceAll(".root", ".pdf");
+    CreateCanvas(pdfName, 0); // what & kLandscape);
+
+    // Title page 
+    fBody->cd();
+
+    TLatex* ltx = new TLatex(.5, .7, "Raw P(#it{N}_{ch}) #rightarrow "
+                            "corrected P(#it{N}_{ch})");
+    ltx->SetNDC();
+    ltx->SetTextSize(0.07);
+    ltx->SetTextAlign(22);
+    ltx->Draw();
+
+    Double_t y = .6;
+    
+    Double_t save = fParName->GetTextSize();
+    fParName->SetTextSize(0.03);
+    fParVal->SetTextSize(0.03);
+
+    TObject* method = file->Get("method");
+    TObject* sys    = file->Get("sys");
+    TObject* sNN    = file->Get("sNN");
+    TObject* trig   = file->Get("trigger");
+    Double_t regP;
+    Double_t minIpZ;
+    Double_t maxIpZ;
+    UShort_t maxN;
+    GetParameter(file, "regParam", regP);
+    GetParameter(file, "minIpZ", minIpZ);
+    GetParameter(file, "maxIpZ", maxIpZ);
+    GetParameter(file, "maxN", maxN);
+
+    DrawParameter(y, "System", sys->GetTitle());
+    DrawParameter(y, "#sqrt{s_{NN}}", sNN->GetTitle());
+    DrawParameter(y, "Trigger", trig->GetTitle());
+    DrawParameter(y, "Method", method->GetTitle());
+    DrawParameter(y, "Reg. param.", Form("%g", regP));
+    DrawParameter(y, "Max #it{N}_{ch}", Form("%d", maxN));
+    DrawParameter(y, "IP_{z} range", Form("%+5.2fcm - %+5.2fcm", 
+                                         minIpZ, maxIpZ));
+
+
+    PrintCanvas("Title page");
+    fParName->SetTextSize(save);
+    fParVal->SetTextSize(save);
+
+    
+    
+    // Loop over all keys in the file 
+    TIter next(file->GetListOfKeys());
+    TKey* k  = 0;
+    while ((k = static_cast<TKey*>(next()))) {
+      file->cd();
+
+      TString clName = k->GetClassName();
+      TClass* cl     = gROOT->GetClass(clName);
+      if (!cl) { 
+       // Warning("Run", "Ignoring object %s of unknown type %s", 
+       //         k->GetName(), clName.Data());
+       continue;
+      }
+      if (!cl->InheritsFrom(TDirectory::Class())) {
+       // Warning("Run", "Ignoring object %s of type %s", 
+       //         k->GetName(), clName.Data());
+       continue;
+      }
+      file->cd(k->GetName());
+      ProcessType(gDirectory);
+    }
+    file->cd();
+    CloseCanvas();
+  }
+  void ProcessType(TDirectory* d) 
+  { 
+    Printf(" ProcessType: %s", d->GetPath());
+    // d->ls();
+    MakeChapter(d->GetName()); 
+    static TRegexp regex("[pm][0-9]d[0-9]*_[pm][0-9]d[0-9]*");
+    
+    TIter next(d->GetListOfKeys());
+    TKey* k  = 0;
+    while ((k = static_cast<TKey*>(next()))) {
+      d->cd();
+      TString clName = k->GetClassName();
+      TClass* cl     = gROOT->GetClass(clName);
+      if (!cl) { 
+       // Warning("ProcessType", "Ignoring object %s of unknown type %s", 
+       // k->GetName(), clName.Data());
+       continue;
+      }
+      if (!cl->InheritsFrom(TDirectory::Class())) {
+       // Warning("ProcessType", "Ignoring object %s of type %s", 
+       //  k->GetName(), clName.Data());
+       continue;
+      }
+      TString n(k->GetName());
+      if (n.Index(regex) == kNPOS) { 
+       Warning("ProcessType", "Ignoring non-bin directory %s", n.Data());
+       continue;
+      }
+      d->cd(n);
+      ProcessBin(gDirectory);
+    }
+    d->cd();
+    DrawResults(d);
+  }
+  void DrawResults(TDirectory* d) 
+  { 
+    DrawInPad(fBody, 0, GetStack(d, "corrected"), "nostack", kLogy);
+    TObject* alice = d->Get("alice");
+    TObject* cms   = d->Get("cms");
+    if (cms)   cms->Draw();
+    if (alice) alice->Draw();
+
+    if (alice || cms) { 
+      TObject* dummy = 0;
+      TLegend* l = new TLegend(0.15, 0.12, .4, .3, "", "NDC");
+      l->SetFillColor(0);
+      l->SetFillStyle(0);
+      l->SetBorderSize(0);
+      TLegendEntry* e = l->AddEntry(dummy, "This work", "F");
+      e->SetFillColor(kRed+2);
+      e->SetFillStyle(1001);
+      if (alice) { 
+       e = l->AddEntry(dummy, "ALICE", "F");
+       e->SetFillColor(kPink+1);
+       e->SetFillStyle(1001);
+      }
+      if (cms) { 
+       e = l->AddEntry(dummy, "CMS", "F");
+       e->SetFillColor(kGreen+2);
+       e->SetFillStyle(1001);
+      }
+      l->Draw();
+    }
+    PrintCanvas(" Results");
+  }
+  void ProcessBin(TDirectory* d) 
+  {
+    Printf("  ProcessBin: %s", d->GetPath());
+
+    TString tmp(d->GetName());
+    tmp.ReplaceAll("p", "+");
+    tmp.ReplaceAll("m", "-");
+    tmp.ReplaceAll("d", ".");
+    tmp.ReplaceAll("_", " ");
+    TObjArray* tokens = tmp.Tokenize(" ");
+    if (!tokens || tokens->GetEntriesFast() < 2) { 
+      Error("Other2Stack", "Failed to decode eta range from %s", 
+           d->GetName());
+      if (tokens) tokens->Delete();
+      return;
+    }
+    Double_t eta1 = static_cast<TObjString*>(tokens->At(0))->String().Atof();
+    Double_t eta2 = static_cast<TObjString*>(tokens->At(1))->String().Atof();
+    tokens->Delete();
+
+    fBody->Divide(1,5,0,0);
+    TVirtualPad* p = fBody->cd(5);
+    p->SetRightMargin(0.15);
+    p->SetTopMargin(0.05);
+    p = fBody->cd(4);
+    p->SetBottomMargin(0.15);
+    THStack* all = GetStack(d,"all");
+    DrawInPad(fBody,1,all, "nostack", kLogy);
+    DrawInPad(fBody,2,GetH1(d,"ratioCorrTruth"), "", 0);
+    DrawInPad(fBody,3,GetH1(d,"ratioUnfAcc"), "", 0);
+    DrawInPad(fBody,4,GetH1(d,"triggerVertex"), "", 0);
+    DrawInPad(fBody,5,GetH2(d,"response"), "colz", kLogz);
+
+    PrintCanvas(Form(" %+5.2f<eta<%+5.2f", eta1, eta2));
+    
+    DrawSteps(all, eta1, eta2);
+  }
+  void DrawSteps(THStack* stack, Double_t e1, Double_t e2) 
+  {
+    fBody->Divide(2, 3, 0, 0);
+    TList* hists = stack->GetHists();
+    Int_t  nHist = hists->GetEntries(); 
+    
+    // Make a background stack 
+    THStack* bg = static_cast<THStack*>(stack->Clone("bg"));
+    bg->SetTitle();
+    for (Int_t i = 0; i < nHist; i++) { 
+      // Loop over histograms and set the saved color
+      TH1* h        = static_cast<TH1*>(bg->GetHists()->At(i));
+      h->SetMarkerColor(kGray+1);
+      h->SetFillColor(kGray);
+      h->SetLineColor(kGray);
+
+      TList* lf = h->GetListOfFunctions();
+      if (lf) { 
+       TObject* ll = lf->FindObject("legend");
+       if (ll) lf->Remove(ll);
+      }
+    }
+    const char* txt[] = { "MC 'truth'", 
+                         "Selected MC 'truth'", 
+                         "Measured", 
+                         "Unfolded", 
+                         "Corrected" };
+    // Now loop again, this time drawing the stack
+    for (Int_t i = 0; i < nHist; i++) { 
+      
+      DrawInPad(fBody, i+1, bg, "nostack", kLogy);
+      DrawInPad(fBody, i+1, hists->At(i), "same hist p e", kLogy);
+      gPad->SetGridx();
+      gPad->SetGridy();
+
+      TLatex* l = new TLatex(.95, .95, Form("Step %d", i));
+      l->SetNDC();
+      l->SetTextAlign(33);
+      l->SetTextSize(0.06);
+      l->Draw();
+      l->DrawLatex(.95, .88, txt[i]);
+      
+    }
+    PrintCanvas(Form(" %+5.2f<eta<%+5.2f - Steps", e1, e2));
+  }
+};
+//
+// EOF
+//
+
index 6123a19b8f4cdb6eba388cedcc7363001022af98..4ceb8853bf218d85f035fea8d08d6c842e3dc0ed 100644 (file)
-#ifndef UNFOLDMULTDISTS_H
-#define UNFOLDMULTDISTS_H
 #include <TFile.h>
 #include <TList.h>
 #include <TH1.h>
 #include <TH2.h>
-#include <TString.h>
-#include <TGraph.h>
-#include <TGraphAsymmErrors.h>
+#include <THStack.h>
 #include <TLegend.h>
 #include <TLegendEntry.h>
-#include <TLatex.h>
+#include <TClass.h>
+#include <TRegexp.h>
+#include <TMath.h>
 #include <TParameter.h>
-#include <TCanvas.h>
-#include <THStack.h>
 #include <TMultiGraph.h>
-#include <TRegexp.h>
-#include <TSystem.h>
+#include <TGraphAsymmErrors.h>
 #include "RooUnfold.h"
 #include "RooUnfoldResponse.h"
+#include <fstream>
 
-/** 
- * Structure to do unfolding of multiplcity distribution
+/**
+ * Class to do unfolding of raw histograms produced by AliForwardMultDist
  * 
- * It takes as input 2 files: the output file (forward_multdists.root)
- * of the analysis train MakeMultDistsTrain.C run on real and MC data.
- *
- * It generates the output file forward_unfolded.root with histograms
- * and stacks of the input data and the unfolding results.
  */
-struct Unfolder 
+struct Unfolder
 {
-  enum { 
-    kMeasuredColor = kBlue+1,
-    kUnfoldedColor = kRed+1,
-    kTruthColor    = kGreen+2,
-    kCMSColor      = kYellow+1,
-    kALICEColor    = kMagenta+1,
-    kSysColor      = kBlue-10
+  /**
+   * Colours used 
+   * 
+   */
+  enum {
+    kColorMeasured = kOrange-2, 
+    kColorTruth    = kBlue-3, 
+    kColorAccepted = kMagenta-3,
+    kColorTrgVtx   = kBlack, 
+    kColorUnfolded = kOrange+2,
+    kColorCorrected= kRed+2, 
+    kColorError    = kBlue-10,
+    kColorALICE    = kPink+1, 
+    kColorCMS      = kGreen+2
   };
+
   /** 
-   * Constructor
+   * Constructor 
    */
   Unfolder() {}
-  virtual ~Unfolder() {}
   /** 
-   * Run this code 
+   * Get a top collection from a file
    * 
-   * @param realFile Output file of MakeMultDistsTrain on real data
-   * @param mcFile   Output file of MakeMultDistsTrain on MC data
+   * @param fileName Name of file 
+   * @param results  Wheter it's the results collection or not 
+   * 
+   * @return Collection or null
    */
-  void Run(const TString& method="Bayes", 
-          Double_t       regParam=-1e30, 
-          const TString& realFile="forward_multdists.root", 
-          const TString& mcFile="forward_mcmultdists.root")
+  static TCollection* GetTop(const TString& fileName, Bool_t results=false)
   {
-    if (!gROOT->GetClass("RooUnfold")) gSystem->Load("libRooUnfold.so");
-
-    // --- Open files ------------------------------------------------
-    TFile* realF = TFile::Open(realFile,"READ");
-    if (!realF) { 
-      Error("Run", "Couldn't open %s", realFile.Data());
-      return;
-    }
-    TFile* mcF = TFile::Open(mcFile,"READ");
-    if (!mcF) { 
-      Error("Run", "Couldn't open %s", mcFile.Data());
-      return;
-    }
-    TFile* outF = TFile::Open("forward_unfolded.root", "RECREATE");
-    
-    // --- Get top-level containers ----------------------------------
-    TCollection* realTop = GetCollection(realF, "ForwardMultResults");
-    if (!realTop) {
-      Error("Run", "Failed to get real collection");
-      return;
-    }
-    TCollection* mcTop   = GetCollection(mcF,   "ForwardMultResults");
-    if (!mcTop) { 
-      Error("Run", "Failed to get MC collection");
-      return;
-    }
-
-    // --- Decode the method -----------------------------------------
-    struct Method { 
-      UInt_t  id;
-      TString name;
-    };
-    const Method methods[] = { {RooUnfold::kNone,    "None"},
-                              {RooUnfold::kBayes,   "Bayes"},
-                              {RooUnfold::kSVD,     "SVD"},
-                              {RooUnfold::kBinByBin,"BinByBin"},
-                              {RooUnfold::kTUnfold, "TUnfold"},
-                              {RooUnfold::kInvert,  "Invert"},
-                              {RooUnfold::kDagostini,"Dagostini"}, 
-                              {0xDeadBeef,           "unknown"} };
-    const Method* pMethod = methods;
-    while (pMethod->id != 0xDeadBeef) {
-      if (method.BeginsWith(pMethod->name, TString::kIgnoreCase)) break;
-      pMethod++;
-    }
-    if (pMethod->id == 0xDeadBeef) {
-      Error("Run", "Unknown unfolding method: %s", method.Data());
-      return;
-    }
-    TNamed* methText = new TNamed("method", pMethod->name.Data());
-    methText->SetUniqueID(pMethod->id);
-    outF->cd();
-    methText->Write();
-    
-    // --- Loop over the kinds of bins we have -----------------------
-    const char*  types[] = { "symmetric", "negative", "positive", "other", 0 };
-    const char** ptype   = types;
-    Double_t     regParm = regParam;
-    while ((*ptype)) { 
-      TCollection* realType = GetCollection(realTop, *ptype);
-      TCollection* mcType   = GetCollection(mcTop, *ptype);
-      TDirectory*  outType  = outF->mkdir(*ptype);
-      if(realType && mcType) {
-       // Restore default 
-       regParm = regParam; 
-       // regParm is possibly modified here.
-       ScanType(pMethod->id, regParm, realType, mcType, outType);
-      }
-      // Info("run", "%s gave regularisation parameter %f", *ptype, regParm);
-      ptype++;
-      outF->cd();
+    TFile* file = TFile::Open(fileName, "READ");
+    if (!file) { 
+      Error("GetTop", "Failed to open %s", fileName.Data());
+      return 0;
     }
-    // --- Write used regularisation parameter -----------------------
-    // We get the value that's possibly modified
-    // Info("Run", "Regularisation parameter was %f, now %f", regParam,regParm);
-    TParameter<double>* regP = new TParameter<double>("regParam", regParm);
-    regP->Write();
-
-    // --- Close all files -------------------------------------------
-    // outF->ls();
-    outF->Write();
-    outF->Close();
-    realF->Close();
-    mcF->Close();
+    TCollection* ret = 0;
+    TString cName(Form("ForwardMult%s", results ? "Results" : "Sums"));
+    file->GetObject(cName, ret);
+    if (!ret) 
+      Error("GetTop", "Failed to get collection %s from %s", 
+           cName.Data(), fileName.Data());
+    file->Close();
+    return ret;
   }
   /** 
    * Get an object from a collection 
@@ -148,41 +74,23 @@ struct Unfolder
    * 
    * @return Pointer to object or null
    */
-  TObject* GetObject(TCollection* c, const TString& name, TClass* cl) const
+  static TObject* GetObject(TCollection* c, const TString& name, 
+                           TClass* cl, Bool_t verbose=true)
   {
     TObject* o = c->FindObject(name);
     if (!o) { 
-      Warning("GetObject", "%s not found in %s", name.Data(), c->GetName());
+      if (verbose)
+       Warning("GetObject", "%s not found in %s", name.Data(), c->GetName());
       return 0;
     }
     if (cl && !o->IsA()->InheritsFrom(cl)) {
-      Warning("GetCollection", "%s is not a %s but a %s", 
-             name.Data(), cl->GetName(), o->ClassName());
+      if (verbose) 
+       Warning("GetCollection", "%s is not a %s but a %s", 
+               name.Data(), cl->GetName(), o->ClassName());
       return 0;
     }
     return o;
   }
-  /** 
-   * Get a collection
-   * 
-   * @param d 
-   * @param name 
-   * 
-   * @return 
-   */   
-  TCollection* GetCollection(TDirectory* d, const TString& name) const
-  {
-    TObject* o = d->Get(name);
-    if (!o) { 
-      Warning("GetCollection", "%s not found in %s", name.Data(), d->GetName());
-      return 0;
-    }
-    if (!o->IsA()->InheritsFrom(TCollection::Class())) { 
-      Warning("GetCollection", "%s is not a collection", name.Data());
-      return 0;
-    }
-    return static_cast<TCollection*>(o);
-  }
   /** 
    * Get a collection 
    * 
@@ -191,9 +99,13 @@ struct Unfolder
    * 
    * @return 
    */
-  TCollection* GetCollection(TCollection* c, const TString& name) const
+  static TCollection* GetCollection(TCollection*   c, 
+                                   const TString& name, 
+                                   Bool_t         verbose=-true)
   {
-    return static_cast<TCollection*>(GetObject(c, name, TCollection::Class()));
+    return static_cast<TCollection*>(GetObject(c, name, 
+                                              TCollection::Class(),
+                                              verbose));
   }
   /** 
    * Get a 1D histogram from a collection
@@ -203,9 +115,9 @@ struct Unfolder
    * 
    * @return Pointer to object or null
    */
-  TH1* GetH1(TCollection* c, const TString& name) 
+  static TH1* GetH1(TCollection* c, const TString& name, Bool_t verbose=true) 
   {
-    return static_cast<TH1*>(GetObject(c, name, TH1::Class()));
+    return static_cast<TH1*>(GetObject(c, name, TH1::Class(), verbose));
   }
   /** 
    * Get a 2D histogram from a collection
@@ -215,98 +127,387 @@ struct Unfolder
    * 
    * @return Pointer to object or null
    */
-  TH2* GetH2(TCollection* c, const TString& name) 
+  static TH2* GetH2(TCollection* c, const TString& name, Bool_t verbose=true) 
   {
-    return static_cast<TH2*>(GetObject(c, name, TH2::Class()));
+    return static_cast<TH2*>(GetObject(c, name, TH2::Class(), verbose));
   }
-  TH1* Ratio(const TH1* num, const TGraph* denom,
-            Double_t etaMin, Double_t etaMax) 
+  /** 
+   * Get an unsigned short parameter from the collection 
+   * 
+   * @param c    Collection
+   * @param name Parameter name 
+   * 
+   * @return Value 
+   */
+  static void GetParameter(TCollection* c, const TString& name, UShort_t& v)
   {
-    TH1* ret = static_cast<TH1*>(num->Clone(Form("%s_%s",
-                                                num->GetName(), 
-                                                denom->GetName())));
-    ret->SetTitle(Form("%+5.1f<#eta<%+5.1f to %s", etaMin, etaMax,  
-                      denom->GetTitle()));
-    ret->SetDirectory(0);
-    ret->SetMarkerColor(denom->GetMarkerColor());
-    for (Int_t i = 1; i <= ret->GetNbinsX(); i++) { 
-      Double_t x    = ret->GetXaxis()->GetBinCenter(i);
-      Double_t numY = ret->GetBinContent(i);
-      Double_t numE = ret->GetBinError(i);
-      Double_t denY = denom->Eval(x);
-      
-      if (denY <= 0) {
-       ret->SetBinContent(i,0);
-       ret->SetBinError(i,0);
-       continue;
+    TObject* o = GetObject(c, name, TParameter<int>::Class(), true);
+    v = (!o ? 0 : o->GetUniqueID());
+  }
+  /** 
+   * Get an unsigned short parameter from the collection 
+   * 
+   * @param c    Collection
+   * @param name Parameter name 
+   * 
+   * @return Value 
+   */
+  static void GetParameter(TCollection* c, const TString& name, ULong_t& v)
+  {
+    TObject* o = GetObject(c, name, TParameter<long>::Class(), true);
+    v = (!o ? 0 : o->GetUniqueID());
+  }
+  /** 
+   * Get an unsigned short parameter from the collection 
+   * 
+   * @param c    Collection
+   * @param name Parameter name 
+   * 
+   * @return Value 
+   */
+  static void GetParameter(TCollection* c, const TString& name, Double_t& v)
+  {
+    TObject* o = GetObject(c, name, TParameter<double>::Class(), true);
+    v = (!o ? 0 : static_cast<TParameter<double>*>(o)->GetVal());
+  }
+  /** 
+   * Get the method identifier 
+   * 
+   * @param method 
+   * 
+   * @return 
+   */    
+  static UInt_t MethodId(TString& method) 
+  {
+    struct Method { 
+      UInt_t  id;
+      TString name;
+    };
+    const Method methods[] = { {RooUnfold::kNone,    "None"},
+                              {RooUnfold::kBayes,   "Bayes"},
+                              {RooUnfold::kSVD,     "SVD"},
+                              {RooUnfold::kBinByBin,"BinByBin"},
+                              {RooUnfold::kTUnfold, "TUnfold"},
+                              {RooUnfold::kInvert,  "Invert"},
+                              {RooUnfold::kDagostini,"Dagostini"}, 
+                              {0xDeadBeef,           "unknown"} };
+    const Method* pMethod = methods;
+    while (pMethod->id != 0xDeadBeef) {
+      if (method.BeginsWith(pMethod->name, TString::kIgnoreCase)) {
+       method = pMethod->name;
+       break;
       }
-      
-      ret->SetBinContent(i,(numY-denY)/denY);
-      ret->SetBinError(i,numE/denY);
+      pMethod++;
     }
-    return ret;
+    if (pMethod->id == 0xDeadBeef) 
+      Error("MethodId", "Unknown unfolding method: %s", method.Data());
+
+    return pMethod->id;
   }
+  
   /** 
-   * Scan type (symmetric, negative, ...) list for bins 
+   * Run the unfolding and correction task. 
+   *
+   * The @a measuredFile is assumed to have the structure 
+   *
+   * @verbatim
+   *     /- ForwardMultSums         (TCollection)
+   *          |- [type]             (TCollection)
+   *          |    |- [bin]         (TCollection) 
+   *          |    |    `- rawDist  (TH1)
+   *          |    |- [bin]
+   *          |    ...
+   *          |- [type]
+   *          ...
+   * @endverbatim 
    * 
-   * @param real Real list
-   * @param mc   MC list
-   * @param dir  Output directory 
-   */
-  void ScanType(UShort_t     method,
-               Double_t&    regParam,
-               TCollection* real, 
-               TCollection* mc, 
-               TDirectory*  dir) 
+   * and @a corrFile is assumed to have the structure 
+   *
+   * @verbatim
+   *     /- ForwardMultResults            (TCollection)
+   *          |- [type]                   (TCollection)
+   *          |    |- [bin]               (TCollection) 
+   *          |    |    |- truth          (TH1)
+   *          |    |    |- truthAccepted  (TH1)
+   *          |    |    |- triggerVertex  (TH1)
+   *          |    |    `- response       (TH2)
+   *          |    |- [bin]
+   *          |    ...
+   *          |- [type]
+   *          ...
+   * @endverbatim 
+   *
+   * where @c [type] is one of <i>symmetric</i>, <i>positive</i>,
+   * <i>negative</i>, or <i>other</i>, and [bin] is the @f$ \eta@f$
+   * bin named like
+   *
+   * @verbatim
+   *   [bin]          := [eta_spec] _ [eta_spec]
+   *   [eta_spec]     := [sign_char] [integer_part] d [decimal_part]
+   *   [sign_part]    := p          positive eta 
+   *                  |  m          negative eta 
+   *   [integer_part] := [number]
+   *   [decimal_part] := [number] [number]
+   *   [number]       := 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 
+   * @endverbatim
+   *
+   * That is, the bin @f$ -3\lt\eta\gt3@f$ is labeled
+   * <b>m3d00_p3d00</b>, @f$ 0\lt\eta\gt2.5@f$ is <b>p0d00_p2d50</b> 
+   *
+   * @a measuredFile and @a corrFile can point to the same file.  If
+   * @a corrFile is not specified, it is assumed that @a measuredFile
+   * has the expected @a corrFile <emph>in addition</emph> to the
+   * expected content of that file.
+   * 
+   * @param measuredFile Name of file containing measured data
+   * @param corrFile     Name of file containing correction data
+   * @param method       Unfolding method to use 
+   * @param regParam     Regularization parameter 
+   */  
+  void Run(const TString& measuredFile, const TString& corrFile,
+          const TString& method="Bayes", Double_t regParam=4)
   {
-    // --- Create container stack ------------------------------------
-    TString tit(real->GetName());
-    tit[0] = toupper(tit[0]);
-    tit.Append(" bins");
-    
-    THStack* stack = new THStack("all", tit);
-    dir->Add(stack);
+    // Get the input collections 
+    if (measuredFile.IsNull()) {
+      Error("Run", "No measurements given");
+      return;
+    }
+    TCollection* mTop = GetTop(measuredFile, false);
+    TCollection* cTop = GetTop(corrFile.IsNull() ? measuredFile : corrFile, 
+                              true);
+    if (!mTop || !cTop) return;
+
+    // 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); 
+    if (sys == 0 || sNN == 0) 
+      Warning("Run", "System (%d) and/or collision energy (%d) unknown", 
+             sys, sNN);
     
-    THStack* ratios = 0;
+    // Open the output file 
+    TFile* out = TFile::Open("forward_unfolded.root", "RECREATE");
+    if (!out) { 
+      Error("Run", "Failed to open output file");
+      return;
+    }
 
-    TMultiGraph* mg = 0;
+    // Decode method option and store in file 
+    TString meth(method);
+    UInt_t  mId = MethodId(meth);
+    if (mId == 0xDeadBeef) return;
+
+    // Store information 
+    SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,maxN);
+
+    // Load other data 
+    TString savPath(gROOT->GetMacroPath());
+    gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2/scripts",
+                             gROOT->GetMacroPath()));
+    // Always recompile 
+    if (!gROOT->GetClass("OtherPNch"))
+      gROOT->LoadMacro("OtherPNchData.C++");
+    gROOT->SetMacroPath(savPath);
+
+    // Loop over the input 
+    const char*  inputs[] = { "symmetric", "positive", "negative", 0 };
+    const char** pinput   = inputs;
+    while (*pinput) { 
+      TCollection* mInput = GetCollection(mTop, *pinput, false);
+      TCollection* cInput = GetCollection(cTop, *pinput, false);
+      if (mInput && cInput)
+       ProcessType(mInput, cInput, mId, regParam, out,
+                   sys, sNN);
+      pinput++;
+    }      
+
+    out->Write();
+    // out->ls();
+    out->Close();
+
+    SaveSummarize();
+  }
+  static void AppendAnd(TString& trg, const TString& what)
+  {
+    if (!trg.IsNull()) trg.Append(" & ");
+    trg.Append(what);
+  }
+  /** 
+   * Store some information on the output
+   * 
+   * @param dir      Where to store
+   * @param method   Method used
+   * @param mId      Method identifier 
+   * @param regParam Regularization parameter 
+   * @param sys      Collision system
+   * @param sNN      Center of mass energy 
+   * @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,
+                      Int_t          mId,
+                      Double_t       regParam, 
+                      UShort_t       sys, 
+                      UShort_t       sNN,
+                      UInt_t         trigger,
+                      Double_t       minIpZ, 
+                      Double_t       maxIpZ, 
+                      UShort_t       maxN) const
+  {
+    dir->cd();
 
-    // --- Create list of entries for legend -------------------------
-    TList* l = new TList;
-    l->SetName("legend");
+    TNamed* outMeth = new TNamed("method", method.Data());
+    outMeth->SetUniqueID(mId);
+    outMeth->Write();
 
-    // --- Get the bins container ------------------------------------
-    TList* le = static_cast<TList*>(real->FindObject("lowEdges"));
-    TList* he = static_cast<TList*>(real->FindObject("highEdges"));
-    if (!le || !he) 
-      Warning("ScanType", "didn't get the bin low/high edge containers");
-    else {
-      // Copy to output 
-      dir->cd();
-      le->Clone()->Write("lowEdges", TObject::kSingleKey);
-      he->Clone()->Write("highEdges", TObject::kSingleKey);
-    }
+    TParameter<double>* pR = new TParameter<double>("regParam", regParam);
+    pR->Write();
+    
+    TString tS = (sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "?");
+    TNamed* pS = new TNamed("sys", tS.Data()); pS->SetUniqueID(sys);
+    pS->Write();
+    
+    TString tE;
+    if      (sNN <  1000) tE = Form("%dGeV", sNN);
+    else if (sNN <  3000) tE = Form("%4.2fTeV", float(sNN)/1000);
+    else                  tE = Form("%dTeV", sNN/1000);
+    TNamed* pE = new TNamed("sNN", tE.Data()); pE->SetUniqueID(sNN);
+    pE->Write();
+    
+    TString tT;
+    /** 
+     * Bits of the trigger pattern
+     */
+    enum { 
+      /** In-elastic collision */
+      kInel        = 0x0001, 
+      /** In-elastic collision with at least one SPD tracklet */
+      kInelGt0     = 0x0002, 
+      /** Non-single diffractive collision */
+      kNSD         = 0x0004, 
+      /** Empty bunch crossing */
+      kEmpty       = 0x0008, 
+      /** A-side trigger */
+      kA           = 0x0010, 
+      /** B(arrel) trigger */
+      kB           = 0x0020, 
+      /** C-side trigger */
+      kC           = 0x0080,  
+      /** Empty trigger */
+      kE           = 0x0100,
+      /** pileup from SPD */
+      kPileUp      = 0x0200,    
+      /** true NSD from MC */
+      kMCNSD       = 0x0400,    
+      /** Offline MB triggered */
+      kOffline     = 0x0800,
+      /** At least one SPD cluster */ 
+      kNClusterGt0 = 0x1000,
+      /** V0-AND trigger */
+      kV0AND       = 0x2000, 
+      /** Satellite event */
+      kSatellite   = 0x4000
+    };
+    if ((trigger & kInel)        != 0x0) AppendAnd(tT, "INEL");
+    if ((trigger & kInelGt0)     != 0x0) AppendAnd(tT, "INEL>0");
+    if ((trigger & kNSD)         != 0x0) AppendAnd(tT, "NSD");
+    if ((trigger & kV0AND)       != 0x0) AppendAnd(tT, "V0AND");
+    if ((trigger & kA)           != 0x0) AppendAnd(tT, "A");
+    if ((trigger & kB)           != 0x0) AppendAnd(tT, "B");
+    if ((trigger & kC)           != 0x0) AppendAnd(tT, "C");
+    if ((trigger & kE)           != 0x0) AppendAnd(tT, "E");
+    if ((trigger & kMCNSD)       != 0x0) AppendAnd(tT, "MCNSD");
+    if ((trigger & kNClusterGt0) != 0x0) AppendAnd(tT, "NCluster>0");
+    if ((trigger & kSatellite)   != 0x0) AppendAnd(tT, "Satellite");
+    TNamed* pT = new TNamed("trigger", tT.Data()); pT->SetUniqueID(trigger);
+    pT->Write();
+    
+    TParameter<double>* pY = new TParameter<double>("minIpZ", minIpZ);
+    pY->Write();
 
-    // --- Setup for markers -----------------------------------------
-    const Int_t   nMarkers   = 7;
-    const Int_t   cMarkers[] = { 20,  21,  22,  23,  29,  33,  34  };
-    const Int_t   oMarkers[] = { 24,  25,  26,  32,  30,  27,  28  };
-    const Float_t sMarkers[] = { 1.1, 1.0, 1.2, 1.2, 1.2, 1.2, 1.0 };
-    Int_t iMarker            = 0;
+    TParameter<double>* pZ = new TParameter<double>("maxIpZ", maxIpZ);
+    pZ->Write();
 
-    // --- Containers that allow us to stack the objects properly ----
-    TList* mcTruth      = new TList;
-    TList* mcMeasured   = new TList;
-    TList* realMeasured = new TList;
-    TList* mcUnfolded   = new TList;
-    TList* realUnfolded = new TList;
-    // --- Loop over the contained objects ---------------------------
+    TParameter<int>* pN = new TParameter<int>("maxN", maxN);
+    pN->Write();
+  }
+  /** 
+   * Save a script to do a summary of this step 
+   * 
+   */                 
+  void SaveSummarize()
+  {
+    std::ofstream f("SummarizeUnfold.C");
+    f << "void SummarizeUnfold(const char* file=\"forward_unfolded.root\")\n"
+      << "{\n"
+      << "  const char* fwd=\"$ALICE_ROOT/PWGLF/FORWARD/analysis2\";\n"
+      << "  gROOT->LoadMacro(Form(\"%s/DrawUnfoldedSummary.C\",fwd));\n"
+      << "  DrawUnfoldedSummary(file);\n"
+      << "}\n"
+      << "// EOF" << std::endl;
+    f.close();
+  }
+  /** 
+   * Process a single type - i.e., one of <i>symmetric</i>,
+   * <i>positive</i>, <i>negative</i>, or <i>other</i> - by looping
+   * over all contained objects and call ProcessBin for each found
+   * bin.
+   * 
+   * @param measured     Input collection of measured data
+   * @param corrections  Input collection of correction data
+   * @param method       Unfolding method to use 
+   * @param regParam     Regularisation parameter
+   * @param out          Output directory. 
+   */
+  void ProcessType(TCollection* measured, 
+                  TCollection* corrections,
+                  UInt_t       method,
+                  Double_t     regParam,
+                  TDirectory*  out,
+                  UShort_t     sys, 
+                  UShort_t     sNN)
+  {
+    Printf("  Processing %s ...", measured->GetName());
+    TDirectory* dir = out->mkdir(measured->GetName());
+    
+    // Make some summary stacks 
+    THStack*  allMeasured  = new THStack("measured",      
+                                        "Measured P(#it{N}_{ch})");
+    THStack*  allTruth     = new THStack("truth",        
+                                        "MC 'truth' P(#it{N}_{ch})");
+    THStack*  allTruthA    = new THStack("truthAccepted",
+                                        "MC 'truth' accepted P(#it{N}_{ch})");
+    THStack*  allUnfolded  = new THStack("unfolded",
+                                        "Unfolded P(#it{N}_{ch})");
+    THStack*  allCorrected = new THStack("corrected",
+                                        "Corrected P(#it{N}_{ch})");
+    TMultiGraph* allALICE  = (sys != 1 ? 0 : 
+                             new TMultiGraph("alice", "ALICE Published"));
+    TMultiGraph* allCMS    = (sys != 1 ? 0 : 
+                             new TMultiGraph("cms", "CMS Published"));
+
+    // Loop over the list of objects. 
     static TRegexp regex("[pm][0-9]d[0-9]*_[pm][0-9]d[0-9]*");
-    TIter          next(real);
+    TIter          next(measured);
     TObject*       o = 0;
-    Int_t          f = 1;
+    Int_t          i = 0;
     Double_t       r = regParam;
     while ((o = next())) {
+      // Go back to where we where 
+      dir->cd();
+      
       // if not a collection, don't bother 
       if (!o->IsA()->InheritsFrom(TCollection::Class())) continue;
     
@@ -317,2085 +518,419 @@ struct Unfolder
        //         n.Data(), real->GetName());
        continue;
       }
-      // Cast object and find corresponding MC object 
-      TCollection* realBin = static_cast<TCollection*>(o);
-      TCollection* mcBin   = GetCollection(mc, n.Data());
-      if (!mcBin) { 
-       Warning("ScanType", "No corresponding MC bin for %s found in %s", 
-               n.Data(), mc->GetName());
-       continue;
-      }
-      TDirectory* outBin = dir->mkdir(realBin->GetName());
-      
-      // Restore default
-      r = regParam;
-      // Now do the unfolding 
-      THStack* bin = UnfoldEtaBin(method, r, realBin, mcBin, outBin);
-      if (!bin) { dir->cd(); continue; }
-
-      // Loop over histograms and set properties 
-      TIter   nextH(bin->GetHists());
-      TH1*    hist    = 0;
-      Int_t   cMarker = cMarkers[iMarker % nMarkers];
-      Int_t   oMarker = oMarkers[iMarker % nMarkers];
-      Float_t sMarker = sMarkers[iMarker % nMarkers];
-      TH1*    realH   = 0;
-      while ((hist = static_cast<TH1*>(nextH()))) {
-       TH1* out = static_cast<TH1*>(hist->Clone());
-       if (out->GetMarkerColor() == kUnfoldedColor) {
-          if (out->GetMarkerStyle() != 24)  {
-           realUnfolded->Add(out);
-           realH = out;
-         }
-         else 
-           mcUnfolded->Add(out);
-       }
-       else if (out->GetMarkerColor() == kMeasuredColor) {
-         if (out->GetMarkerStyle() != 24)  
-           realMeasured->Add(out);
-         else 
-           mcMeasured->Add(out);
-       }
-       else if (out->GetMarkerColor() == kTruthColor) 
-         mcTruth->Add(out);
-       else 
-         Warning("", "Unknown color for %s", out->GetName());
-
-       out->SetDirectory(0);
-       out->Scale(f);
-       out->SetMarkerStyle(out->GetMarkerStyle() == 24 ? oMarker : cMarker);
-       out->SetMarkerSize(out->GetMarkerSize() * sMarker);
-       out->SetOption(nextH.GetOption());
-       // stack->Add(out, nextH.GetOption());
-      }
-      TString nn(bin->GetTitle());
-      nn.Append(Form(" (#times%d)", f));
-      TObjString* lee = new TObjString(nn);
-      lee->SetUniqueID(cMarker);
-      l->Add(lee);
-
-      // Now try to get external data and make a multigraph 
-      nn = o->GetName();
-      nn.ReplaceAll("p", "+");
-      nn.ReplaceAll("m", "-");
-      nn.ReplaceAll("d", ".");
-      nn.ReplaceAll("_", " ");
-      TObjArray*  tokens = nn.Tokenize(" ");
-      TObjString* sMin   = static_cast<TObjString*>(tokens->At(0));
-      TObjString* sMax   = static_cast<TObjString*>(tokens->At(1));
-      Double_t    etaMin = sMin->String().Atof();
-      Double_t    etaMax = sMax->String().Atof();
-      if (TMath::Abs(etaMax + etaMin) < 1e-6) { 
-       // Symmetric bin
-       Double_t aEta = TMath::Abs(etaMin);
-       TGraphAsymmErrors* g1 = GetOther(0, aEta, 900, f, cMarker); 
-       TGraphAsymmErrors* g2 = GetOther(1, aEta, 900, f, cMarker); 
-       if (g1 || g2) {
-         if (!mg) mg = new TMultiGraph("other", "Other results");
-         if (!ratios) ratios = new THStack("ratios",tit);
-       }
-       if (g1) {
-         g1->SetTitle("CMS");
-         mg->Add(g1);
-         if (realH) ratios->Add(Ratio(realH, g1, etaMin, etaMax));
-       }
-       if (g2) {
-         g2->SetTitle("ALICE");
-         mg->Add(g2);
-         if (realH) ratios->Add(Ratio(realH, g2, etaMin, etaMax));
-       }
-      }
-
-      // Increment scaling and marker 
-      f *= 10;
-      iMarker++;
-
-      dir->cd();
+      TCollection* mBin = static_cast<TCollection*>(o);
+      TCollection* cBin = GetCollection(corrections, n.Data());
+      if (!cBin) continue;
+
+      THStack* binS = ProcessBin(mBin, cBin, method, r, dir);
+      if (!binS) continue;
+
+      Bin2Stack(binS, i, allMeasured, allTruth, allTruthA, 
+               allUnfolded, allCorrected);
+      Other2Stack(o->GetName(), i, sNN, allALICE, allCMS);
+      i++;
+    }
+    dir->Add(allMeasured);
+    dir->Add(allTruth);
+    dir->Add(allTruthA);
+    dir->Add(allUnfolded);
+    dir->Add(allCorrected);
+    if (allALICE && allALICE->GetListOfGraphs()) {
+      if (allALICE->GetListOfGraphs()->GetEntries() > 0)
+       dir->Add(allALICE);
+      else 
+       delete allALICE;
+    }
+    if (allCMS && allCMS->GetListOfGraphs()) {
+      if (allCMS->GetListOfGraphs()->GetEntries() > 0) 
+       dir->Add(allCMS);
+      else 
+       delete allCMS;
     }
-    regParam = r;
-    // Info("ScanType", "Regularisation parameter was %f, now %f", regParam, r);
-    TH1*  tmp = 0;
-    TIter nextMCT(mcTruth);
-    while ((tmp = static_cast<TH1*>(nextMCT()))) 
-      stack->Add(tmp, tmp->GetOption());
-    TIter nextMCM(mcMeasured);
-    while ((tmp = static_cast<TH1*>(nextMCM()))) 
-      stack->Add(tmp, tmp->GetOption());
-    TIter nextRM(realMeasured);
-    while ((tmp = static_cast<TH1*>(nextRM()))) 
-      stack->Add(tmp, tmp->GetOption());
-    TIter nextMCU(mcUnfolded);
-    while ((tmp = static_cast<TH1*>(nextMCU()))) 
-      stack->Add(tmp, tmp->GetOption());
-    TIter nextRU(realUnfolded);
-    while ((tmp = static_cast<TH1*>(nextRU()))) 
-      stack->Add(tmp, tmp->GetOption());
-    
-    dir->cd();
-    if (mg) mg->Write();
-    if (ratios) dir->Add(ratios);
-    l->Write(l->GetName(), TObject::kSingleKey);
   }
   /** 
-   * Do unfolding in an @f$\eta@f$ bin 
-   * 
-   * @param real Real list
-   * @param mc   MC list
-   * @param dir  Output directory
+   * Process a single eta bin 
    * 
-   * @return Stack of histograms on success
-   */  
-  THStack* UnfoldEtaBin(UInt_t       method, 
-                       Double_t&    regParam,
-                       TCollection* real, 
-                       TCollection* mc, 
-                       TDirectory* dir)
+   * @param measured     Input collection of measured data
+   * @param corrections  Input collection of correction data
+   * @param method       Unfolding method to use 
+   * @param regParam     Regularisation parameter
+   * @param out          Output directory. 
+   *
+   * @return Stack of histograms or null 
+   */
+  THStack* ProcessBin(TCollection* measured, 
+                     TCollection* corrections, 
+                     UInt_t       method,
+                     Double_t     regParam, 
+                     TDirectory*  out)
   {
-    TH1*  realRaw  = GetH1(real, "rawDist");
-    TH1*  mcRaw    = GetH1(mc,   "rawDist");
-    TH1*  mcTruth  = GetH1(mc,   "truth");
-    TH2*  response = GetH2(mc,   "response");
-    
-    if (!realRaw) { 
-      Warning("UnfoldEtaBin", "Real raw distribution not found in %s", 
-             real->GetName());
+    Printf("   Processing %s ...", measured->GetName());
+    // Try to get the data 
+    TH1* inRaw    = GetH1(measured, "rawDist");
+    TH1* inTruth  = GetH1(corrections, "truth");
+    TH1* inTruthA = GetH1(corrections, "truthAccepted");
+    TH1* inTrgVtx = GetH1(corrections, "triggerVertex");
+    TH2* inResp   = GetH2(corrections, "response");
+    if (!inRaw || !inTruth || !inTruthA || !inTrgVtx || !inResp) 
       return 0;
-    }
-    if (!mcRaw) { 
-      Warning("UnfoldEtaBin", "MC raw distribution not found in %s", 
-             mc->GetName());
-      return 0;
-    }
-    if (!mcTruth) { 
-      Warning("UnfoldEtaBin", "MC true distribution not found in %s", 
-             mc->GetName());
-      return 0;
-    }
-    if (!response) { 
-      Warning("UnfoldEtaBin", "Response matrix not found in %s", 
-             mc->GetName());
-      return 0;
-    }
     
-    Int_t    mN = realRaw->GetNbinsX();
-    Double_t mL = realRaw->GetXaxis()->GetXmin();
-    Double_t mH = realRaw->GetXaxis()->GetXmax();
-    Int_t    tN = mcRaw->GetNbinsX();
-    Double_t tL = mcRaw->GetXaxis()->GetXmin();
-    Double_t tH = mcRaw->GetXaxis()->GetXmax();
+    // Make output directory
+    TDirectory* dir = out->mkdir(measured->GetName());
+    dir->cd();
+
+    // Copy the input to the output 
+    TH1* outRaw    = static_cast<TH1*>(inRaw    ->Clone("measured"));
+    TH1* outTruth  = static_cast<TH1*>(inTruth  ->Clone("truth"));
+    TH1* outTruthA = static_cast<TH1*>(inTruthA ->Clone("truthAccepted"));
+    TH1* outTrgVtx = static_cast<TH1*>(inTrgVtx ->Clone("triggerVertex"));
+    TH2* outResp   = static_cast<TH2*>(inResp   ->Clone("response"));
+
+    // Make our response matrix 
+    RooUnfoldResponse matrix(0, 0, inResp);
     
-    RooUnfoldResponse matrix(mN, mL, mH, tN, tL, tH);
-    for (Int_t i = 1; i <= mN; i++) { 
-      Double_t mX = response->GetYaxis()->GetBinCenter(i);
-      for (Int_t j = 1; j <= tN; j++) { 
-       Double_t tX = response->GetXaxis()->GetBinCenter(j);
-       matrix.Fill(mX, tX, response->GetBinContent(j, i));
+    // Store regularization parameter 
+    Double_t             r        = regParam;
+    RooUnfold::Algorithm algo     = (RooUnfold::Algorithm)method;
+    RooUnfold*           unfolder = RooUnfold::New(algo, &matrix, inRaw, r);
+    unfolder->SetVerbose(1);
+
+    // Do the unfolding and get the result
+    TH1* res = unfolder->Hreco();
+    res->SetDirectory(0);
+
+    // Make a copy to store on the output 
+    TH1* outUnfold = static_cast<TH1*>(res->Clone("unfolded"));
+    TString tit(outUnfold->GetTitle());
+    tit.ReplaceAll("Unfold Reponse matrix", "Unfolded P(#it{N}_{ch})");
+    outUnfold->SetTitle(tit);
+
+    // Clone the unfolded results and divide by the trigger/vertex
+    // bias correction
+    TH1* outCorr   = static_cast<TH1*>(outUnfold->Clone("corrected"));
+    outCorr->Divide(inTrgVtx);
+    tit.ReplaceAll("Unfolded", "Corrected");
+    outCorr->SetTitle(tit);
+
+    // Now normalize the output to integral=1 
+    TH1*  hists[] = { outRaw, outUnfold, outCorr, 0 };
+    TH1** phist   = hists;
+    while (*phist) { 
+      TH1* h = *phist;
+      if (h) { 
+       Double_t intg = h->Integral(1, h->GetXaxis()->GetXmax());
+       h->Scale(1. / intg, "width");
       }
+      phist++;
     }
-
-    Double_t regParm = regParam;
-    RooUnfold::Algorithm alg = (RooUnfold::Algorithm)method;
-    RooUnfold* realUnfold = RooUnfold::New(alg, &matrix, realRaw, regParm);
-    realUnfold->SetVerbose(0);
-    TH1* resReal = realUnfold->Hreco();
-    TH1* outReal = static_cast<TH1*>(resReal->Clone("realUnfolded"));
-    resReal->SetDirectory(0);
-    regParam = realUnfold->GetRegParm();
-    // Info("UnfoldEtaBin", "Used regularization parameter: %f", regParam);
-    delete resReal;
-    delete realUnfold;
-
-    RooUnfold* mcUnfold = RooUnfold::New(alg, &matrix, mcRaw, regParm);
-    mcUnfold->SetVerbose(0);
-    TH1* resMC = mcUnfold->Hreco();
-    TH1* outMC = static_cast<TH1*>(resMC->Clone("mcUnfolded"));
-    resMC->SetDirectory(0);
-    delete resMC;
-    delete mcUnfold;
-
-    TH1* outRealRaw  = static_cast<TH1*>(realRaw->Clone("realRaw"));
-    TH1* outMcRaw    = static_cast<TH1*>(mcRaw->Clone("mcRaw"));
-    TH1* outMcTruth  = static_cast<TH1*>(mcTruth->Clone("mcTruth"));
-    TH2* outResponse = static_cast<TH2*>(matrix.Hresponse()->Clone("response"));
-
-    outRealRaw ->SetTitle("Measured");
-    outMcRaw   ->SetTitle("Measured (MC)");
-    outMcTruth ->SetTitle("Truth (MC)");
-    outResponse->SetTitle("Response matrix");
-    outReal    ->SetTitle("Unfolded");
-    outMC      ->SetTitle("Unfolded (MC)");
-
-    outRealRaw ->SetDirectory(dir);
-    outMcRaw   ->SetDirectory(dir);
-    outMcTruth ->SetDirectory(dir);
-    outResponse->SetDirectory(dir);
-    outReal    ->SetDirectory(dir);
-    outMC      ->SetDirectory(dir);
     
-    outMcRaw  ->SetMarkerStyle(24); // Open circle
-    outRealRaw->SetMarkerStyle(20); // Circle
-    outMcTruth->SetMarkerStyle(24); // Open Circle
-    outMC     ->SetMarkerStyle(24); // Open Circle
-    outReal   ->SetMarkerStyle(20); // Circle
+    // And make ratios
+    TH1* ratioTrue = static_cast<TH1*>(outCorr->Clone("ratioCorrTruth"));
+    tit = ratioTrue->GetTitle();
+    tit.ReplaceAll("Corrected", "Corrected/MC 'truth'");
+    ratioTrue->SetTitle(tit);
+    ratioTrue->Divide(outTruth);
+    ratioTrue->SetYTitle("P_{corrected}(#it{N}_{ch})/P_{truth}(#it{N}_{ch})");
+
+    TH1* ratioAcc  = static_cast<TH1*>(outUnfold->Clone("ratioUnfAcc"));
+    tit = ratioAcc->GetTitle();
+    tit.ReplaceAll("Unfolded", "Unfolded/MC selected");
+    ratioAcc->SetTitle(tit);
+    ratioAcc->Divide(outTruthA);
+    ratioAcc->SetYTitle("P_{unfolded}(#it{N}_{ch})/P_{MC}(#it{N}_{ch})");
     
-    outMcRaw  ->SetMarkerSize(1.2); // Open circle
-    outMC     ->SetMarkerSize(1.2); // Open circle
-    outMcTruth->SetMarkerSize(1.6); // Open square
-
-    outMcRaw  ->SetMarkerColor(kMeasuredColor);
-    outRealRaw->SetMarkerColor(kMeasuredColor);
-    outMcTruth->SetMarkerColor(kTruthColor);
-    outMC     ->SetMarkerColor(kUnfoldedColor);
-    outReal   ->SetMarkerColor(kUnfoldedColor);
-
-    outMcRaw  ->SetFillColor(kSysColor);
-    outRealRaw->SetFillColor(kSysColor);
-    outMcTruth->SetFillColor(kSysColor);
-    outMC     ->SetFillColor(kSysColor);
-    outReal   ->SetFillColor(kSysColor);
-
-    outMcRaw  ->SetFillStyle(1001);
-    outRealRaw->SetFillStyle(0);
-    outMcTruth->SetFillStyle(1001);
-    outMC     ->SetFillStyle(0);
-    outReal   ->SetFillStyle(0);
 
-    outMcRaw  ->SetLineColor(kBlack);
-    outRealRaw->SetLineColor(kBlack);
-    outMcTruth->SetLineColor(kBlack);
-    outMC     ->SetLineColor(kBlack);
-    outReal   ->SetLineColor(kBlack);
-
-    TString tit(real->GetName());
+    // Make a stack 
+    tit = measured->GetName();
     tit.ReplaceAll("m", "-");
     tit.ReplaceAll("p", "+");
     tit.ReplaceAll("d", ".");
     tit.ReplaceAll("_", "<#it{#eta}<");
-    
     THStack* stack = new THStack("all", tit);
-    stack->Add(outMcRaw,      "E2");
-    stack->Add(outRealRaw,    "E1");
-    stack->Add(outMcTruth,    "E2");
-    stack->Add(outMC,         "E1");
-    stack->Add(outReal,       "E1");
-
+    stack->Add(outTruth,  "E2");
+    stack->Add(outTruthA, "E2");
+    stack->Add(outRaw,    "E1");
+    stack->Add(outUnfold, "E1");
+    stack->Add(outCorr,   "E1");
     dir->Add(stack);
 
+    // Rest of the function is devoted to making the output look nice 
+    outRaw   ->SetDirectory(dir); 
+    outTruth ->SetDirectory(dir);  
+    outTruthA->SetDirectory(dir);  
+    outTrgVtx->SetDirectory(dir);  
+    outResp  ->SetDirectory(dir);  
+    outUnfold->SetDirectory(dir);   
+    outCorr  ->SetDirectory(dir); 
+
+    outRaw   ->SetMarkerStyle(20);  // Measured is closed
+    outTruth ->SetMarkerStyle(24);  // MC is open
+    outTruthA->SetMarkerStyle(24);  // MC is open
+    outTrgVtx->SetMarkerStyle(20);  // Derived is closed
+    outUnfold->SetMarkerStyle(20);  // Derived is closed   
+    outCorr  ->SetMarkerStyle(20);  // Derived is closed 
+
+    outRaw   ->SetMarkerSize(0.9); 
+    outTruth ->SetMarkerSize(1.6);  
+    outTruthA->SetMarkerSize(1.4);  
+    outTrgVtx->SetMarkerSize(1.0);  
+    outUnfold->SetMarkerSize(0.9);   
+    outCorr  ->SetMarkerSize(1.0);
+    outRaw   ->SetMarkerColor(kColorMeasured); 
+    outTruth ->SetMarkerColor(kColorTruth);  
+    outTruthA->SetMarkerColor(kColorAccepted);  
+    outTrgVtx->SetMarkerColor(kColorTrgVtx);  
+    outUnfold->SetMarkerColor(kColorUnfolded);   
+    outCorr  ->SetMarkerColor(kColorCorrected); 
+
+    outRaw   ->SetFillColor(kColorError);     
+    outTruth ->SetFillColor(kColorError);  
+    outTruthA->SetFillColor(kColorError);  
+    outTrgVtx->SetFillColor(kColorError);  
+    outUnfold->SetFillColor(kColorError);   
+    outCorr  ->SetFillColor(kColorError); 
+
+    outRaw   ->SetFillStyle(0); 
+    outTruth ->SetFillStyle(1001);  
+    outTruthA->SetFillStyle(1001);  
+    outTrgVtx->SetFillStyle(0);  
+    outUnfold->SetFillStyle(0);   
+    outCorr  ->SetFillStyle(0);
+
+    outRaw   ->SetLineColor(kBlack); 
+    outTruth ->SetLineColor(kBlack);  
+    outTruthA->SetLineColor(kBlack);  
+    outTrgVtx->SetLineColor(kBlack);  
+    outUnfold->SetLineColor(kBlack);   
+    outCorr  ->SetLineColor(kBlack); 
+
+    // Legend 
+    TLegend* l = StackLegend(stack);
+    l->AddEntry(outRaw,     "Raw",                 "lp");
+    l->AddEntry(outTruth,   "MC 'truth'",          "fp");
+    l->AddEntry(outTruthA,  "MC 'truth' accepted", "fp");
+    l->AddEntry(outUnfold,  "Unfolded",            "lp");
+    l->AddEntry(outCorr,    "Corrected",           "lp");
+
     return stack;
   }
+  static void BinAttributes(Int_t i,
+                           Int_t&    open, 
+                           Int_t&    closed, 
+                           Float_t&  size,
+                           Double_t& factor) 
+  {
+    // --- Setup for markers -----------------------------------------
+    const Int_t   nMarkers = 7;
+    const Int_t   aClosed[] = { 20,  21,  22,  23,  29,  33,  34  };
+    const Int_t   aOpen[]   = { 24,  25,  26,  32,  30,  27,  28  };
+    const Float_t aSize[]   = { 1.1, 1.0, 1.2, 1.2, 1.2, 1.2, 1.0 };
+    Int_t         j         = i % nMarkers;
+
+    open   = aOpen[j];
+    closed = aClosed[j];
+    size   = aSize[j];
+    factor = TMath::Power(10, i);
+  }
   /** 
-   * Drwa output 
+   * Add the bin histograms to our summary stacks 
    * 
-   * @param which 
+   * @param bin       Bin stack
+   * @param i         Current off-set in the stacks 
+   * @param measured  All measured @f$ P(N_{ch})@f$ 
+   * @param truth     All MC truth @f$ P(N_{ch})@f$ 
+   * @param accepted  All MC accepted @f$ P(N_{ch})@f$ 
+   * @param unfolded  All unfolded @f$ P(N_{ch})@f$ 
+   * @param corrected All corrected @f$ P(N_{ch})@f$ 
    */
-  void DrawType(const char* which="symmetric") 
+  void Bin2Stack(const THStack* bin, Int_t i, 
+                THStack* measured, 
+                THStack* truth, 
+                THStack* accepted, 
+                THStack* unfolded,
+                THStack* corrected)
   {
-    const char* fname = "forward_unfolded.root";
-    TFile*   outF = 0;
-    TObject* outO = gROOT->GetListOfFiles()->FindObject(fname);
-    if (outO) outF = static_cast<TFile*>(outO);
-    else      outF = TFile::Open(fname, "READ");
-    
-    if (!outF) { 
-      Warning("DrawType", "Failed to open file %s", fname);
-      return;
-    }
-    
-    TNamed* method = 0;
-    outF->GetObject("method", method);
-    TParameter<double>* regParam = 0;
-    outF->GetObject("regParam", regParam);
-    if (!method) Warning("DrawType", "Didn't find the method string");
-    if (!regParam) 
-      Warning("DrawType", "Didn't find the regularization parameter");
-    
-
-    THStack* stack = 0;
-    outF->GetObject(Form("/%s/all", which), stack);
-    if (!stack) { 
-      Warning("DrawType", "Failed to get /%s/all from file %s", which, fname);
-      return;
-    }
-    
-    TList* leg = 0;
-    outF->GetObject(Form("/%s/legend", which), leg);
-    if (!leg) { 
-      Warning("DrawType", "Failed to get /%s/legend from file %s", 
-             which, fname);
-      return;
-    }
+    Int_t open, closed;
+    Double_t factor; 
+    Float_t  size;
+    BinAttributes(i, open, closed, size, factor);
 
-    TMultiGraph* other = 0;
-    outF->GetObject(Form("/%s/other", which), other);
-    // if (!other) Warning("DrawType", "No other data found for %s", which);
-
-    THStack* ratios = 0;
-    outF->GetObject(Form("/%s/ratios", which), ratios);
-    // if (!ratios) Warning("DrawType", "No ratios data found for %s", which);
-
-    TCanvas* c = new TCanvas(which, Form("%s P(#it{N}_{ch})", which));
-    c->SetLogy();
-    c->SetTopMargin(0.01);
-    c->SetRightMargin(0.02);
-    c->cd();
-    
-    THStack* drawn = static_cast<THStack*>(stack->Clone(which));
-    drawn->Draw("nostack");
-    drawn->GetXaxis()->SetTitle("#it{#eta}");
-    drawn->GetYaxis()->SetTitle("P(#it{N}_{ch})");
-    drawn->GetHistogram()->SetDirectory(0);
-    TIter nextH(drawn->GetHists());
+    TIter next(bin->GetHists());
     TH1*  h = 0;
-    while ((h = static_cast<TH1*>(nextH()))) h->SetDirectory(0);
-
-    Bool_t        hasCMS   = false;
-    Bool_t        hasALICE = false;
-    if (other) {
-      TIter nextG(other->GetListOfGraphs());
-      TGraph* g = 0;
-      while ((g = static_cast<TGraph*>(nextG()))) {
-       g->Draw("p same");
-       if (g->GetMarkerColor() == kALICEColor) hasALICE = true;
-       if (g->GetMarkerColor() == kCMSColor)   hasCMS   = true;
+    while ((h = static_cast<TH1*>(next()))) {
+      THStack* tmp = 0;
+      Int_t    col = h->GetMarkerColor();
+      Int_t    sty = 0;
+      switch (col) { 
+      case kColorMeasured:  tmp = measured;   sty = closed;  break;
+      case kColorTruth:     tmp = truth;      sty = open;    break;
+      case kColorAccepted:  tmp = accepted;   sty = open;    break;
+      case kColorUnfolded:  tmp = unfolded;   sty = closed;  break;
+      case kColorCorrected: tmp = corrected;  sty = closed;  break;
+      default: continue; 
       }
-    }
-    TLegend* l = new TLegend(.65, .75, .975, .975);
-    l->SetBorderSize(0);
-    l->SetFillColor(0);
-    l->SetFillStyle(0);
-    
-    TIter         next(leg);
-    TObject*      o = 0;
-    TLegendEntry* e = 0;
-    Int_t         d = 0;
-    while ((o = next())) { 
-      e = l->AddEntry(Form("dummy%02d", d++), o->GetName(), "p");
-      e->SetMarkerStyle(o->GetUniqueID());
-      e->SetLineColor(kBlack);
-      e->SetMarkerSize(1.3);
-    }
-    l->Draw();
+      // Now clone, and add to the appropriate stack 
+      TH1* cln = static_cast<TH1*>(h->Clone(h->GetName()));
+      cln->SetDirectory(0);
+      cln->SetMarkerStyle(sty);
+      cln->SetMarkerSize(size);
+      cln->Scale(factor); // Scale by 10^i
 
-    Double_t x1 = other ? .12 : .7;
-    Double_t y1 = .1;
-    Double_t w1 = other ? .20 : .15;
-    Double_t h1 = other ? .3 : .25;
-    TLegend* l2 = new TLegend(x1, y1, x1+w1, y1+h1);
-    l2->SetBorderSize(0);
-    l2->SetFillColor(0);
-    l2->SetFillStyle(0);
+      // Make sure we do not get the old legend 
+      TObject* tst = cln->FindObject("legend");
+      if (tst) cln->GetListOfFunctions()->Remove(tst);
 
-    const char* oth[] = { "Measured", 
-                         "Truth", 
-                         "Unfolded", 
-                         "CMS", 
-                         "ALICE" };
-    Int_t       col[] = { kMeasuredColor, 
-                         kTruthColor, 
-                         kUnfoldedColor, 
-                         kCMSColor, 
-                         kALICEColor };
-    for (Int_t i = 0; i < 5; i++) {
-      if (i == 3 && !hasCMS)   continue;
-      if (i == 4 && !hasALICE) continue;
-      e = l2->AddEntry(Form("dummy%02d", d++), oth[i], "f");
-      e->SetFillColor(col[i]);
-      e->SetFillStyle(1001);
+      tmp->Add(cln, next.GetOption());
     }
-    e = l2->AddEntry(Form("dummy%02d", d++), "Real", "p");
-    e->SetMarkerStyle(20);
-    e->SetMarkerSize(1.3);
-    e->SetMarkerColor(kBlack);
-    e = l2->AddEntry(Form("dummy%02d", d++), "MC", "p");
-    e->SetMarkerStyle(24);
-    e->SetMarkerSize(1.3);
-    e->SetMarkerColor(kBlack);
-    l2->Draw();
-
-
-    TString post;
-    if (method && regParam) {
-      TString mes = method->GetTitle();
-      if      (method->GetUniqueID() == RooUnfold::kBayes || 
-              method->GetUniqueID() == RooUnfold::kDagostini)
-       mes.Append(Form(" (N_{iter}=%d)", int(regParam->GetVal())));
-      else if (method->GetUniqueID() == RooUnfold::kSVD) 
-       mes.Append(Form(" (k=%d)", int(regParam->GetVal())));
-      else if (method->GetUniqueID() == RooUnfold::kTUnfold) 
-       mes.Append(Form(" (#tau=%f)", regParam->GetVal()));
+    
+    // Add entries to our stacks 
+    TString   txt      = bin->GetTitle();
+    if      (i == 0) txt.Append(" (#times1)");
+    else if (i == 1) txt.Append(" (#times10)");
+    else             txt.Append(Form(" (#times10^{%d})", i));
+    THStack*  stacks[] = { measured, truth, accepted, unfolded, corrected, 0 };
+    THStack** pstack   = stacks;
+    while (*pstack) { 
+      TLegend* leg = StackLegend(*pstack);
+      pstack++;
+      if (!leg) continue;
       
-      TLatex* meth = new TLatex(.97, .749, mes);
-      meth->SetNDC();
-      meth->SetTextAlign(33);
-      meth->SetTextFont(42);
-      meth->SetTextSize(.03);
-      meth->Draw();
-
-      post.Form("_%s", method->GetTitle());
-      post.ToLower();
-    }
-    c->Modified();
-    c->Update();
-    c->cd();
-    c->SaveAs(Form("pnch_%s%s.pdf", which, post.Data()));
-
-    if (ratios) {
-      c = new TCanvas(Form("r%s", which), 
-                     Form("%s P(#it{N}_{ch})", which));
-      c->SetTopMargin(0.01);
-      c->SetRightMargin(0.02);
-      c->cd();
-      ratios->Draw("nostack");
-      c->BuildLegend();
+      TObject* dummy = 0;
+      TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
+      e->SetMarkerStyle(closed);
+      e->SetMarkerSize(1.2*size);
+      e->SetMarkerColor(kBlack);
+      e->SetFillColor(0);
+      e->SetFillStyle(0);
+      e->SetLineColor(kBlack);
     }
-
-    outF->cd();
   }
-  void DrawAll()
-  {
-    DrawType("symmetric");
-    DrawType("negative");
-    DrawType("positive");
-    DrawType("other");
-  }
-  TGraphAsymmErrors* GetOther(UShort_t type, Double_t eta, UShort_t sNN,
-                             Int_t factor, Int_t marker)
+  /** 
+   * Add distributions from other experiments to stacks 
+   * 
+   * @param name     Name of current bin 
+   * @param i        Index of current bin
+   * @param sNN      Center of mass energy 
+   * @param allALICE Stack of ALICE data 
+   * @param allCMS   Stack of CMS data 
+   */
+  void Other2Stack(const TString& name, Int_t i,
+                  UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS) 
   {
-    TGraphAsymmErrors* g = 0;
-    Int_t color = kBlack;
-    switch (type) { 
-    case 0: g = GetCMS(eta, sNN);   color = kCMSColor; break;
-    case 1: g = GetALICE(eta, sNN); color = kALICEColor; break;
-    }
-    if (!g) {
-      // Warning("GetOther", "No other data found for type=%d, eta=%f, sNN=%d", 
-      //         type, eta, sNN);
-      return g;
+    if (!allALICE && !allCMS) return;
+
+    TString tmp(name);
+    tmp.ReplaceAll("p", "+");
+    tmp.ReplaceAll("m", "-");
+    tmp.ReplaceAll("_", " ");
+    TObjArray* tokens = tmp.Tokenize(" ");
+    if (!tokens || tokens->GetEntriesFast() < 2) { 
+      Error("Other2Stack", "Failed to decode eta range from %s", name.Data());
+      if (tokens) tokens->Delete();
+      return;
     }
+    Double_t eta1 = static_cast<TObjString*>(tokens->At(0))->String().Atof();
+    Double_t eta2 = static_cast<TObjString*>(tokens->At(1))->String().Atof();
+    tokens->Delete();
     
-    g->SetMarkerStyle(marker);
-    g->SetMarkerColor(color);
-    for (Int_t j = 0; j < g->GetN(); j++) { 
-      g->SetPoint(j, g->GetX()[j], g->GetY()[j]*factor);
-      g->SetPointError(j, g->GetEXlow()[j], g->GetEXhigh()[j], 
-                      g->GetEYlow()[j]*factor, g->GetEYhigh()[j]*factor);
-    }
-    return g;
-  }    
-
-  TGraphAsymmErrors* GetCMS(Double_t eta, UShort_t sNN)
-  {
-    Double_t aEta = TMath::Abs(eta);
-    TGraphAsymmErrors* g = 0;
-    switch (sNN) { 
-    case 900: 
-      if      (aEta <= 0.51)  g = CMSsqrts900eta05();
-      else if (aEta <= 1.01)  g = CMSsqrts900eta10();
-      else if (aEta <= 1.51)  g = CMSsqrts900eta15();
-      else if (aEta <= 2.01)  g = CMSsqrts900eta20();
-      else if (aEta <= 2.41)  g = CMSsqrts900eta24();
-      break;
-    default: 
-      return g;
+    if (TMath::Abs(eta2-eta1) > 1e3) 
+      // Not symmetric bin 
+      return;
+    Double_t aEta = TMath::Abs(eta1);
+
+    Int_t open, closed;
+    Double_t factor; 
+    Float_t  size;
+    BinAttributes(i, open, closed, size, factor);
+
+    if (allALICE) {
+      TGraphAsymmErrors* g = GetOther(1, aEta, sNN, factor);
+      if (g) {
+       g->SetMarkerStyle(closed);
+       g->SetMarkerColor(kColorALICE);
+       g->SetMarkerSize(size);
+       allALICE->Add(g, "p same");
+      }
     }
-    if (g) {
-      g->SetMarkerColor(kCMSColor);
-      g->SetLineColor(kBlack);
+    if (allCMS) {
+      TGraphAsymmErrors* g = GetOther(0, aEta, sNN, factor);
+      if (g) {
+       g->SetMarkerStyle(closed);
+       g->SetMarkerColor(kColorCMS);
+       g->SetMarkerSize(size);
+       allCMS->Add(g, "p same");
+      }
     }
-    return g;
+
+    
   }
   /** 
-   * CMS pp, @f$\sqrt{s} = 900GeV@f$, NSD, @f$|\eta|\lt0.5@f$ 
+   * 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
+   * stack.
    * 
-   * @return 
+   * @param stack Stack to get the legend from/modify 
+   * 
+   * @return The legend object or null
    */
-  TGraphAsymmErrors* CMSsqrts900eta05()
-  {
-    // Plot: p8068_d2x1y1
-    double x[] = { 0.0,        1.0,
-                  2.0, 3.0,
-                  4.0, 5.0,
-                  6.0, 7.0,
-                  8.0, 9.0,
-                  10.0,        11.0,
-                  12.0,        13.0,
-                  14.0,        15.0,
-                  16.0,        17.0,
-                  18.0,        19.0,
-                  20.0,        22.0,
-                  26.0 };
-    double xem[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.5,
-                    2.5 };
-    double xep[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.5,
-                    2.5 };
-    double y[] = { 0.193,      0.1504,
-                  0.13507,     0.11647,
-                  0.09566,     0.07567,
-                  0.05817,     0.04427,
-                  0.0337,      0.02572,
-                  0.01948,     0.01457,
-                  0.01069,     0.00769,
-                  0.0055,      0.00396,
-                  0.00289,     0.002112,
-                  0.001544,    0.001118,
-                  7.52E-4,     3.35E-4,
-                  9.8E-5 };
-    double yem[] = { 0.02189155088156159,      0.012180722474467595,
-                    0.007909551188278638,      0.006072832946821442,
-                    0.004838605170914444,      0.003828276374558138,
-                    0.003016836091006603,      0.002365797962633327,
-                    0.0018847015678881366,     0.0015610893632332517,
-                    0.001308166656049603,      0.0011140017953306899,
-                    9.219544457292888E-4,      7.203471385380802E-4,
-                    5.508175741568165E-4,      4.3011626335213136E-4,
-                    3.3837848631377264E-4,     2.746961958236772E-4,
-                    2.2377220560203628E-4,     1.9070920271449933E-4,
-                    1.6011246047700347E-4,     8.134494452638099E-5,
-                    2.5059928172283337E-5 };
-    double yep[] = { 0.031663227883461285,     0.015861904047118684,
-                    0.010948538715280683,      0.008490583018850945,
-                    0.007046453008429135,      0.005675746646917919,
-                    0.004395201929377079,      0.003365070578754627,
-                    0.0026049952015310893,     0.0020535335400231475,
-                    0.0016522106403240478,     0.0013497407158413798,
-                    0.001118033988749895,      9.060905032059435E-4,
-                    7.061161377563892E-4,      5.457105459856901E-4,
-                    4.148493702538308E-4,      3.257759966602819E-4,
-                    2.6689698387205504E-4,     2.292705825002414E-4,
-                    1.904205871222962E-4,      9.646242791885347E-5,
-                    3.417601498127012E-5 };
-    double ysm[] = { 0.002,    0.0014,
-                    7.9E-4,    7.2E-4,
-                    6.1E-4,    5.4E-4,
-                    4.7E-4,    4.1E-4,
-                    3.6E-4,    3.1E-4,
-                    2.7E-4,    2.3E-4,
-                    2.0E-4,    1.7E-4,
-                    1.5E-4,    1.3E-4,
-                    1.1E-4,    9.7E-5,
-                    8.5E-5,    7.1E-5,
-                    5.6E-5,    2.9E-5,
-                    1.2E-5 };
-    double ysp[] = { 0.002,    0.0014,
-                    7.9E-4,    7.2E-4,
-                    6.1E-4,    5.4E-4,
-                    4.7E-4,    4.1E-4,
-                    3.6E-4,    3.1E-4,
-                    2.7E-4,    2.3E-4,
-                    2.0E-4,    1.7E-4,
-                    1.5E-4,    1.3E-4,
-                    1.1E-4,    9.7E-5,
-                    8.5E-5,    7.1E-5,
-                    5.6E-5,    2.9E-5,
-                    1.2E-5 };
-    int np = 23;
-    TGraphAsymmErrors* g = new TGraphAsymmErrors(np, 
-                                                x, y, 
-                                                xem, 
-                                                xep, 
-                                                yem, 
-                                                yep);
-    g->SetName("/HepData/8068/d2x1y1");
-    g->SetTitle("/HepData/8068/d2x1y1");
-
-    return g;
-  }
-  TGraphAsymmErrors* CMSsqrts900eta10()
+  TLegend* StackLegend(THStack* stack) 
   {
+    TList* hists = stack->GetHists();
+    if (!hists) return 0;
+    
+    TObject* first = hists->First();
+    if (!first) return 0;
 
-    // Plot: p8068_d3x1y1
-    double x[] = { 0.0,        1.0,
-                  2.0, 3.0,
-                  4.0, 5.0,
-                  6.0, 7.0,
-                  8.0, 9.0,
-                  10.0,        11.0,
-                  12.0,        13.0,
-                  14.0,        15.0,
-                  16.0,        17.0,
-                  18.0,        19.0,
-                  20.0,        21.0,
-                  22.0,        23.0,
-                  24.0,        25.0,
-                  26.0,        27.0,
-                  28.0,        29.0,
-                  30.0,        31.0,
-                  32.0,        33.0,
-                  34.0,        36.0,
-                  40.0,        45.0,
-                  50.0,        55.0 };
-    double xem[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.5,
-                    2.5,       2.5,
-                    2.5,       2.5 };
-    double xep[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.5,
-                    2.5,       2.5,
-                    2.5,       2.5 };
-    double y[] = { 0.1076,     0.0655,
-                  0.07574,     0.08168,
-                  0.0813,      0.07622,
-                  0.06843,     0.06034,
-                  0.05259,     0.0458,
-                  0.03987,     0.03472,
-                  0.03035,     0.02661,
-                  0.02328,     0.02018,
-                  0.01728,     0.0147,
-                  0.0125,      0.01067,
-                  0.00916,     0.00786,
-                  0.00671,     0.0057,
-                  0.00478,     0.00394,
-                  0.00319,     0.00251,
-                  0.00196,     0.00156,
-                  0.001301,    0.001121,
-                  9.84E-4,     8.25E-4,
-                  6.39E-4,     3.78E-4,
-                  1.96E-4,     6.0E-5,
-                  1.25E-5,     1.6E-6 };
-    double yem[] = { 0.015802847844613326,     0.009464142856064674,
-                    0.008214195030555337,      0.00676426640516176,
-                    0.0058916636020736966,     0.005200240379059414,
-                    0.004441733895676326,      0.003671511950137164,
-                    0.0029149442533262966,     0.0024150569351466646,
-                    0.0020554318281081475,     0.001767059704707229,
-                    0.0015500322577288513,     0.0013800362314084368,
-                    0.0012684636376341263,     0.0011763077828527702,
-                    0.001084158659975559,      9.82344135219425E-4,
-                    8.805679985100526E-4,      7.692203845452875E-4,
-                    6.705221845696084E-4,      5.913543776789007E-4,
-                    5.375872022286245E-4,      5.154609587543949E-4,
-                    5.124451190127582E-4,      4.904079934095691E-4,
-                    4.4922154890432404E-4,     3.7E-4,
-                    3.0083217912982643E-4,     2.4166091947189146E-4,
-                    2.1108529081866412E-4,     1.8252944967867517E-4,
-                    1.662077013859466E-4,      1.4663219291819924E-4,
-                    1.2229881438509532E-4,     7.826237921249263E-5,
-                    4.74236228055175E-5,       2.1095023109728988E-5,
-                    6.161980201201558E-6,      1.3601470508735443E-6 };
-    double yep[] = { 0.02725949375905576,      0.011154371340420759,
-                    0.009740395269186974,      0.008248078564126314,
-                    0.007285581651453781,      0.006454339625399333,
-                    0.005595158621522718,      0.004804039966528172,
-                    0.00405504623894722,       0.0034544174617437305,
-                    0.0029644561052577585,     0.002545584412271571,
-                    0.0021980445855350615,     0.0019087430418995638,
-                    0.0016690416411821486,     0.001468911161370898,
-                    0.0013179529581893276,     0.0011866338946785568,
-                    0.001065129100156408,      9.339164844888434E-4,
-                    8.052328855678958E-4,      6.868041933477111E-4,
-                    5.94810894318522E-4,       5.440588203494177E-4,
-                    5.220153254455275E-4,      5.192301994298867E-4,
-                    4.972926703662542E-4,      4.368065933568311E-4,
-                    3.5735136770411276E-4,     2.8792360097775935E-4,
-                    2.4200206610688264E-4,     2.0310588371585893E-4,
-                    1.670748335327616E-4,      1.364587849865299E-4,
-                    1.1970797801316335E-4,     9.470480452437458E-5,
-                    5.758472019555187E-5,      2.195449840010015E-5,
-                    7.481310045707236E-6,      1.7029386365926402E-6 };
-    double ysm[] = { 0.0018,   0.0011,
-                    6.3E-4,    6.8E-4,
-                    6.1E-4,    5.6E-4,
-                    5.3E-4,    4.8E-4,
-                    4.5E-4,    4.1E-4,
-                    3.8E-4,    3.6E-4,
-                    3.5E-4,    3.3E-4,
-                    3.1E-4,    2.9E-4,
-                    2.7E-4,    2.5E-4,
-                    2.3E-4,    2.1E-4,
-                    2.0E-4,    1.9E-4,
-                    1.7E-4,    1.6E-4,
-                    1.5E-4,    1.4E-4,
-                    1.3E-4,    1.2E-4,
-                    1.1E-4,    1.0E-4,
-                    9.4E-5,    8.6E-5,
-                    8.3E-5,    7.5E-5,
-                    6.1E-5,    3.5E-5,
-                    2.0E-5,    1.1E-5,
-                    4.6E-6,    1.1E-6 };
-    double ysp[] = { 0.0018,   0.0011,
-                    6.3E-4,    6.8E-4,
-                    6.1E-4,    5.6E-4,
-                    5.3E-4,    4.8E-4,
-                    4.5E-4,    4.1E-4,
-                    3.8E-4,    3.6E-4,
-                    3.5E-4,    3.3E-4,
-                    3.1E-4,    2.9E-4,
-                    2.7E-4,    2.5E-4,
-                    2.3E-4,    2.1E-4,
-                    2.0E-4,    1.9E-4,
-                    1.7E-4,    1.6E-4,
-                    1.5E-4,    1.4E-4,
-                    1.3E-4,    1.2E-4,
-                    1.1E-4,    1.0E-4,
-                    9.4E-5,    8.6E-5,
-                    8.3E-5,    7.5E-5,
-                    6.1E-5,    3.5E-5,
-                    2.0E-5,    1.1E-5,
-                    4.6E-6,    1.1E-6 };
-    int np = 40;
-    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
-                                                x,
-                                                y,
-                                                xem,
-                                                xep,
-                                                yem,
-                                                yep);
-    g->SetName("/HepData/8068/d3x1y1");
-    g->SetTitle("/HepData/8068/d3x1y1");
-    return g;
-  }
-  TGraphAsymmErrors* CMSsqrts900eta15()
-  {
-    // Plot: p8068_d4x1y1
-    double x[] = { 0.0,        1.0,
-                  2.0, 3.0,
-                  4.0, 5.0,
-                  6.0, 7.0,
-                  8.0, 9.0,
-                  10.0,        11.0,
-                  12.0,        13.0,
-                  14.0,        15.0,
-                  16.0,        17.0,
-                  18.0,        19.0,
-                  20.0,        21.0,
-                  22.0,        23.0,
-                  24.0,        25.0,
-                  26.0,        27.0,
-                  28.0,        29.0,
-                  30.0,        31.0,
-                  32.0,        33.0,
-                  34.0,        35.0,
-                  36.0,        37.0,
-                  38.0,        39.0,
-                  40.0,        41.0,
-                  42.0,        43.0,
-                  44.0,        45.0,
-                  46.5,        49.0,
-                  53.0,        58.0,
-                  63.0,        68.0 };
-    double xem[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    1.0,       1.5,
-                    2.5,       2.5,
-                    2.5,       2.5 };
-    double xep[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    1.0,       1.5,
-                    2.5,       2.5,
-                    2.5,       2.5 };
-    double y[] = { 0.0753,     0.03557,
-                  0.04347,     0.05068,
-                  0.05566,     0.05786,
-                  0.05835,     0.05688,
-                  0.05348,     0.0489,
-                  0.04422,     0.04002,
-                  0.03647,     0.03337,
-                  0.03055,     0.02797,
-                  0.02556,     0.02329,
-                  0.02122,     0.01936,
-                  0.01765,     0.01604,
-                  0.01451,     0.01304,
-                  0.01168,     0.01045,
-                  0.00937,     0.0084,
-                  0.00751,     0.00667,
-                  0.00592,     0.00526,
-                  0.00466,     0.00413,
-                  0.00362,     0.00316,
-                  0.00272,     0.00233,
-                  0.00203,     0.00178,
-                  0.00156,     0.00136,
-                  0.0012,      0.001053,
-                  9.13E-4,     7.84E-4,
-                  5.85E-4,     3.73E-4,
-                  2.13E-4,     9.6E-5,
-                  3.6E-5,      7.4E-6 };
-    double yem[] = { 0.014524806367039803,     0.0061507154055443014,
-                    0.006900384047283165,      0.0068944397886992964,
-                    0.00601252858621063,       0.004686747272896204,
-                    0.0042880881520789655,     0.0039679339712248235,
-                    0.0035288099977187778,     0.0031197756329582422,
-                    0.0027127108213003464,     0.0023229722340140013,
-                    0.002033937068839643,      0.0017866449003649271,
-                    0.001591257364476281,      0.0014504137340772805,
-                    0.0013412307780542466,     0.0012490796611905903,
-                    0.0011594826432508594,     0.0010893117092916976,
-                    0.0010384603988597735,     9.972462083156796E-4,
-                    9.588013350011566E-4,      8.984430978086481E-4,
-                    8.160882305241266E-4,      7.178439941937245E-4,
-                    6.390618123468182E-4,      5.88727441181401E-4,
-                    5.571355310873647E-4,      5.255473337388365E-4,
-                    4.846648326421054E-4,      4.6615448083226656E-4,
-                    4.4384682042344296E-4,     4.3081318457076036E-4,
-                    4.0853396431630995E-4,     3.992492955535426E-4,
-                    3.676955262170047E-4,      3.2695565448543625E-4,
-                    2.9068883707497264E-4,     2.505992817228334E-4,
-                    2.1954498400100151E-4,     1.8601075237738274E-4,
-                    1.8601075237738274E-4,     1.7042300314218148E-4,
-                    1.7998055450520203E-4,     1.6058953888718904E-4,
-                    1.4985659811966905E-4,     8.561541917201597E-5,
-                    6.324555320336759E-5,      4.657252408878007E-5,
-                    2.4166091947189146E-5,     5.590169943749474E-6 };
-    double yep[] = { 0.03405304685340212,      0.00724989655098609,
-                    0.007817985674072318,      0.007751728839426725,
-                    0.006919566460407762,      0.00571192611997039,
-                    0.005232991496266739,      0.004813003220443552,
-                    0.00435332057170156,       0.0038838769290491172,
-                    0.003445692963686695,      0.0030549959083442323,
-                    0.0027252339349127445,     0.0024367396249907374,
-                    0.002189794510907359,      0.001969568480657629,
-                    0.0017906702655709678,     0.0016395731151735808,
-                    0.0015002999700059986,     0.0013710215169719256,
-                    0.0012712198865656563,     0.0011812705024675763,
-                    0.0011035397591387453,     0.0010139033484509261,
-                    9.121403400793104E-4,      8.132035415564789E-4,
-                    7.242237223399962E-4,      6.545991139621257E-4,
-                    6.135144660071188E-4,      5.724508712544685E-4,
-                    5.692099788303082E-4,      5.126402247190518E-4,
-                    4.531004303683677E-4,      4.4011362169330773E-4,
-                    4.1785164831552356E-4,     4.0853396431630995E-4,
-                    3.8626415831655934E-4,     3.54682957019364E-4,
-                    3.1780497164141406E-4,     2.7730849247724094E-4,
-                    2.459674775249769E-4,      2.3706539182259396E-4,
-                    2.1954498400100151E-4,     2.0068881383873889E-4,
-                    1.5864425612041553E-4,     1.6058953888718904E-4,
-                    1.4338758663147937E-4,     1.0245974819410792E-4,
-                    6.609841147864296E-5,      4.657252408878007E-5,
-                    2.5079872407968906E-5,     5.941380311005179E-6 };
-    double ysm[] = { 0.0019,   9.3E-4,
-                    5.3E-4,    5.8E-4,
-                    5.2E-4,    5.0E-4,
-                    4.9E-4,    4.7E-4,
-                    4.5E-4,    4.3E-4,
-                    4.2E-4,    3.9E-4,
-                    3.7E-4,    3.6E-4,
-                    3.6E-4,    3.4E-4,
-                    3.3E-4,    3.1E-4,
-                    3.0E-4,    2.9E-4,
-                    2.8E-4,    2.7E-4,
-                    2.7E-4,    2.6E-4,
-                    2.4E-4,    2.3E-4,
-                    2.2E-4,    2.1E-4,
-                    2.0E-4,    1.9E-4,
-                    1.8E-4,    1.8E-4,
-                    1.7E-4,    1.6E-4,
-                    1.5E-4,    1.5E-4,
-                    1.4E-4,    1.3E-4,
-                    1.3E-4,    1.2E-4,
-                    1.1E-4,    1.1E-4,
-                    1.1E-4,    1.0E-4,
-                    8.8E-5,    7.5E-5,
-                    5.6E-5,    3.3E-5,
-                    2.0E-5,    1.2E-5,
-                    1.0E-5,    4.1E-6 };
-    double ysp[] = { 0.0019,   9.3E-4,
-                    5.3E-4,    5.8E-4,
-                    5.2E-4,    5.0E-4,
-                    4.9E-4,    4.7E-4,
-                    4.5E-4,    4.3E-4,
-                    4.2E-4,    3.9E-4,
-                    3.7E-4,    3.6E-4,
-                    3.6E-4,    3.4E-4,
-                    3.3E-4,    3.1E-4,
-                    3.0E-4,    2.9E-4,
-                    2.8E-4,    2.7E-4,
-                    2.7E-4,    2.6E-4,
-                    2.4E-4,    2.3E-4,
-                    2.2E-4,    2.1E-4,
-                    2.0E-4,    1.9E-4,
-                    1.8E-4,    1.8E-4,
-                    1.7E-4,    1.6E-4,
-                    1.5E-4,    1.5E-4,
-                    1.4E-4,    1.3E-4,
-                    1.3E-4,    1.2E-4,
-                    1.1E-4,    1.1E-4,
-                    1.1E-4,    1.0E-4,
-                    8.8E-5,    7.5E-5,
-                    5.6E-5,    3.3E-5,
-                    2.0E-5,    1.2E-5,
-                    1.0E-5,    4.1E-6 };
-    int np = 52;
-    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
-                                                x,
-                                                y,
-                                                xem,
-                                                xep,
-                                                yem,
-                                                yep);
-    g->SetName("/HepData/8068/d4x1y1");
-    g->SetTitle("/HepData/8068/d4x1y1");
-    return g;
-  }
-  TGraphAsymmErrors* CMSsqrts900eta24()
-  {
-    // Plot: p8068_d5x1y1
-    double x[] = { 0.0,        1.0,
-                  2.0, 3.0,
-                  4.0, 5.0,
-                  6.0, 7.0,
-                  8.0, 9.0,
-                  10.0,        11.0,
-                  12.0,        13.0,
-                  14.0,        15.0,
-                  16.0,        17.0,
-                  18.0,        19.0,
-                  20.0,        21.0,
-                  22.0,        23.0,
-                  24.0,        25.0,
-                  26.0,        27.0,
-                  28.0,        29.0,
-                  30.0,        31.0,
-                  32.0,        33.0,
-                  34.0,        35.0,
-                  36.0,        37.0,
-                  38.0,        39.0,
-                  40.0,        41.0,
-                  42.0,        43.0,
-                  44.0,        45.0,
-                  46.0,        47.0,
-                  48.0,        49.0,
-                  50.0,        51.0,
-                  52.0,        53.0,
-                  54.0,        55.0,
-                  56.0,        57.5,
-                  60.0,        66.5,
-                  76.5,        86.5 };
-    double xem[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.0,
-                    1.5,       5.0,
-                    5.0,       5.0 };
-    double xep[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.0,
-                    1.5,       5.0,
-                    5.0,       5.0 };
-    double y[] = { 0.0561,     0.02321,
-                  0.02645,     0.031,
-                  0.0362,      0.04081,
-                  0.04402,     0.04548,
-                  0.04592,     0.04515,
-                  0.04342,     0.04101,
-                  0.03836,     0.03573,
-                  0.03316,     0.03079,
-                  0.02868,     0.02679,
-                  0.02506,     0.02343,
-                  0.02187,     0.02039,
-                  0.01901,     0.0178,
-                  0.01668,     0.01556,
-                  0.01445,     0.01334,
-                  0.01229,     0.01132,
-                  0.01045,     0.00967,
-                  0.00896,     0.00829,
-                  0.00763,     0.00699,
-                  0.00637,     0.00578,
-                  0.00527,     0.00482,
-                  0.00441,     0.00412,
-                  0.00363,     0.00328,
-                  0.00292,     0.00261,
-                  0.00233,     0.00209,
-                  0.00185,     0.00164,
-                  0.00146,     0.00131,
-                  0.00117,     0.00106,
-                  9.7E-4,      8.9E-4,
-                  8.2E-4,      6.63E-4,
-                  4.53E-4,     2.12E-4,
-                  7.7E-5,      1.12E-5 };
-    double yem[] = { 0.01767766952966369,      0.003660737630587584,
-                    0.005079094407470687,      0.006069810540700591,
-                    0.006147235150862541,      0.005458580401532985,
-                    0.004473712105176193,      0.0035970126494078384,
-                    0.003456848275524976,      0.003228761372415125,
-                    0.0029797483115189447,     0.002739361239413305,
-                    0.0025401181074902798,     0.002313114783144148,
-                    0.0020573040611440983,     0.001820109886792553,
-                    0.0016324827717314506,     0.001484318025222358,
-                    0.0013727709204379296,     0.001287866452703851,
-                    0.001217538500417954,      0.0011279184367674819,
-                    0.0010412012293500234,     9.741663102366044E-4,
-                    9.360555539069249E-4,      9.139474820797965E-4,
-                    9.013878188659973E-4,      8.697700845625814E-4,
-                    8.190848551890091E-4,      7.433034373659253E-4,
-                    6.741661516273269E-4,      6.239390995922599E-4,
-                    5.961543424315552E-4,      5.8309518948453E-4,
-                    5.738466694161429E-4,      5.515432893255072E-4,
-                    5.292447448959697E-4,      4.976946855251721E-4,
-                    4.609772228646444E-4,      4.295346318982906E-4,
-                    4.0249223594996216E-4,     4.386342439892262E-4,
-                    3.712142238654117E-4,      3.6674241641784496E-4,
-                    3.6249137920783716E-4,     3.5341194094144583E-4,
-                    3.3105890714493693E-4,     3.130495168499705E-4,
-                    2.8178005607210744E-4,     2.5553864678361276E-4,
-                    2.3323807579381201E-4,     2.1633307652783938E-4,
-                    1.8601075237738274E-4,     1.6401219466856724E-4,
-                    1.562049935181331E-4,      1.562049935181331E-4,
-                    1.5480633061990717E-4,     1.395886814895821E-4,
-                    1.0756393447619885E-4,     6.676076692189808E-5,
-                    2.9966648127543395E-5,     6.551335741663679E-6 };
-    double yep[] = { 0.09323352401362935,      0.004710169848317574,
-                    0.005736898116578331,      0.006618164398078972,
-                    0.006715772479767314,      0.006096630216767292,
-                    0.00519042387479096,       0.004362247127341595,
-                    0.004172217635742412,      0.003913693907295255,
-                    0.003644283194264683,      0.0033638668225719043,
-                    0.0031145144083789367,     0.0028765604460883488,
-                    0.002629087294100369,      0.0023903974564912843,
-                    0.002181604913819182,      0.0020026232796010335,
-                    0.0018514858897652987,     0.0017280335644888382,
-                    0.00160822883943797,       0.0014787156589419076,
-                    0.0013612494260788505,     0.0012539936203984452,
-                    0.0011763077828527702,     0.0011157060544785082,
-                    0.001084158659975559,      0.001033247308247159,
-                    9.630160954002794E-4,      9.244457799135652E-4,
-                    8.35224520712844E-4,       6.989277502002622E-4,
-                    6.519202405202648E-4,      6.203224967708329E-4,
-                    6.10982814815605E-4,       5.88727441181401E-4,
-                    5.664803615307418E-4,      5.348831648126533E-4,
-                    4.976946855251721E-4,      4.477722635447622E-4,
-                    4.204759208325728E-4,      3.9357337308308855E-4,
-                    3.981205847478877E-4,      3.8483762809787716E-4,
-                    3.807886552931954E-4,      3.716180835212409E-4,
-                    3.584689665786984E-4,      3.4928498393145964E-4,
-                    3.1780497164141406E-4,     2.9068883707497264E-4,
-                    2.505992817228334E-4,      2.2472205054244234E-4,
-                    2.0248456731316588E-4,     1.8027756377319944E-4,
-                    1.6401219466856724E-4,     1.6401219466856724E-4,
-                    1.5640012787718557E-4,     1.4223923509355636E-4,
-                    1.0288342918079667E-4,     6.86804193347711E-5,
-                    3.361547262794322E-5,      7.071067811865475E-6 };
-    double ysm[] = { 0.0025,   8.1E-4,
-                    4.4E-4,    4.9E-4,
-                    4.6E-4,    4.5E-4,
-                    4.6E-4,    4.4E-4,
-                    4.3E-4,    4.3E-4,
-                    4.2E-4,    4.0E-4,
-                    3.9E-4,    3.9E-4,
-                    3.9E-4,    3.8E-4,
-                    3.7E-4,    3.6E-4,
-                    3.4E-4,    3.1E-4,
-                    3.0E-4,    2.9E-4,
-                    2.9E-4,    2.9E-4,
-                    2.9E-4,    2.8E-4,
-                    2.7E-4,    2.6E-4,
-                    2.5E-4,    2.5E-4,
-                    2.4E-4,    2.3E-4,
-                    2.3E-4,    2.2E-4,
-                    2.2E-4,    2.1E-4,
-                    2.0E-4,    1.9E-4,
-                    1.9E-4,    1.8E-4,
-                    1.8E-4,    1.8E-4,
-                    1.7E-4,    1.6E-4,
-                    1.5E-4,    1.5E-4,
-                    1.4E-4,    1.4E-4,
-                    1.3E-4,    1.3E-4,
-                    1.2E-4,    1.2E-4,
-                    1.1E-4,    1.0E-4,
-                    1.0E-4,    1.0E-4,
-                    9.4E-5,    6.6E-5,
-                    3.7E-5,    1.9E-5,
-                    1.3E-5,    3.4E-6 };
-    double ysp[] = { 0.0025,   8.1E-4,
-                    4.4E-4,    4.9E-4,
-                    4.6E-4,    4.5E-4,
-                    4.6E-4,    4.4E-4,
-                    4.3E-4,    4.3E-4,
-                    4.2E-4,    4.0E-4,
-                    3.9E-4,    3.9E-4,
-                    3.9E-4,    3.8E-4,
-                    3.7E-4,    3.6E-4,
-                    3.4E-4,    3.1E-4,
-                    3.0E-4,    2.9E-4,
-                    2.9E-4,    2.9E-4,
-                    2.9E-4,    2.8E-4,
-                    2.7E-4,    2.6E-4,
-                    2.5E-4,    2.5E-4,
-                    2.4E-4,    2.3E-4,
-                    2.3E-4,    2.2E-4,
-                    2.2E-4,    2.1E-4,
-                    2.0E-4,    1.9E-4,
-                    1.9E-4,    1.8E-4,
-                    1.8E-4,    1.8E-4,
-                    1.7E-4,    1.6E-4,
-                    1.5E-4,    1.5E-4,
-                    1.4E-4,    1.4E-4,
-                    1.3E-4,    1.3E-4,
-                    1.2E-4,    1.2E-4,
-                    1.1E-4,    1.0E-4,
-                    1.0E-4,    1.0E-4,
-                    9.4E-5,    6.6E-5,
-                    3.7E-5,    1.9E-5,
-                    1.3E-5,    3.4E-6 };
-    int np = 62;
-    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
-                                                x,
-                                                y,
-                                                xem,
-                                                xep,
-                                                yem,
-                                                yep);
-    g->SetName("/HepData/8068/d5x1y1");
-    g->SetTitle("/HepData/8068/d5x1y1");
-    return g;
-  }
-  TGraphAsymmErrors* CMSsqrts900eta20()
-  {
-    // Plot: p8068_d6x1y1
-    double x[] = { 0.0,        1.0,
-                  2.0, 3.0,
-                  4.0, 5.0,
-                  6.0, 7.0,
-                  8.0, 9.0,
-                  10.0,        11.0,
-                  12.0,        13.0,
-                  14.0,        15.0,
-                  16.0,        17.0,
-                  18.0,        19.0,
-                  20.0,        21.0,
-                  22.0,        23.0,
-                  24.0,        25.0,
-                  26.0,        27.0,
-                  28.0,        29.0,
-                  30.0,        31.0,
-                  32.0,        33.0,
-                  34.0,        35.0,
-                  36.0,        37.0,
-                  38.0,        39.0,
-                  40.0,        41.0,
-                  42.0,        43.0,
-                  44.0,        45.0,
-                  46.0,        47.0,
-                  48.0,        49.0,
-                  50.0,        51.0,
-                  52.0,        53.0,
-                  54.0,        55.0,
-                  56.0,        57.0,
-                  58.0,        59.0,
-                  60.0,        61.0,
-                  62.5,        64.5,
-                  67.5,        74.5,
-                  84.5,        94.5 };
-    double xem[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    1.0,       1.0,
-                    2.0,       5.0,
-                    5.0,       5.0 };
-    double xep[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    1.0,       1.0,
-                    2.0,       5.0,
-                    5.0,       5.0 };
-    double y[] = { 0.0494,     0.01789,
-                  0.01895,     0.02174,
-                  0.0255,      0.02961,
-                  0.0335,      0.03646,
-                  0.03826,     0.03906,
-                  0.0392,      0.03868,
-                  0.03753,     0.03596,
-                  0.0341,      0.03206,
-                  0.03,        0.02802,
-                  0.02616,     0.02448,
-                  0.023,       0.0217,
-                  0.02058,     0.01959,
-                  0.01866,     0.01776,
-                  0.01685,     0.01595,
-                  0.01504,     0.01415,
-                  0.01326,     0.01241,
-                  0.01161,     0.01087,
-                  0.01015,     0.00948,
-                  0.00887,     0.00829,
-                  0.00777,     0.00734,
-                  0.00681,     0.00624,
-                  0.00578,     0.00535,
-                  0.00494,     0.00457,
-                  0.00422,     0.00388,
-                  0.00356,     0.00326,
-                  0.00298,     0.00274,
-                  0.00253,     0.00234,
-                  0.00212,     0.00189,
-                  0.00167,     0.0015,
-                  0.00136,     0.00123,
-                  0.00112,     0.001041,
-                  9.17E-4,     7.64E-4,
-                  4.92E-4,     2.45E-4,
-                  9.9E-5,      1.64E-5 };
-    double yem[] = { 0.031045933711196384,     0.0025678006153126453,
-                    0.003808017857100988,      0.004928792955683978,
-                    0.005475372133471842,      0.005444713031923721,
-                    0.004928792955683978,      0.0041622710147226115,
-                    0.003405994715204356,      0.0030985964564621835,
-                    0.0029811742652854096,     0.002811547616527239,
-                    0.0026404734423962684,     0.002492308167141455,
-                    0.0023525518060183073,     0.0021654098919142305,
-                    0.001967053634245899,      0.0017886866690396059,
-                    0.0016107451691686056,     0.0014552319402761885,
-                    0.001326989073052224,      0.001225275479229059,
-                    0.001164817582284883,      0.0011168258592994702,
-                    0.0010756393447619887,     0.001027861858422619,
-                    9.7082439194738E-4,        9.202716990106781E-4,
-                    8.825531145489205E-4,      8.509406559801923E-4,
-                    8.228000972289684E-4,      8.099382692526635E-4,
-                    7.912016177940992E-4,      7.725283166331187E-4,
-                    7.316419889536139E-4,      6.815423684555495E-4,
-                    6.407807737440318E-4,      5.950630218724736E-4,
-                    5.636488268416782E-4,      5.636488268416782E-4,
-                    5.322593352868506E-4,      4.919349550499538E-4,
-                    4.7853944456021594E-4,     4.7423622805517506E-4,
-                    4.651881339845203E-4,      4.428317965096906E-4,
-                    4.248529157249601E-4,      3.9824615503479755E-4,
-                    3.671511950137164E-4,      3.3615472627943223E-4,
-                    3.138470965295043E-4,      3.138470965295043E-4,
-                    3.176476034853718E-4,      3.3541019662496844E-4,
-                    3.3105890714493693E-4,     2.996664812754339E-4,
-                    2.641968962724581E-4,      2.3853720883753126E-4,
-                    2.1633307652783938E-4,     1.9416487838947602E-4,
-                    1.6401219466856724E-4,     1.5181897114655994E-4,
-                    1.5966527487215246E-4,     1.3961733416735901E-4,
-                    9.740636529508736E-5,      6.800735254367721E-5,
-                    3.733630940518893E-5,      9.437160589923222E-6 };
-    double yep[] = { 0.027627884464793896,     0.0038148263394288343,
-                    0.004395600072800072,      0.005347317084295637,
-                    0.005874325493194942,      0.0058935897380119695,
-                    0.005446999173857106,      0.004729587719875803,
-                    0.004031935019317648,      0.0037237615390892046,
-                    0.0035759474269066098,     0.00339607125955861,
-                    0.0032050585018061684,     0.003026549190084311,
-                    0.0028468403537957655,     0.00266865134478073,
-                    0.002469412885687608,      0.0022703523955544874,
-                    0.0020813697413001857,     0.0019043371550227129,
-                    0.00175524927004685,       0.001643471934655411,
-                    0.0015533190271158081,     0.001475127113166862,
-                    0.0014046351839534703,     0.0013267252918370102,
-                    0.0012490796611905903,     0.0011691449867317568,
-                    0.0011208925015361642,     0.0011472575996697516,
-                    9.932774033471214E-4,      9.52102935611481E-4,
-                    9.139474820797965E-4,      8.75956619930462E-4,
-                    8.34865258589672E-4,       7.749193506423748E-4,
-                    7.244998274671983E-4,      6.685057965343307E-4,
-                    6.092618484691126E-4,      5.727128425310541E-4,
-                    5.412947441089743E-4,      5.280151512977634E-4,
-                    5.147815070493501E-4,      5.015974481593781E-4,
-                    4.924428900898052E-4,      4.792702786528704E-4,
-                    4.7010637094172637E-4,     4.428317965096906E-4,
-                    4.204759208325728E-4,      3.8910152916687437E-4,
-                    3.6674241641784496E-4,     3.488552708502481E-4,
-                    3.2649655434629013E-4,     3.0886890422961E-4,
-                    3.04138126514911E-4,       2.996664812754339E-4,
-                    2.9068883707497264E-4,     2.641968962724581E-4,
-                    2.418677324489565E-4,      2.1095023109728988E-4,
-                    1.972308292331602E-4,      1.739482681718907E-4,
-                    1.4475496537252186E-4,     1.2083045973594572E-4,
-                    1.097679370308106E-4,      6.992138442565336E-5,
-                    3.54682957019364E-5,       9.078546139112804E-6 };
-    double ysm[] = { 0.0063,   8.0E-4,
-                    3.7E-4,    4.3E-4,
-                    4.1E-4,    4.0E-4,
-                    4.3E-4,    4.3E-4,
-                    4.2E-4,    4.2E-4,
-                    4.3E-4,    4.2E-4,
-                    4.0E-4,    4.0E-4,
-                    3.9E-4,    3.9E-4,
-                    3.8E-4,    3.7E-4,
-                    3.6E-4,    3.6E-4,
-                    3.5E-4,    3.3E-4,
-                    3.2E-4,    3.2E-4,
-                    3.1E-4,    3.1E-4,
-                    3.1E-4,    3.0E-4,
-                    3.0E-4,    2.9E-4,
-                    2.9E-4,    2.8E-4,
-                    2.8E-4,    2.8E-4,
-                    2.7E-4,    2.6E-4,
-                    2.5E-4,    2.5E-4,
-                    2.4E-4,    2.4E-4,
-                    2.3E-4,    2.2E-4,
-                    2.1E-4,    2.0E-4,
-                    2.0E-4,    1.9E-4,
-                    1.9E-4,    1.9E-4,
-                    1.8E-4,    1.7E-4,
-                    1.6E-4,    1.6E-4,
-                    1.5E-4,    1.5E-4,
-                    1.4E-4,    1.3E-4,
-                    1.3E-4,    1.3E-4,
-                    1.2E-4,    1.1E-4,
-                    1.0E-4,    9.3E-5,
-                    7.3E-5,    5.8E-5,
-                    3.2E-5,    2.0E-5,
-                    1.3E-5,    4.1E-6 };
-    double ysp[] = { 0.0063,   8.0E-4,
-                    3.7E-4,    4.3E-4,
-                    4.1E-4,    4.0E-4,
-                    4.3E-4,    4.3E-4,
-                    4.2E-4,    4.2E-4,
-                    4.3E-4,    4.2E-4,
-                    4.0E-4,    4.0E-4,
-                    3.9E-4,    3.9E-4,
-                    3.8E-4,    3.7E-4,
-                    3.6E-4,    3.6E-4,
-                    3.5E-4,    3.3E-4,
-                    3.2E-4,    3.2E-4,
-                    3.1E-4,    3.1E-4,
-                    3.1E-4,    3.0E-4,
-                    3.0E-4,    2.9E-4,
-                    2.9E-4,    2.8E-4,
-                    2.8E-4,    2.8E-4,
-                    2.7E-4,    2.6E-4,
-                    2.5E-4,    2.5E-4,
-                    2.4E-4,    2.4E-4,
-                    2.3E-4,    2.2E-4,
-                    2.1E-4,    2.0E-4,
-                    2.0E-4,    1.9E-4,
-                    1.9E-4,    1.9E-4,
-                    1.8E-4,    1.7E-4,
-                    1.6E-4,    1.6E-4,
-                    1.5E-4,    1.5E-4,
-                    1.4E-4,    1.3E-4,
-                    1.3E-4,    1.3E-4,
-                    1.2E-4,    1.1E-4,
-                    1.0E-4,    9.3E-5,
-                    7.3E-5,    5.8E-5,
-                    3.2E-5,    2.0E-5,
-                    1.3E-5,    4.1E-6 };
-    int np = 68;
-    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
-                                                x,
-                                                y,
-                                                xem,
-                                                xep,
-                                                yem,
-                                                yep);
-    g->SetName("/HepData/8068/d6x1y1");
-    g->SetTitle("/HepData/8068/d6x1y1");
-    return   g;
+    TH1*    hist = static_cast<TH1*>(first);
+    TList*  list = hist->GetListOfFunctions();
+    TObject* o   = list->FindObject("legend");
+    if (o) return static_cast<TLegend*>(o);
+    
+    TLegend* l = new TLegend(0.65, 0.65, 0.9, 0.9, "", "NDC");
+    l->SetName("legend");
+    l->SetBorderSize(0);
+    l->SetFillColor(0);
+    l->SetFillStyle(0);
+    l->SetTextFont(42);
+    list->Add(l);
+    
+    return l;
   }
-  TGraphAsymmErrors* GetALICE(Double_t eta, UShort_t sNN)
+  
+  /* =================================================================
+   *
+   * Measurements from other sources, such as published ALICE, CMS, ...
+   */
+  TGraphAsymmErrors* GetOther(UShort_t type, Double_t eta, UShort_t sNN,
+                             Int_t factor)
   {
-    Double_t aEta = TMath::Abs(eta);
-    TGraphAsymmErrors* g = 0;
-    switch (sNN) { 
-    case 900:
-      if      (aEta <= 0.51) g = ALICEsqrts900eta05();
-      else if (aEta <= 1.01) g = ALICEsqrts900eta10();
-      else if (aEta <= 1.31) g = ALICEsqrts900eta13();
-      break;
+    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);
+      return 0;
     }
-    if (g) { 
-      g->SetMarkerColor(kALICEColor);
+
+
+    for (Int_t j = 0; j < g->GetN(); j++) { 
+      g->SetPoint(j, g->GetX()[j], g->GetY()[j]*factor);
+      g->SetPointError(j, g->GetEXlow()[j], g->GetEXhigh()[j], 
+                      g->GetEYlow()[j]*factor, g->GetEYhigh()[j]*factor);
     }
     return g;
-  }
-  TGraphAsymmErrors* ALICEsqrts900eta05()
-  {
-    // Plot: p7742_d8x1y1
-    double x[] = { 0.0,        1.0,
-                  2.0, 3.0,
-                  4.0, 5.0,
-                  6.0, 7.0,
-                  8.0, 9.0,
-                  10.0,        11.0,
-                  12.0,        13.0,
-                  14.0,        15.0,
-                  16.0,        17.0,
-                  18.0,        19.0,
-                  20.0,        21.5,
-                  23.5,        25.5 };
-    double xem[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.0,
-                    1.0,       1.0 };
-    double xep[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.0,
-                    1.0,       1.0 };
-    double y[] = { 0.179767,   0.155985,
-                  0.140806,    0.120933,
-                  0.097273,    0.075258,
-                  0.057395,    0.043794,
-                  0.033212,    0.024939,
-                  0.018642,    0.013927,
-                  0.010337,    0.007693,
-                  0.005665,    0.004134,
-                  0.002994,    0.002144,
-                  0.001527,    0.001093,
-                  8.03E-4,     4.82E-4,
-                  2.56E-4,     1.02E-4 };
-    double yem[] = { 0.028627623478032542,     0.006826681770816624,
-                    0.004736912179891031,      0.005286866936097409,
-                    0.004605827830043151,      0.003753897441326814,
-                    0.003172958871463669,      0.0023132784095305087,
-                    0.0017387768689512751,     0.0012885437516824954,
-                    9.891637882575362E-4,      7.64476945368531E-4,
-                    5.769662035162892E-4,      4.500277769204919E-4,
-                    3.7107950630558943E-4,     3.1458703088334716E-4,
-                    2.7004073766748605E-4,     2.210701246211256E-4,
-                    1.7819090885900997E-4,     1.3576450198781714E-4,
-                    1.2447088012864694E-4,     9.762171889492625E-5,
-                    6.095900261651268E-5,      3.4828149534536E-5 };
-    double yep[] = { 0.028627623478032542,     0.006826681770816624,
-                    0.004736912179891031,      0.005286866936097409,
-                    0.004605827830043151,      0.003753897441326814,
-                    0.003172958871463669,      0.0023132784095305087,
-                    0.0017387768689512751,     0.0012885437516824954,
-                    9.891637882575362E-4,      7.64476945368531E-4,
-                    5.769662035162892E-4,      4.500277769204919E-4,
-                    3.7107950630558943E-4,     3.1458703088334716E-4,
-                    2.7004073766748605E-4,     2.210701246211256E-4,
-                    1.7819090885900997E-4,     1.3576450198781714E-4,
-                    1.2447088012864694E-4,     9.762171889492625E-5,
-                    6.095900261651268E-5,      3.4828149534536E-5 };
-    double ysm[] = { 0.003595, 0.00312,
-                    0.002816,  0.002419,
-                    0.001945,  0.001505,
-                    0.001148,  8.76E-4,
-                    6.64E-4,   4.99E-4,
-                    3.97E-4,   3.57E-4,
-                    3.03E-4,   2.66E-4,
-                    2.34E-4,   2.01E-4,
-                    1.71E-4,   1.46E-4,
-                    1.26E-4,   9.6E-5,
-                    9.7E-5,    7.1E-5,
-                    4.6E-5,    2.7E-5 };
-    double ysp[] = { 0.003595, 0.00312,
-                    0.002816,  0.002419,
-                    0.001945,  0.001505,
-                    0.001148,  8.76E-4,
-                    6.64E-4,   4.99E-4,
-                    3.97E-4,   3.57E-4,
-                    3.03E-4,   2.66E-4,
-                    2.34E-4,   2.01E-4,
-                    1.71E-4,   1.46E-4,
-                    1.26E-4,   9.6E-5,
-                    9.7E-5,    7.1E-5,
-                    4.6E-5,    2.7E-5 };
-    int np = 24;
-    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
-                                                x,
-                                                y,
-                                                xem,
-                                                xep,
-                                                yem,
-                                                yep);
-    g->SetName("/HepData/7742/d8x1y1");
-    g->SetTitle("/HepData/7742/d8x1y1");
-    return g;
-  }
-  TGraphAsymmErrors* ALICEsqrts900eta10()
-  {
-    // Plot: p7742_d9x1y1
-    double x[] = { 0.0,        1.0,
-                  2.0, 3.0,
-                  4.0, 5.0,
-                  6.0, 7.0,
-                  8.0, 9.0,
-                  10.0,        11.0,
-                  12.0,        13.0,
-                  14.0,        15.0,
-                  16.0,        17.0,
-                  18.0,        19.0,
-                  20.0,        21.0,
-                  22.0,        23.0,
-                  24.0,        25.0,
-                  26.0,        27.0,
-                  28.0,        29.0,
-                  30.0,        31.0,
-                  32.0,        33.0,
-                  34.0,        35.0,
-                  36.0,        37.0,
-                  38.0,        39.0,
-                  40.0,        41.5 };
-    double xem[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.0 };
-    double xep[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.0 };
-    double y[] = { 0.092396,   0.076087,
-                  0.073233,    0.078506,
-                  0.080885,    0.077084,
-                  0.070058,    0.062258,
-                  0.054437,    0.047195,
-                  0.040978,    0.035352,
-                  0.03048,     0.026336,
-                  0.022788,    0.019725,
-                  0.017071,    0.014674,
-                  0.012509,    0.010593,
-                  0.008935,    0.007527,
-                  0.006361,    0.005416,
-                  0.004665,    0.004041,
-                  0.003519,    0.00306,
-                  0.002637,    0.002234,
-                  0.001856,    0.001512,
-                  0.001213,    9.62E-4,
-                  7.61E-4,     6.02E-4,
-                  4.84E-4,     3.93E-4,
-                  3.27E-4,     2.76E-4,
-                  2.37E-4,     1.7E-4 };
-    double yem[] = { 0.021888152137629163,     0.0066076761421849355,
-                    0.004451948000594796,      0.0034909631908686747,
-                    0.002952765483407038,      0.0028889466938661224,
-                    0.0029957812002881653,     0.0032637882897026274,
-                    0.002922475833946279,      0.0026227742945209753,
-                    0.0023487871338203465,     0.002091141554271255,
-                    0.0017347841364273539,     0.0014762689456870655,
-                    0.0012594490859101847,     0.0010506721658062519,
-                    9.019456746389995E-4,      7.836938177630343E-4,
-                    6.432239112470866E-4,      5.679656679765071E-4,
-                    4.967293025381128E-4,      4.475187146924696E-4,
-                    4.040519768544636E-4,      3.3646693745448454E-4,
-                    2.9206163733020465E-4,     2.820017730440715E-4,
-                    2.6982401672201087E-4,     2.5470178640912594E-4,
-                    2.3519353732617736E-4,     2.1914607000811127E-4,
-                    1.9727392123643713E-4,     1.7402586014727813E-4,
-                    1.5060212481900778E-4,     1.3009611831257687E-4,
-                    1.1269871339105873E-4,     9.693812459502196E-5,
-                    8.48528137423857E-5,       7.427651041883969E-5,
-                    6.726812023536856E-5,      6.030754513325841E-5,
-                    5.544366510251645E-5,      4.601086828130936E-5 };
-    double yep[] = { 0.021888152137629163,     0.0066076761421849355,
-                    0.004451948000594796,      0.0034909631908686747,
-                    0.002952765483407038,      0.0028889466938661224,
-                    0.0029957812002881653,     0.0032637882897026274,
-                    0.002922475833946279,      0.0026227742945209753,
-                    0.0023487871338203465,     0.002091141554271255,
-                    0.0017347841364273539,     0.0014762689456870655,
-                    0.0012594490859101847,     0.0010506721658062519,
-                    9.019456746389995E-4,      7.836938177630343E-4,
-                    6.432239112470866E-4,      5.679656679765071E-4,
-                    4.967293025381128E-4,      4.475187146924696E-4,
-                    4.040519768544636E-4,      3.3646693745448454E-4,
-                    2.9206163733020465E-4,     2.820017730440715E-4,
-                    2.6982401672201087E-4,     2.5470178640912594E-4,
-                    2.3519353732617736E-4,     2.1914607000811127E-4,
-                    1.9727392123643713E-4,     1.7402586014727813E-4,
-                    1.5060212481900778E-4,     1.3009611831257687E-4,
-                    1.1269871339105873E-4,     9.693812459502196E-5,
-                    8.48528137423857E-5,       7.427651041883969E-5,
-                    6.726812023536856E-5,      6.030754513325841E-5,
-                    5.544366510251645E-5,      4.601086828130936E-5 };
-    double ysm[] = { 0.001848, 0.001522,
-                    0.001465,  0.00157,
-                    0.001618,  0.001542,
-                    0.001401,  0.001245,
-                    0.001089,  9.44E-4,
-                    8.2E-4,    7.07E-4,
-                    6.1E-4,    5.27E-4,
-                    4.56E-4,   3.94E-4,
-                    3.41E-4,   3.0E-4,
-                    2.76E-4,   2.56E-4,
-                    2.38E-4,   2.23E-4,
-                    2.07E-4,   1.91E-4,
-                    1.8E-4,    1.7E-4,
-                    1.59E-4,   1.47E-4,
-                    1.3E-4,    1.25E-4,
-                    1.14E-4,   1.02E-4,
-                    9.1E-5,    8.2E-5,
-                    7.4E-5,    6.6E-5,
-                    6.0E-5,    5.4E-5,
-                    5.0E-5,    4.6E-5,
-                    4.3E-5,    3.4E-5 };
-    double ysp[] = { 0.001848, 0.001522,
-                    0.001465,  0.00157,
-                    0.001618,  0.001542,
-                    0.001401,  0.001245,
-                    0.001089,  9.44E-4,
-                    8.2E-4,    7.07E-4,
-                    6.1E-4,    5.27E-4,
-                    4.56E-4,   3.94E-4,
-                    3.41E-4,   3.0E-4,
-                    2.76E-4,   2.56E-4,
-                    2.38E-4,   2.23E-4,
-                    2.07E-4,   1.91E-4,
-                    1.8E-4,    1.7E-4,
-                    1.59E-4,   1.47E-4,
-                    1.3E-4,    1.25E-4,
-                    1.14E-4,   1.02E-4,
-                    9.1E-5,    8.2E-5,
-                    7.4E-5,    6.6E-5,
-                    6.0E-5,    5.4E-5,
-                    5.0E-5,    4.6E-5,
-                    4.3E-5,    3.4E-5 };
-    int np = 42;
-    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
-                                                x,
-                                                y,
-                                                xem,
-                                                xep,
-                                                yem,
-                                                yep);
-    g->SetName("/HepData/7742/d9x1y1");
-    g->SetTitle("/HepData/7742/d9x1y1");
-    return g;
-  }
-  TGraphAsymmErrors* ALICEsqrts900eta13()
-  {
-    // Plot: p7742_d10x1y1
-    double x[] = { 0.0,        1.0,
-                  2.0, 3.0,
-                  4.0, 5.0,
-                  6.0, 7.0,
-                  8.0, 9.0,
-                  10.0,        11.0,
-                  12.0,        13.0,
-                  14.0,        15.0,
-                  16.0,        17.0,
-                  18.0,        19.0,
-                  20.0,        21.0,
-                  22.0,        23.0,
-                  24.0,        25.0,
-                  26.0,        27.0,
-                  28.0,        29.0,
-                  30.0,        31.0,
-                  32.0,        33.0,
-                  34.0,        35.0,
-                  36.0,        37.0,
-                  38.0,        39.0,
-                  40.0,        41.5,
-                  43.5,        45.5,
-                  47.5,        49.5,
-                  51.5,        53.5 };
-    double xem[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.0,
-                    1.0,       1.0,
-                    1.0,       1.0,
-                    1.0,       1.0 };
-    double xep[] = { 0.5,      0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       0.5,
-                    0.5,       1.0,
-                    1.0,       1.0,
-                    1.0,       1.0,
-                    1.0,       1.0 };
-    double y[] = { 0.068695,   0.056382,
-                  0.049828,    0.056116,
-                  0.063345,    0.066198,
-                  0.065598,    0.061117,
-                  0.055042,    0.04927,
-                  0.044132,    0.039615,
-                  0.035906,    0.03262,
-                  0.029489,    0.026657,
-                  0.023839,    0.021302,
-                  0.018884,    0.016732,
-                  0.014707,    0.012944,
-                  0.01137,     0.009982,
-                  0.008753,    0.007661,
-                  0.006711,    0.005891,
-                  0.005179,    0.004572,
-                  0.004064,    0.003636,
-                  0.003278,    0.002968,
-                  0.002677,    0.002392,
-                  0.002111,    0.001829,
-                  0.001561,    0.001319,
-                  0.001115,    7.93E-4,
-                  5.24E-4,     3.51E-4,
-                  2.2E-4,      1.56E-4,
-                  1.27E-4,     8.5E-5 };
-    double yem[] = { 0.021322725646595934,     0.007630831737104416,
-                    0.004038016592338372,      0.0036354708085748677,
-                    0.0035003511252444373,     0.0030936373413831173,
-                    0.0026374775828431228,     0.002920906195001818,
-                    0.00276560481631053,       0.002864668392676542,
-                    0.0026843833556330957,     0.0022856071403458645,
-                    0.002325611317481922,      0.002006865217198205,
-                    0.0017203360136903488,     0.0014989629748596194,
-                    0.0012394615766533467,     0.0011637718848640398,
-                    9.151568171630477E-4,      9.078199160626517E-4,
-                    7.98723982361867E-4,       6.967302203866285E-4,
-                    6.308502199413107E-4,      5.846990678973245E-4,
-                    5.4516236113657E-4,        5.460009157501477E-4,
-                    4.930608481719068E-4,      4.3796118549478783E-4,
-                    3.4394330928221294E-4,     3.2323366161339073E-4,
-                    3.0869078379504626E-4,     2.9443165590676553E-4,
-                    2.756320010448714E-4,      2.630304164920856E-4,
-                    2.5628889948649746E-4,     2.4142286552851616E-4,
-                    2.245551157288562E-4,      2.0678491240900531E-4,
-                    1.8828170383762728E-4,     1.7255723688098393E-4,
-                    1.5911316727411343E-4,     1.3023056476879765E-4,
-                    1.0761505470890214E-4,     8.62670273047588E-5,
-                    6.129437168288782E-5,      6.037383539249432E-5,
-                    3.535533905932738E-5,      3.8078865529319545E-5 };
-    double yep[] = { 0.021322725646595934,     0.007630831737104416,
-                    0.004038016592338372,      0.0036354708085748677,
-                    0.0035003511252444373,     0.0030936373413831173,
-                    0.0026374775828431228,     0.002920906195001818,
-                    0.00276560481631053,       0.002864668392676542,
-                    0.0026843833556330957,     0.0022856071403458645,
-                    0.002325611317481922,      0.002006865217198205,
-                    0.0017203360136903488,     0.0014989629748596194,
-                    0.0012394615766533467,     0.0011637718848640398,
-                    9.151568171630477E-4,      9.078199160626517E-4,
-                    7.98723982361867E-4,       6.967302203866285E-4,
-                    6.308502199413107E-4,      5.846990678973245E-4,
-                    5.4516236113657E-4,        5.460009157501477E-4,
-                    4.930608481719068E-4,      4.3796118549478783E-4,
-                    3.4394330928221294E-4,     3.2323366161339073E-4,
-                    3.0869078379504626E-4,     2.9443165590676553E-4,
-                    2.756320010448714E-4,      2.630304164920856E-4,
-                    2.5628889948649746E-4,     2.4142286552851616E-4,
-                    2.245551157288562E-4,      2.0678491240900531E-4,
-                    1.8828170383762728E-4,     1.7255723688098393E-4,
-                    1.5911316727411343E-4,     1.3023056476879765E-4,
-                    1.0761505470890214E-4,     8.62670273047588E-5,
-                    6.129437168288782E-5,      6.037383539249432E-5,
-                    3.535533905932738E-5,      3.8078865529319545E-5 };
-    double ysm[] = { 0.012402, 0.001128,
-                    9.97E-4,   0.001122,
-                    0.001267,  0.001324,
-                    0.001312,  0.001222,
-                    0.001101,  9.85E-4,
-                    8.83E-4,   7.92E-4,
-                    7.18E-4,   6.52E-4,
-                    5.9E-4,    5.33E-4,
-                    4.77E-4,   4.26E-4,
-                    3.94E-4,   3.81E-4,
-                    3.58E-4,   3.33E-4,
-                    3.16E-4,   3.08E-4,
-                    2.91E-4,   2.71E-4,
-                    2.55E-4,   2.39E-4,
-                    2.24E-4,   2.12E-4,
-                    2.07E-4,   1.99E-4,
-                    1.82E-4,   1.72E-4,
-                    1.72E-4,   1.62E-4,
-                    1.49E-4,   1.38E-4,
-                    1.27E-4,   1.2E-4,
-                    1.14E-4,   9.6E-5,
-                    8.5E-5,    7.1E-5,
-                    5.1E-5,    5.4E-5,
-                    2.5E-5,    3.3E-5 };
-    double ysp[] = { 0.012402, 0.001128,
-                    9.97E-4,   0.001122,
-                    0.001267,  0.001324,
-                    0.001312,  0.001222,
-                    0.001101,  9.85E-4,
-                    8.83E-4,   7.92E-4,
-                    7.18E-4,   6.52E-4,
-                    5.9E-4,    5.33E-4,
-                    4.77E-4,   4.26E-4,
-                    3.94E-4,   3.81E-4,
-                    3.58E-4,   3.33E-4,
-                    3.16E-4,   3.08E-4,
-                    2.91E-4,   2.71E-4,
-                    2.55E-4,   2.39E-4,
-                    2.24E-4,   2.12E-4,
-                    2.07E-4,   1.99E-4,
-                    1.82E-4,   1.72E-4,
-                    1.72E-4,   1.62E-4,
-                    1.49E-4,   1.38E-4,
-                    1.27E-4,   1.2E-4,
-                    1.14E-4,   9.6E-5,
-                    8.5E-5,    7.1E-5,
-                    5.1E-5,    5.4E-5,
-                    2.5E-5,    3.3E-5 };
-    int np = 48;
-    TGraphAsymmErrors* g = new TGraphAsymmErrors(np,
-                                                x,
-                                                y,
-                                                xem,
-                                                xep,
-                                                yem,
-                                                yep);
-    g->SetName("/HepData/7742/d10x1y1");
-    g->SetTitle("/HepData/7742/d10x1y1");
-    return g;
-  }
-  ClassDef(Unfolder,1);
+  }    
 };
-void UnfoldMultDists(const char* method="Bayes", 
-                    Double_t    regParam=4,
-                    const char* realFile="forward_multdists.root",
-                     const char* mcFile="forward_mcmultdists.root")
+
+void
+UnfoldMultDists(const TString& method="Bayes", 
+               Double_t       regParam=1e-30,
+               const TString& measured="forward_multdists.root", 
+               const TString& corrections="")
 {
-  Unfolder* u = new Unfolder;
-  u->Run(method, regParam, realFile, mcFile);
-  
-  u->DrawAll();
+  Unfolder m;
+  m.Run(measured, corrections, method, regParam);
 }
-#endif