AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UnfoldMultDists.C
Go to the documentation of this file.
1 
11 #include <TFile.h>
12 #include <TList.h>
13 #include <TH1.h>
14 #include <TH2.h>
15 #include <THStack.h>
16 #include <TLegend.h>
17 #include <TLegendEntry.h>
18 #include <TClass.h>
19 #include <TRegexp.h>
20 #include <TMath.h>
21 #include <TParameter.h>
22 #include <TMultiGraph.h>
23 #include <TGraphAsymmErrors.h>
24 #include "RooUnfold.h"
25 #include "RooUnfoldResponse.h"
26 #include <fstream>
27 
33 struct Unfolder
34 {
39  enum {
40  kColorMeasured = kOrange-2,
41  kColorTruth = kBlue-3,
42  kColorAccepted = kMagenta-3,
43  kColorTrgVtx = kBlack,
44  kColorUnfolded = kOrange+2,
45  kColorCorrected= kRed+2,
46  kColorError = kBlue-10,
47  kColorALICE = kPink+1,
48  kColorCMS = kGreen+2
49  };
50 
54  Unfolder() {}
63  static TCollection* GetTop(const TString& fileName, Bool_t results=false)
64  {
65  TFile* file = TFile::Open(fileName, "READ");
66  if (!file) {
67  Error("GetTop", "Failed to open %s", fileName.Data());
68  return 0;
69  }
70  TCollection* ret = 0;
71  TString cName(Form("ForwardMult%s", results ? "Results" : "Sums"));
72  file->GetObject(cName, ret);
73  if (!ret)
74  Error("GetTop", "Failed to get collection %s from %s",
75  cName.Data(), fileName.Data());
76  file->Close();
77  return ret;
78  }
89  static TObject* GetObject(TCollection* c, const TString& name,
90  TClass* cl, Bool_t verbose=true)
91  {
92  TObject* o = c->FindObject(name);
93  if (!o) {
94  if (verbose)
95  Warning("GetObject", "%s not found in %s", name.Data(), c->GetName());
96  return 0;
97  }
98  if (cl && !o->IsA()->InheritsFrom(cl)) {
99  if (verbose)
100  Warning("GetCollection", "%s is not a %s but a %s",
101  name.Data(), cl->GetName(), o->ClassName());
102  return 0;
103  }
104  return o;
105  }
116  const TString& name,
117  Bool_t verbose=-true)
118  {
119  return static_cast<TCollection*>(GetObject(c, name,
120  TCollection::Class(),
121  verbose));
122  }
132  static TH1* GetH1(TCollection* c, const TString& name, Bool_t verbose=true)
133  {
134  return static_cast<TH1*>(GetObject(c, name, TH1::Class(), verbose));
135  }
145  static TH2* GetH2(TCollection* c, const TString& name, Bool_t verbose=true)
146  {
147  return static_cast<TH2*>(GetObject(c, name, TH2::Class(), verbose));
148  }
158  static void GetParameter(TCollection* c, const TString& name, UShort_t& v)
159  {
160  TObject* o = GetObject(c, name, TParameter<int>::Class(), true);
161  v = (!o ? 0 : o->GetUniqueID());
162  }
172  static void GetParameter(TCollection* c, const TString& name, ULong_t& v)
173  {
174  TObject* o = GetObject(c, name, TParameter<long>::Class(), true);
175  v = (!o ? 0 : o->GetUniqueID());
176  }
186  static void GetParameter(TCollection* c, const TString& name, Double_t& v)
187  {
188  TObject* o = GetObject(c, name, TParameter<double>::Class(), true);
189  v = (!o ? 0 : static_cast<TParameter<double>*>(o)->GetVal());
190  }
199  {
200  struct Method {
201  UInt_t id;
202  TString name;
203  };
204  const Method methods[] = { {RooUnfold::kNone, "None"},
205  {RooUnfold::kBayes, "Bayes"},
206  {RooUnfold::kSVD, "SVD"},
207  {RooUnfold::kBinByBin,"BinByBin"},
208  {RooUnfold::kTUnfold, "TUnfold"},
209  {RooUnfold::kInvert, "Invert"},
210  {RooUnfold::kDagostini,"Dagostini"},
211  {0xDeadBeef, "unknown"} };
212  const Method* pMethod = methods;
213  while (pMethod->id != 0xDeadBeef) {
214  if (method.BeginsWith(pMethod->name, TString::kIgnoreCase)) {
215  method = pMethod->name;
216  break;
217  }
218  pMethod++;
219  }
220  if (pMethod->id == 0xDeadBeef)
221  Error("MethodId", "Unknown unfolding method: %s", method.Data());
222 
223  return pMethod->id;
224  }
225 
285  void Run(const TString& measuredFile, const TString& corrFile,
286  const TString& method="Bayes", Double_t regParam=4)
287  {
288  // Get the input collections
289  if (measuredFile.IsNull()) {
290  Error("Run", "No measurements given");
291  return;
292  }
293  TCollection* mTop = GetTop(measuredFile, false);
294  TCollection* cTop = GetTop(corrFile.IsNull() ? measuredFile : corrFile,
295  true);
296  if (!mTop || !cTop) return;
297 
298  // Get some info from the input collection
299  UShort_t sys;
300  UShort_t sNN;
301  ULong_t trig;
302  Double_t minZ;
303  Double_t maxZ;
304  GetParameter(mTop, "sys", sys);
305  GetParameter(mTop, "sNN", sNN);
306  GetParameter(mTop, "trigger", trig);
307  GetParameter(mTop, "minIpZ", minZ);
308  GetParameter(mTop, "maxIpZ", maxZ);
309  if (sys == 0 || sNN == 0)
310  Warning("Run", "System (%d) and/or collision energy (%d) unknown",
311  sys, sNN);
312 
313  // Open the output file
314  TFile* out = TFile::Open("forward_unfolded.root", "RECREATE");
315  if (!out) {
316  Error("Run", "Failed to open output file");
317  return;
318  }
319 
320  // Decode method option and store in file
321  TString meth(method);
322  UInt_t mId = MethodId(meth);
323  if (mId == 0xDeadBeef) return;
324 
325  // Store information
326  SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,
327  corrFile.IsNull());
328 
329  // Load other data
330  TString savPath(gROOT->GetMacroPath());
331  gROOT->SetMacroPath(Form("%s:$(ALICE_PHYSICS)/PWGLF/FORWARD/analysis2/scripts",
332  gROOT->GetMacroPath()));
333  // Always recompile
334  if (!gROOT->GetClass("OtherPNch"))
335  gROOT->LoadMacro("OtherPNchData.C++");
336  gROOT->SetMacroPath(savPath);
337 
338  // Loop over the input
339  const char* inputs[] = { "symmetric", "positive", "negative", 0 };
340  const char** pinput = inputs;
341  while (*pinput) {
342  TCollection* mInput = GetCollection(mTop, *pinput, false);
343  TCollection* cInput = GetCollection(cTop, *pinput, false);
344  if (mInput && cInput)
345  ProcessType(mInput, cInput, mId, regParam, out,
346  sys, sNN);
347  pinput++;
348  }
349 
350  out->Write();
351  // out->ls();
352  out->Close();
353 
354  SaveSummarize();
355  }
362  static void AppendAnd(TString& trg, const TString& what)
363  {
364  if (!trg.IsNull()) trg.Append(" & ");
365  trg.Append(what);
366  }
381  void SaveInformation(TDirectory* dir,
382  const TString& method,
383  Int_t mId,
384  Double_t regParam,
385  UShort_t sys,
386  UShort_t sNN,
387  UInt_t trigger,
388  Double_t minIpZ,
389  Double_t maxIpZ,
390  Bool_t self) const
391  {
392  dir->cd();
393 
394  TParameter<bool>* pM = new TParameter<bool>("self", self);
395  pM->SetBit(BIT(19));
396  pM->Write();
397 
398  TNamed* outMeth = new TNamed("method", method.Data());
399  outMeth->SetUniqueID(mId);
400  outMeth->Write();
401 
402  TParameter<double>* pR = new TParameter<double>("regParam", regParam);
403  pR->SetBit(BIT(19));
404  pR->Write();
405 
406  TString tS = (sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "?");
407  TNamed* pS = new TNamed("sys", tS.Data()); pS->SetUniqueID(sys);
408  pS->Write();
409 
410  TString tE;
411  if (sNN < 1000) tE = Form("%dGeV", sNN);
412  else if (sNN < 3000) tE = Form("%4.2fTeV", float(sNN)/1000);
413  else tE = Form("%dTeV", sNN/1000);
414  TNamed* pE = new TNamed("sNN", tE.Data()); pE->SetUniqueID(sNN);
415  pE->Write();
416 
417  TString tT;
421  enum {
423  kInel = 0x0001,
425  kInelGt0 = 0x0002,
427  kNSD = 0x0004,
429  kEmpty = 0x0008,
431  kA = 0x0010,
433  kB = 0x0020,
435  kC = 0x0080,
437  kE = 0x0100,
439  kPileUp = 0x0200,
441  kMCNSD = 0x0400,
443  kOffline = 0x0800,
445  kNClusterGt0 = 0x1000,
447  kV0AND = 0x2000,
449  kSatellite = 0x4000
450  };
451  if ((trigger & kInel) != 0x0) AppendAnd(tT, "INEL");
452  if ((trigger & kInelGt0) != 0x0) AppendAnd(tT, "INEL>0");
453  if ((trigger & kNSD) != 0x0) AppendAnd(tT, "NSD");
454  if ((trigger & kV0AND) != 0x0) AppendAnd(tT, "V0AND");
455  if ((trigger & kA) != 0x0) AppendAnd(tT, "A");
456  if ((trigger & kB) != 0x0) AppendAnd(tT, "B");
457  if ((trigger & kC) != 0x0) AppendAnd(tT, "C");
458  if ((trigger & kE) != 0x0) AppendAnd(tT, "E");
459  if ((trigger & kMCNSD) != 0x0) AppendAnd(tT, "MCNSD");
460  if ((trigger & kNClusterGt0) != 0x0) AppendAnd(tT, "NCluster>0");
461  if ((trigger & kSatellite) != 0x0) AppendAnd(tT, "Satellite");
462  TNamed* pT = new TNamed("trigger", tT.Data()); pT->SetUniqueID(trigger);
463  pT->Write();
464 
465  TParameter<double>* pY = new TParameter<double>("minIpZ", minIpZ);
466  pY->SetBit(BIT(19));
467  pY->Write();
468 
469  TParameter<double>* pZ = new TParameter<double>("maxIpZ", maxIpZ);
470  pZ->SetBit(BIT(19));
471  pZ->Write();
472  }
478  {
479  std::ofstream f("SummarizeUnfold.C");
480  f << "void SummarizeUnfold(const char* file=\"forward_unfolded.root\")\n"
481  << "{\n"
482  << " const char* fwd=\"$ALICE_PHYSICS/PWGLF/FORWARD/analysis2\";\n"
483  << " gROOT->LoadMacro(Form(\"%s/DrawUnfoldedSummary.C\",fwd));\n"
484  << " DrawUnfoldedSummary(file);\n"
485  << "}\n"
486  << "// EOF" << std::endl;
487  f.close();
488  }
503  void ProcessType(TCollection* measured,
504  TCollection* corrections,
505  UInt_t method,
506  Double_t regParam,
507  TDirectory* out,
508  UShort_t sys,
509  UShort_t sNN)
510  {
511  Printf(" Processing %s ...", measured->GetName());
512  TDirectory* dir = out->mkdir(measured->GetName());
513 
514  // Make some summary stacks
515  THStack* allMeasured = new THStack("measured",
516  "Measured P(#it{N}_{ch})");
517  THStack* allTruth = new THStack("truth",
518  "MC 'truth' P(#it{N}_{ch})");
519  THStack* allTruthA = new THStack("truthAccepted",
520  "MC 'truth' accepted P(#it{N}_{ch})");
521  THStack* allUnfolded = new THStack("unfolded",
522  "Unfolded P(#it{N}_{ch})");
523  THStack* allCorrected = new THStack("corrected",
524  "Corrected P(#it{N}_{ch})");
525  THStack* allRatio = (sys != 1 ? 0 :
526  new THStack("ratios", "Ratios to other"));
527  TMultiGraph* allALICE = (sys != 1 ? 0 :
528  new TMultiGraph("alice", "ALICE Published"));
529  TMultiGraph* allCMS = (sys != 1 ? 0 :
530  new TMultiGraph("cms", "CMS Published"));
531 
532  // Loop over the list of objects.
533  static TRegexp regex("[pm][0-9]d[0-9]*_[pm][0-9]d[0-9]*");
534  TIter next(measured);
535  TObject* o = 0;
536  Int_t i = 0;
537  Double_t r = regParam;
538  while ((o = next())) {
539  // Go back to where we where
540  dir->cd();
541 
542  // if not a collection, don't bother
543  if (!o->IsA()->InheritsFrom(TCollection::Class())) continue;
544 
545  // If it doesn't match our regular expression, don't bother
546  TString n(o->GetName());
547  if (n.Index(regex) == kNPOS) {
548  // Warning("ScanType", "%s in %s doesn't match eta range regexp",
549  // n.Data(), real->GetName());
550  continue;
551  }
552  TCollection* mBin = static_cast<TCollection*>(o);
553  TCollection* cBin = GetCollection(corrections, n.Data());
554  if (!cBin) continue;
555 
556  THStack* binS = ProcessBin(mBin, cBin, method, r, dir);
557  if (!binS) continue;
558 
559  TH1* result = 0;
560  Bin2Stack(binS, i, allMeasured, allTruth, allTruthA,
561  allUnfolded, allCorrected, result);
562 
563  TGraph* alice = 0;
564  TGraph* cms = 0;
565  Other2Stack(o->GetName(), i, sNN, allALICE, allCMS, alice, cms);
566  Ratio2Stack(i, result, alice, cms, allRatio);
567  i++;
568  }
569  dir->Add(allMeasured);
570  dir->Add(allTruth);
571  dir->Add(allTruthA);
572  dir->Add(allUnfolded);
573  dir->Add(allCorrected);
574  if (allALICE && allALICE->GetListOfGraphs()) {
575  if (allALICE->GetListOfGraphs()->GetEntries() > 0)
576  dir->Add(allALICE);
577  else
578  delete allALICE;
579  }
580  if (allCMS && allCMS->GetListOfGraphs()) {
581  if (allCMS->GetListOfGraphs()->GetEntries() > 0)
582  dir->Add(allCMS);
583  else
584  delete allCMS;
585  }
586  if (allRatio && allRatio->GetHists()) {
587  if (allRatio->GetHists()->GetEntries() > 0)
588  dir->Add(allRatio);
589  else
590  delete allRatio;
591  }
592  }
604  THStack* ProcessBin(TCollection* measured,
605  TCollection* corrections,
606  UInt_t method,
607  Double_t regParam,
608  TDirectory* out)
609  {
610  Printf(" Processing %s ...", measured->GetName());
611  // Try to get the data
612  TH1* inRaw = GetH1(measured, "rawDist");
613  TH1* inTruth = GetH1(corrections, "truth");
614  TH1* inTruthA = GetH1(corrections, "truthAccepted");
615  TH1* inTrgVtx = GetH1(corrections, "triggerVertex");
616  TH2* inResp = GetH2(corrections, "response");
617  if (!inRaw || !inTruth || !inTruthA || !inTrgVtx || !inResp)
618  return 0;
619 
620  // Make output directory
621  TDirectory* dir = out->mkdir(measured->GetName());
622  dir->cd();
623 
624  // Copy the input to the output
625  TH1* outRaw = static_cast<TH1*>(inRaw ->Clone("measured"));
626  TH1* outTruth = static_cast<TH1*>(inTruth ->Clone("truth"));
627  TH1* outTruthA = static_cast<TH1*>(inTruthA ->Clone("truthAccepted"));
628  TH1* outTrgVtx = static_cast<TH1*>(inTrgVtx ->Clone("triggerVertex"));
629  TH2* outResp = static_cast<TH2*>(inResp ->Clone("response"));
630 
631  // Make our response matrix
632  RooUnfoldResponse matrix(0, 0, inResp);
633 
634  // Store regularization parameter
635  Double_t r = regParam;
636  RooUnfold::Algorithm algo = (RooUnfold::Algorithm)method;
637  RooUnfold* unfolder = RooUnfold::New(algo, &matrix, inRaw, r);
638  unfolder->SetVerbose(0);
639 
640  // Do the unfolding and get the result
641  TH1* res = unfolder->Hreco();
642  res->SetDirectory(0);
643 
644  // Make a copy to store on the output
645  TH1* outUnfold = static_cast<TH1*>(res->Clone("unfolded"));
646  TString tit(outUnfold->GetTitle());
647  tit.ReplaceAll("Unfold Reponse matrix", "Unfolded P(#it{N}_{ch})");
648  outUnfold->SetTitle(tit);
649 
650  // Clone the unfolded results and divide by the trigger/vertex
651  // bias correction
652  TH1* outCorr = static_cast<TH1*>(outUnfold->Clone("corrected"));
653  outCorr->Divide(inTrgVtx);
654  tit.ReplaceAll("Unfolded", "Corrected");
655  outCorr->SetTitle(tit);
656 
657  // Now normalize the output to integral=1
658  TH1* hists[] = { outRaw, outUnfold, outCorr, 0 };
659  TH1** phist = hists;
660  while (*phist) {
661  TH1* h = *phist;
662  if (h) {
663  Double_t intg = h->Integral(1, h->GetXaxis()->GetXmax());
664  h->Scale(1. / intg, "width");
665  }
666  phist++;
667  }
668 
669  // And make ratios
670  TH1* ratioTrue = static_cast<TH1*>(outCorr->Clone("ratioCorrTruth"));
671  tit = ratioTrue->GetTitle();
672  tit.ReplaceAll("Corrected", "Corrected/MC 'truth'");
673  ratioTrue->SetTitle(tit);
674  ratioTrue->Divide(outTruth);
675  ratioTrue->SetYTitle("P_{corrected}(#it{N}_{ch})/P_{truth}(#it{N}_{ch})");
676 
677  TH1* ratioAcc = static_cast<TH1*>(outUnfold->Clone("ratioUnfAcc"));
678  tit = ratioAcc->GetTitle();
679  tit.ReplaceAll("Unfolded", "Unfolded/MC selected");
680  ratioAcc->SetTitle(tit);
681  ratioAcc->Divide(outTruthA);
682  ratioAcc->SetYTitle("P_{unfolded}(#it{N}_{ch})/P_{MC}(#it{N}_{ch})");
683 
684 
685  // Make a stack
686  tit = measured->GetName();
687  tit.ReplaceAll("m", "-");
688  tit.ReplaceAll("p", "+");
689  tit.ReplaceAll("d", ".");
690  tit.ReplaceAll("_", "<#it{#eta}<");
691  THStack* stack = new THStack("all", tit);
692  stack->Add(outTruth, "E2");
693  stack->Add(outTruthA, "E2");
694  stack->Add(outRaw, "E1");
695  stack->Add(outUnfold, "E1");
696  stack->Add(outCorr, "E1");
697  dir->Add(stack);
698 
699  // Rest of the function is devoted to making the output look nice
700  outRaw ->SetDirectory(dir);
701  outTruth ->SetDirectory(dir);
702  outTruthA->SetDirectory(dir);
703  outTrgVtx->SetDirectory(dir);
704  outResp ->SetDirectory(dir);
705  outUnfold->SetDirectory(dir);
706  outCorr ->SetDirectory(dir);
707 
708  outRaw ->SetMarkerStyle(20); // Measured is closed
709  outTruth ->SetMarkerStyle(24); // MC is open
710  outTruthA->SetMarkerStyle(24); // MC is open
711  outTrgVtx->SetMarkerStyle(20); // Derived is closed
712  outUnfold->SetMarkerStyle(20); // Derived is closed
713  outCorr ->SetMarkerStyle(20); // Derived is closed
714 
715  outRaw ->SetMarkerSize(0.9);
716  outTruth ->SetMarkerSize(1.6);
717  outTruthA->SetMarkerSize(1.4);
718  outTrgVtx->SetMarkerSize(1.0);
719  outUnfold->SetMarkerSize(0.9);
720  outCorr ->SetMarkerSize(1.0);
721 
722  outRaw ->SetMarkerColor(kColorMeasured);
723  outTruth ->SetMarkerColor(kColorTruth);
724  outTruthA->SetMarkerColor(kColorAccepted);
725  outTrgVtx->SetMarkerColor(kColorTrgVtx);
726  outUnfold->SetMarkerColor(kColorUnfolded);
727  outCorr ->SetMarkerColor(kColorCorrected);
728 
729  outRaw ->SetFillColor(kColorError);
730  outTruth ->SetFillColor(kColorError);
731  outTruthA->SetFillColor(kColorError);
732  outTrgVtx->SetFillColor(kColorError);
733  outUnfold->SetFillColor(kColorError);
734  outCorr ->SetFillColor(kColorError);
735 
736  outRaw ->SetFillStyle(0);
737  outTruth ->SetFillStyle(1001);
738  outTruthA->SetFillStyle(1001);
739  outTrgVtx->SetFillStyle(0);
740  outUnfold->SetFillStyle(0);
741  outCorr ->SetFillStyle(0);
742 
743  outRaw ->SetLineColor(kBlack);
744  outTruth ->SetLineColor(kBlack);
745  outTruthA->SetLineColor(kBlack);
746  outTrgVtx->SetLineColor(kBlack);
747  outUnfold->SetLineColor(kBlack);
748  outCorr ->SetLineColor(kBlack);
749 
750  // Legend
751  TLegend* l = StackLegend(stack);
752  l->AddEntry(outRaw, "Raw", "lp");
753  l->AddEntry(outTruth, "MC 'truth'", "fp");
754  l->AddEntry(outTruthA, "MC 'truth' accepted", "fp");
755  l->AddEntry(outUnfold, "Unfolded", "lp");
756  l->AddEntry(outCorr, "Corrected", "lp");
757 
758  return stack;
759  }
760  static void BinAttributes(Int_t i,
761  Int_t& open,
762  Int_t& closed,
763  Float_t& size,
764  Double_t& factor)
765  {
766  // --- Setup for markers -----------------------------------------
767  const Int_t nMarkers = 7;
768  const Int_t aClosed[] = { 20, 21, 22, 23, 29, 33, 34 };
769  const Int_t aOpen[] = { 24, 25, 26, 32, 30, 27, 28 };
770  const Float_t aSize[] = { 1.1, 1.0, 1.2, 1.2, 1.2, 1.2, 1.0 };
771  Int_t j = i % nMarkers;
772 
773  open = aOpen[j];
774  closed = aClosed[j];
775  size = aSize[j];
776  factor = TMath::Power(10, i);
777  }
790  void Bin2Stack(const THStack* bin, Int_t i,
791  THStack* measured,
792  THStack* truth,
793  THStack* accepted,
794  THStack* unfolded,
795  THStack* corrected,
796  TH1*& result)
797  {
798  Int_t open, closed;
799  Double_t factor;
800  Float_t size;
801  BinAttributes(i, open, closed, size, factor);
802 
803  TIter next(bin->GetHists());
804  TH1* h = 0;
805  while ((h = static_cast<TH1*>(next()))) {
806  THStack* tmp = 0;
807  Int_t col = h->GetMarkerColor();
808  Int_t sty = 0;
809  switch (col) {
810  case kColorMeasured: tmp = measured; sty = closed; break;
811  case kColorTruth: tmp = truth; sty = open; break;
812  case kColorAccepted: tmp = accepted; sty = open; break;
813  case kColorUnfolded: tmp = unfolded; sty = closed; break;
814  case kColorCorrected: tmp = corrected; sty = closed; break;
815  default: continue;
816  }
817  // Now clone, and add to the appropriate stack
818  TH1* cln = static_cast<TH1*>(h->Clone(h->GetName()));
819  cln->SetDirectory(0);
820  cln->SetMarkerStyle(sty);
821  cln->SetMarkerSize(size);
822  cln->Scale(factor); // Scale by 10^i
823  if (col == kColorCorrected) result = cln;
824 
825  // Make sure we do not get the old legend
826  TObject* tst = cln->FindObject("legend");
827  if (tst) cln->GetListOfFunctions()->Remove(tst);
828 
829  tmp->Add(cln, next.GetOption());
830  }
831 
832  // Add entries to our stacks
833  TString txt = bin->GetTitle();
834  if (i == 0) txt.Append(" (#times1)");
835  else if (i == 1) txt.Append(" (#times10)");
836  else txt.Append(Form(" (#times10^{%d})", i));
837  THStack* stacks[] = { measured, truth, accepted, unfolded, corrected, 0 };
838  THStack** pstack = stacks;
839  while (*pstack) {
840  TLegend* leg = StackLegend(*pstack);
841  pstack++;
842  if (!leg) continue;
843 
844  TObject* dummy = 0;
845  TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
846  e->SetMarkerStyle(closed);
847  e->SetMarkerSize(1.2*size);
848  e->SetMarkerColor(kBlack);
849  e->SetFillColor(0);
850  e->SetFillStyle(0);
851  e->SetLineColor(kBlack);
852  }
853  }
865  void Other2Stack(const TString& name, Int_t i,
866  UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS,
867  TGraph*& alice, TGraph*& cms)
868  {
869  if (!allALICE && !allCMS) return;
870 
871  TString tmp(name);
872  tmp.ReplaceAll("p", "+");
873  tmp.ReplaceAll("m", "-");
874  tmp.ReplaceAll("_", " ");
875  tmp.ReplaceAll("d", ".");
876  TObjArray* tokens = tmp.Tokenize(" ");
877  if (!tokens || tokens->GetEntriesFast() < 2) {
878  Error("Other2Stack", "Failed to decode eta range from %s", name.Data());
879  if (tokens) tokens->Delete();
880  return;
881  }
882  Double_t eta1 = static_cast<TObjString*>(tokens->At(0))->String().Atof();
883  Double_t eta2 = static_cast<TObjString*>(tokens->At(1))->String().Atof();
884  tokens->Delete();
885 
886  if (TMath::Abs(eta2+eta1) > 1e-3) {
887  // Not symmetric bin
888  // Info("Other2Stack", "bin [%f,%f] is not symmetric (%f)",
889  // eta1, eta2, TMath::Abs(eta2-eta1));
890  return;
891  }
892  Double_t aEta = TMath::Abs(eta1);
893 
894  Int_t open, closed;
895  Double_t factor;
896  Float_t size;
897  BinAttributes(i, open, closed, size, factor);
898 
899  if (allALICE) {
900  TGraphAsymmErrors* g = GetOther(1, aEta, sNN, factor);
901  if (g) {
902  g->SetMarkerStyle(closed);
903  g->SetMarkerColor(kColorALICE);
904  g->SetMarkerSize(size);
905  allALICE->Add(g, "p same");
906  alice = g;
907  }
908  }
909  if (allCMS) {
910  TGraphAsymmErrors* g = GetOther(0, aEta, sNN, factor);
911  if (g) {
912  g->SetMarkerStyle(closed);
913  g->SetMarkerColor(kColorCMS);
914  g->SetMarkerSize(size);
915  allCMS->Add(g, "p same");
916  cms = g;
917  }
918  }
919  }
929  void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
930  {
931  if (!all || !res || !(alice || cms)) return;
932 
933  Int_t off = 5*ib;
934  TGraph* gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
935  TGraph** pg = gs;
936  while (*pg) {
937  TGraph* g = *pg;
938  const char* n = (g == alice ? "ALICE" : "CMS");
939 
940  TH1* r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
941  TString tit(r->GetTitle());
942  tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
943  r->SetTitle(tit);
944  r->SetMarkerColor(g->GetMarkerColor());
945  r->SetLineColor(g->GetLineColor());
946 
947  TObject* tst = r->FindObject("legend");
948  if (tst) r->GetListOfFunctions()->Remove(tst);
949 
950  for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
951  Double_t c = r->GetBinContent(i);
952  Double_t e = r->GetBinError(i);
953  Double_t o = g->Eval(r->GetBinCenter(i));
954  if (o < 1e-12) {
955  r->SetBinContent(i, 0);
956  r->SetBinError(i, 0);
957  continue;
958  }
959  r->SetBinContent(i, (c - o) / o + off);
960  r->SetBinError(i, e / o);
961  }
962  all->Add(r);
963  pg++;
964  }
965  TLegend* leg = StackLegend(all);
966  if (!leg) return;
967 
968  TString txt = res->GetTitle();
969  txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
970  if (ib == 0) txt.Append(" "); // (#times1)");
971  // else if (ib == 1) txt.Append(" (#times10)");
972  else txt.Append(Form(" (+%d)", off));
973 
974  TObject* dummy = 0;
975  TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
976  e->SetMarkerStyle(res->GetMarkerStyle());
977  e->SetMarkerSize(res->GetMarkerSize());
978  e->SetMarkerColor(kBlack);
979  e->SetFillColor(0);
980  e->SetFillStyle(0);
981  e->SetLineColor(kBlack);
982  }
983 
993  TLegend* StackLegend(THStack* stack)
994  {
995  TList* hists = stack->GetHists();
996  if (!hists) return 0;
997 
998  TObject* first = hists->First();
999  if (!first) return 0;
1000 
1001  TH1* hist = static_cast<TH1*>(first);
1002  TList* list = hist->GetListOfFunctions();
1003  TObject* o = list->FindObject("legend");
1004  if (o) return static_cast<TLegend*>(o);
1005 
1006  TLegend* l = new TLegend(0.65, 0.65, 0.9, 0.9, "", "NDC");
1007  l->SetName("legend");
1008  l->SetBorderSize(0);
1009  l->SetFillColor(0);
1010  l->SetFillStyle(0);
1011  l->SetTextFont(42);
1012  list->Add(l);
1013 
1014  return l;
1015  }
1016 
1017  /* =================================================================
1018  *
1019  * Measurements from other sources, such as published ALICE, CMS, ...
1020  */
1022  Int_t factor)
1023  {
1024  TString oC = Form("OtherPNch::GetData(%hu,%f,%hu);",
1025  type, eta, sNN);
1026  TGraphAsymmErrors* g =
1027  reinterpret_cast<TGraphAsymmErrors*>(gROOT->ProcessLine(oC));
1028  if (!g) {
1029  // Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
1030  // type, eta, sNN);
1031  return 0;
1032  }
1033 
1034 
1035  for (Int_t j = 0; j < g->GetN(); j++) {
1036  g->SetPoint(j, g->GetX()[j], g->GetY()[j]*factor);
1037  g->SetPointError(j, g->GetEXlow()[j], g->GetEXhigh()[j],
1038  g->GetEYlow()[j]*factor, g->GetEYhigh()[j]*factor);
1039  }
1040  return g;
1041  }
1042 };
1043 
1044 void
1046  Double_t regParam=1e-30,
1047  const TString& measured="forward_multdists.root",
1048  const TString& corrections="")
1049 {
1050  Unfolder m;
1051  m.Run(measured, corrections, method, regParam);
1052 }
void Run(const TString &measuredFile, const TString &corrFile, const TString &method="Bayes", Double_t regParam=4)
double Double_t
Definition: External.C:58
void Ratio2Stack(Int_t ib, TH1 *res, TGraph *alice, TGraph *cms, THStack *all)
void Bin2Stack(const THStack *bin, Int_t i, THStack *measured, THStack *truth, THStack *accepted, THStack *unfolded, THStack *corrected, TH1 *&result)
TString fileName
static TCollection * GetTop(const TString &fileName, Bool_t results=false)
TList * list
TDirectory file where lists per trigger are stored in train ouput.
static void GetParameter(TCollection *c, const TString &name, ULong_t &v)
static TH2 * GetH2(TCollection *c, const TString &name, Bool_t verbose=true)
TCanvas * c
Definition: TestFitELoss.C:172
static void GetParameter(TCollection *c, const TString &name, UShort_t &v)
static void AppendAnd(TString &trg, const TString &what)
void UnfoldMultDists(const TString &method="Bayes", Double_t regParam=1e-30, const TString &measured="forward_multdists.root", const TString &corrections="")
void Other2Stack(const TString &name, Int_t i, UShort_t sNN, TMultiGraph *allALICE, TMultiGraph *allCMS, TGraph *&alice, TGraph *&cms)
static TH1 * GetH1(TCollection *c, const TString &name, Bool_t verbose=true)
static TObject * GetObject(TCollection *c, const TString &name, TClass *cl, Bool_t verbose=true)
int Int_t
Definition: External.C:63
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, Bool_t self) const
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Int_t method
unsigned long ULong_t
Definition: External.C:38
static void GetParameter(TCollection *c, const TString &name, Double_t &v)
TLegend * StackLegend(THStack *stack)
static TCollection * GetCollection(TCollection *c, const TString &name, Bool_t verbose=-true)
Definition: External.C:220
TFile * file
TList with histograms for a given trigger.
void ProcessType(TCollection *measured, TCollection *corrections, UInt_t method, Double_t regParam, TDirectory *out, UShort_t sys, UShort_t sNN)
static void BinAttributes(Int_t i, Int_t &open, Int_t &closed, Float_t &size, Double_t &factor)
unsigned short UShort_t
Definition: External.C:28
void SaveSummarize()
THStack * ProcessBin(TCollection *measured, TCollection *corrections, UInt_t method, Double_t regParam, TDirectory *out)
TGraphAsymmErrors * GetOther(UShort_t type, Double_t eta, UShort_t sNN, Int_t factor)
bool Bool_t
Definition: External.C:53
static UInt_t MethodId(TString &method)
Definition: External.C:196
TDirectoryFile * dir