A better way to specify the Nch axis for the MultDists task.
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / scripts / UnfoldMultDists.C
index 4ceb8853bf218d85f035fea8d08d6c842e3dc0ed..168284661e5711800f5afc6f2fab861c8fce5a75 100644 (file)
@@ -280,13 +280,11 @@ struct Unfolder
     // Get some info from the input collection 
     UShort_t sys;
     UShort_t sNN; 
-    UShort_t maxN; 
     ULong_t  trig; 
     Double_t minZ; 
     Double_t maxZ;
     GetParameter(mTop, "sys",     sys);          
     GetParameter(mTop, "snn",     sNN);          
-    GetParameter(mTop, "maxN",    maxN);         
     GetParameter(mTop, "trigger", trig);
     GetParameter(mTop, "minIpZ",  minZ); 
     GetParameter(mTop, "maxIpZ",  maxZ); 
@@ -307,7 +305,8 @@ struct Unfolder
     if (mId == 0xDeadBeef) return;
 
     // Store information 
-    SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,maxN);
+    SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,
+                   corrFile.IsNull());
 
     // Load other data 
     TString savPath(gROOT->GetMacroPath());
@@ -336,6 +335,12 @@ struct Unfolder
 
     SaveSummarize();
   }
+  /** 
+   * Append an & to a string and the next term.
+   * 
+   * @param trg  Output string
+   * @param what Term
+   */
   static void AppendAnd(TString& trg, const TString& what)
   {
     if (!trg.IsNull()) trg.Append(" & ");
@@ -353,7 +358,6 @@ struct Unfolder
    * @param trigger  Trigger mask 
    * @param minIpZ   Least z coordinate of interaction point
    * @param maxIpZ   Largest z coordinate of interaction point
-   * @param maxN     Largest @f$N_{ch}@f$ to consider 
    */
   void SaveInformation(TDirectory* dir, 
                       const TString& method,
@@ -364,15 +368,20 @@ struct Unfolder
                       UInt_t         trigger,
                       Double_t       minIpZ, 
                       Double_t       maxIpZ, 
-                      UShort_t       maxN) const
+                      Bool_t         self) const
   {
     dir->cd();
 
+    TParameter<bool>* pM = new TParameter<bool>("self", self);
+    pM->SetBit(BIT(19));
+    pM->Write();
+
     TNamed* outMeth = new TNamed("method", method.Data());
     outMeth->SetUniqueID(mId);
     outMeth->Write();
 
     TParameter<double>* pR = new TParameter<double>("regParam", regParam);
+    pR->SetBit(BIT(19));
     pR->Write();
     
     TString tS = (sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "?");
@@ -435,13 +444,12 @@ struct Unfolder
     pT->Write();
     
     TParameter<double>* pY = new TParameter<double>("minIpZ", minIpZ);
+    pY->SetBit(BIT(19));
     pY->Write();
 
     TParameter<double>* pZ = new TParameter<double>("maxIpZ", maxIpZ);
+    pZ->SetBit(BIT(19));
     pZ->Write();
-
-    TParameter<int>* pN = new TParameter<int>("maxN", maxN);
-    pN->Write();
   }
   /** 
    * Save a script to do a summary of this step 
@@ -493,6 +501,8 @@ struct Unfolder
                                         "Unfolded P(#it{N}_{ch})");
     THStack*  allCorrected = new THStack("corrected",
                                         "Corrected P(#it{N}_{ch})");
+    THStack*  allRatio     = (sys != 1 ? 0 : 
+                             new THStack("ratios", "Ratios to other"));
     TMultiGraph* allALICE  = (sys != 1 ? 0 : 
                              new TMultiGraph("alice", "ALICE Published"));
     TMultiGraph* allCMS    = (sys != 1 ? 0 : 
@@ -525,9 +535,14 @@ struct Unfolder
       THStack* binS = ProcessBin(mBin, cBin, method, r, dir);
       if (!binS) continue;
 
+      TH1* result = 0;
       Bin2Stack(binS, i, allMeasured, allTruth, allTruthA, 
-               allUnfolded, allCorrected);
-      Other2Stack(o->GetName(), i, sNN, allALICE, allCMS);
+               allUnfolded, allCorrected, result);
+
+      TGraph* alice = 0;
+      TGraph* cms   = 0;
+      Other2Stack(o->GetName(), i, sNN, allALICE, allCMS, alice, cms);
+      Ratio2Stack(i, result, alice, cms, allRatio);
       i++;
     }
     dir->Add(allMeasured);
@@ -547,6 +562,12 @@ struct Unfolder
       else 
        delete allCMS;
     }
+    if (allRatio && allRatio->GetHists()) { 
+      if (allRatio->GetHists()->GetEntries() > 0) 
+       dir->Add(allRatio);
+      else 
+       delete allRatio;
+    }
   }
   /** 
    * Process a single eta bin 
@@ -749,7 +770,8 @@ struct Unfolder
                 THStack* truth, 
                 THStack* accepted, 
                 THStack* unfolded,
-                THStack* corrected)
+                THStack* corrected,
+                TH1*&    result)
   {
     Int_t open, closed;
     Double_t factor; 
@@ -776,6 +798,7 @@ struct Unfolder
       cln->SetMarkerStyle(sty);
       cln->SetMarkerSize(size);
       cln->Scale(factor); // Scale by 10^i
+      if (col == kColorCorrected) result = cln;
 
       // Make sure we do not get the old legend 
       TObject* tst = cln->FindObject("legend");
@@ -816,7 +839,8 @@ struct Unfolder
    * @param allCMS   Stack of CMS data 
    */
   void Other2Stack(const TString& name, Int_t i,
-                  UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS) 
+                  UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS,
+                  TGraph*& alice, TGraph*& cms) 
   {
     if (!allALICE && !allCMS) return;
 
@@ -824,6 +848,7 @@ struct Unfolder
     tmp.ReplaceAll("p", "+");
     tmp.ReplaceAll("m", "-");
     tmp.ReplaceAll("_", " ");
+    tmp.ReplaceAll("d", ".");
     TObjArray* tokens = tmp.Tokenize(" ");
     if (!tokens || tokens->GetEntriesFast() < 2) { 
       Error("Other2Stack", "Failed to decode eta range from %s", name.Data());
@@ -851,6 +876,7 @@ struct Unfolder
        g->SetMarkerColor(kColorALICE);
        g->SetMarkerSize(size);
        allALICE->Add(g, "p same");
+       alice = g;
       }
     }
     if (allCMS) {
@@ -860,11 +886,74 @@ struct Unfolder
        g->SetMarkerColor(kColorCMS);
        g->SetMarkerSize(size);
        allCMS->Add(g, "p same");
+       cms = g;
+      }
+    }
+  }
+  /** 
+   * Create ratios to other data 
+   * 
+   * @param i 
+   * @param res 
+   * @param alice 
+   * @param cms 
+   * @param all 
+   */
+  void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
+  {
+    if (!all || !res || !(alice || cms)) return;
+
+    Int_t        off  = 5*ib;
+    TGraph*      gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
+    TGraph**     pg   = gs;
+    while (*pg) { 
+      TGraph*     g = *pg;
+      const char* n = (g == alice ? "ALICE" : "CMS");
+
+      TH1*    r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
+      TString tit(r->GetTitle());
+      tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
+      r->SetTitle(tit);
+      r->SetMarkerColor(g->GetMarkerColor());
+      r->SetLineColor(g->GetLineColor());
+
+      TObject* tst = r->FindObject("legend");
+      if (tst) r->GetListOfFunctions()->Remove(tst);
+
+      for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
+       Double_t c = r->GetBinContent(i);
+       Double_t e = r->GetBinError(i);
+       Double_t o = g->Eval(r->GetBinCenter(i));
+       if (o < 1e-12) { 
+         r->SetBinContent(i, 0);
+         r->SetBinError(i, 0);
+         continue;
+       }
+       r->SetBinContent(i, (c - o) / o + off);
+       r->SetBinError(i, e / o);
       }
+      all->Add(r);
+      pg++;
     }
