AliPhysics  5eaf189 (5eaf189)
PureMCWeights.C
Go to the documentation of this file.
1 
11 #include "AliTrackletAODUtils.C"
12 #ifndef __CINT__
13 #include "AliTrackletWeights.C"
14 #include <THStack.h>
15 #include <TFile.h>
16 #include <TH1.h>
17 #include <TH2.h>
18 #include <TF1.h>
19 #else
20 class THStack;
21 class TFile;
22 class TCanvas;
23 #endif
24 
32 {
41  const char* CentName(Double_t c1, Double_t c2)
42  {
43  static TString tmp;
44  tmp.Form("cent%06.2f_%06.2f", c1, c2);
45  tmp.ReplaceAll(".", "d");
46  return tmp.Data();
47  }
55  void Run(const char* nFile, const char* dFile, const char* oFile=0) {
56 
57  const Color_t cc[] = { kMagenta+2, // 0
58  kBlue+2, // 1
59  kAzure-1, // 2 // 10,
60  kCyan+2, // 3
61  kGreen+1, // 4
62  kSpring+5, // 5 //+10,
63  kYellow+1, // 6
64  kOrange+5, // 7 //+10,
65  kRed+1, // 8
66  kPink+5, // 9 //+10,
67  kBlack }; // 10
68  TString oN(oFile);
69  if (oN.IsNull()) oN.Form("%c2%c.root", dFile[0], nFile[0]);
70  Printf("********************************************************\n"
71  " Generating pure MC weights\n"
72  "\n"
73  " Numerator (truth) file: %s\n"
74  " Denominator (correction) file: %s\n"
75  " Output (weight) file: %s\n",
76  nFile, dFile, oN.Data());
77 
78 
79  TFile* nF = TFile::Open(nFile, "READ");
80  TFile* dF = TFile::Open(dFile, "READ");
81  if (!nF || !dF) return;
82 
83  Container* nT = GetC(nF, "MidRapidityMCResults");
84  Container* dT = GetC(dF, "MidRapidityMCResults");
85  if (!nT || !dT) return;
86 
87  TH1* nC = GetH1(nT, "cent");
88  TH1* dC = GetH1(dT, "cent");
89  if (!nC || !dC || !CheckConsistency(nC,dC)) return;
90 
91  TList* pdg = new TList;
92  TH2D* pt = 0;
93  pdg->SetOwner();
94 
95  for (Int_t b = 1; b <= nC->GetNbinsX(); b++) {
96  Double_t c1 = nC->GetXaxis()->GetBinLowEdge(b);
97  Double_t c2 = nC->GetXaxis()->GetBinUpEdge (b);
98  const char* cN = CentName(c1, c2);
99  Color_t cC = cc[b%10];
100  Container* nB = GetC(nT, cN);
101  Container* dB = GetC(dT, cN);
102  if (!nB || !dB) continue;
103 
104  Container* nG = GetC(nB, "generated");
105  Container* dG = GetC(dB, "generated");
106  if (!nG || !dG) continue;
107 
108  TH2* nPt2 = GetH2(nG, "etaPt");
109  TH2* dPt2 = GetH2(dG, "etaPt");
110  if (!nPt2 || !dPt2) continue;
111 
112  TH1* nPt = nPt2->ProjectionY(Form("n%s",cN)); nPt->SetDirectory(0);
113  TH1* dPt = dPt2->ProjectionY(Form("d%s",cN)); dPt->SetDirectory(0);
114  TH1* rPt = static_cast<TH1*>(nPt->Clone(Form("r%s", cN)));
115  rPt->Divide(dPt);
116  rPt->SetMarkerColor(cC);
117  rPt->SetDirectory(0);
118 
119  if (!pt) {
120  pt = static_cast<TH2D*>(Make2D(0, "pt", "pt weights",
121  kBlack, 20, *(nC->GetXaxis()),
122  *(nPt->GetXaxis())));
123  pt->SetDirectory(0);
124  }
125  for (Int_t p = 1; p <= rPt->GetNbinsX(); p++) {
126  pt->SetBinContent(b, p, rPt->GetBinContent(p));
127  pt->SetBinError (b, p, rPt->GetBinError (p));
128  }
129  // nPt->SetMarkerColor(kRed); pt->Add(nPt);
130  // dPt->SetMarkerColor(kBlue); pt->Add(dPt);
131  // pt->Add(nPt2->ProjectionY());
132  // pt->Add(dPt2->ProjectionY());
133  delete nPt;
134  delete dPt;
135 
136  TH2* nPdg2 = GetH2(nG, "etaPdg");
137  TH2* dPdg2 = GetH2(dG, "etaPdg");
138  if (!nPdg2 || !dPdg2) continue;
139  if (!CheckConsistency(nPdg2,dPdg2)) continue;
140 
141  for (Int_t p = 1; p <= nPdg2->GetNbinsY(); p++) {
142  TString sPdg = nPdg2->GetYaxis()->GetBinLabel(p);
143  Int_t iPdg = sPdg.Atoi();
144  if (iPdg < 0) continue;
145 
146  TH1* nPdg = nPdg2->ProjectionX("ntmp", p, p);
147  TH1* dPdg = dPdg2->ProjectionX("dtmp", p, p);
148  if (!nPdg || !dPdg) continue;
149 
150  TH1* rPdg =
151  static_cast<TH1*>(nPdg->Clone(Form("r%s_%s",cN,sPdg.Data())));
152  rPdg->Divide(dPdg);
153  rPdg->SetDirectory(0);
154  if (rPdg->GetEntries() < 1) {
155  delete nPdg;
156  delete dPdg;
157  delete rPdg;
158  continue;
159  }
160 
161  TF1* fPdg = new TF1("fPdg", "pol0", -5, .5);
162  rPdg->Fit(fPdg, "NQR", "", -.5, +.5);
163 
164  TH1* hPdg = static_cast<TH1*>(pdg->FindObject(sPdg));
165  if (!hPdg) {
166  TString tmp;
167  Color_t c;
168  Style_t s;
169  PdgAttr(iPdg, tmp, c, s);
170  hPdg = Make1D(pdg,sPdg, tmp, c, s, *(nC->GetXaxis()));
171  hPdg->SetBinContent(0, iPdg);
172  }
173  hPdg->SetBinContent(b, fPdg->GetParameter(0));
174  hPdg->SetBinError (b, fPdg->GetParError (0));
175 
176  delete nPdg;
177  delete dPdg;
178  delete rPdg;
179  delete fPdg;
180  }
181 
182  }
183  pt->Draw("colz");
184  // pdg->Draw("nostack");
185 
187  w->SetPtWeight(pt);
188  TIter next(pdg);
189  TH1D* hPdg = 0;
190  while ((hPdg = static_cast<TH1D*>(next()))) {
191  Int_t iPdg = hPdg->GetBinContent(0);
192  hPdg->SetBinContent(0, 0);
193  w->AddAbundanceWeight(iPdg, hPdg);
194  }
195  TCanvas* c = new TCanvas("c","c");
196  w->Draw();
197 
198  TFile* out = TFile::Open(oN, "RECREATE");
199  w->Write();
200  out->Write();
201 
202  }
203 };
Int_t pdg
static void PdgAttr(Int_t pdg, TString &nme, Color_t &c, Style_t &s)
static TH2 * Make2D(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis, const TAxis &yAxis)
const Color_t cc[]
Definition: DrawKs.C:1
double Double_t
Definition: External.C:58
Bool_t AddAbundanceWeight(Short_t pdg, const TH1D *h, UShort_t mode=0)
static TH1 * GetH1(Container *parent, const char *name, Bool_t verb=true)
const char * CentName(Double_t c1, Double_t c2)
Definition: PureMCWeights.C:41
TCanvas * c
Definition: TestFitELoss.C:172
void Draw(Option_t *option="")
static Bool_t CheckConsistency(const TH1 *h1, const TH1 *h2)
Utilities for midrapidity analysis.
int Int_t
Definition: External.C:63
Encode simulation weights for 2nd pass.
Definition: External.C:228
Definition: External.C:212
Bool_t SetPtWeight(const TH2D *h, UShort_t mode=0)
static TH2 * GetH2(Container *parent, const char *name, Bool_t verb=true)
Definition: External.C:220
static Container * GetC(Container *parent, const char *name, Bool_t verb=true)
static TH1 * Make1D(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis)
void Run(const char *nFile, const char *dFile, const char *oFile=0)
Definition: PureMCWeights.C:55
Definition: External.C:196