+    TLegend* leg = StackLegend(all);
+    if (!leg) return;
+      
+    TString   txt      = res->GetTitle();
+    txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
+    if      (ib == 0) txt.Append(" "); // (#times1)");
+    // else if (ib == 1) txt.Append(" (#times10)");
+    else              txt.Append(Form(" (+%d)", off));
 
-    
+    TObject* dummy = 0;
+    TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
+    e->SetMarkerStyle(res->GetMarkerStyle());
+    e->SetMarkerSize(res->GetMarkerSize());
+    e->SetMarkerColor(kBlack);
+    e->SetFillColor(0);
+    e->SetFillStyle(0);
+    e->SetLineColor(kBlack);
   }
+
   /** 
    * Get or create a stack legend.  This is done by adding a TLegend
    * object to the list of functions for the first histogram in the
@@ -905,16 +994,16 @@ struct Unfolder
   TGraphAsymmErrors* GetOther(UShort_t type, Double_t eta, UShort_t sNN,
                              Int_t factor)
   {
-    TString oC = Form("OtherPNch::GetData(%hu,%f,%hu)", 
+    TString oC = Form("OtherPNch::GetData(%hu,%f,%hu);", 
                       type, eta, sNN);
     TGraphAsymmErrors* g = 
       reinterpret_cast<TGraphAsymmErrors*>(gROOT->ProcessLine(oC));
     if (!g) { 
-      Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
-             type, eta, sNN);
+      // Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
+      //         type, eta, sNN);
       return 0;
     }
-
+  
 
     for (Int_t j = 0; j < g->GetN(); j++) { 
       g->SetPoint(j, g->GetX()[j], g->GetY()[j]*factor);