AliPhysics  95775ff (95775ff)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MakeDeltaWeights.C
Go to the documentation of this file.
1 
11 #include <TSystem.h>
12 #include "AliTrackletAODUtils.C"
13 #include "AliTrackletWeights.C"
14 
21 {
43  AliTrackletDeltaWeights* Run(const char* realFileName,
44  const char* simFileName,
45  Int_t dimen=2,
46  Bool_t scaleToTail=true)
47  {
48  TFile* realFile = 0;
49  TFile* simFile = 0;
50  if (!(realFile = OpenFile(realFileName))) return 0;
51  if (!(simFile = OpenFile(simFileName))) return 0;
52 
53  TString outName(gSystem->BaseName(simFileName));
54  const char* base = "MidRapidity";
55  Container* realTop = GetC(realFile, Form("%sResults", base));
56  Container* simTop = GetC(simFile, Form("%sMCResults", base));
57  if (!realTop) {
58  realTop = GetC(realFile, Form("%sMCResults", base));
59  if (realTop)
60  Warning("Run","\n"
61  "*********************************************\n"
62  "* Testing MC vs. MC: *\n"
63  "* 'Data' file: %23s *\n"
64  "* Simulation file: %23s *\n"
65  "*********************************************\n",
66  realFileName, simFileName);
67  // fRealIsSim = true;
68  outName.Append("_");
69  outName.Append(gSystem->BaseName(realFileName));
70  outName.ReplaceAll(".root","");
71  }
72  outName.ReplaceAll(".root","");
73  outName.Prepend("delta_");
74  outName.Append(".root");
75 
76  TH1* realCent = GetH1(realTop, "cent");
77  TH1* simCent = GetH1(simTop, "cent");
78  TH1* realIPz = GetH1(realTop, "ipz");
79  TH1* simIPz = GetH1(simTop, "ipz");
80 
81  // Check consistency of found histograms
82  if (!CheckConsistency(realCent, simCent)) {
83  Warning("Post", "Centrality bins are incompatible, giving up");
84  return 0;
85  }
86  if (!CheckConsistency(realIPz, simIPz)) {
87  Warning("Post", "IPz bins are incompatible, giving up");
88  return 0;
89  }
90 
91  TFile* out = TFile::Open(outName, "RECREATE");
92  out->cd();
93 
94  AliTrackletDeltaWeights* weights =
95  new AliTrackletDeltaWeights("weights", "w_{#Delta}");
96  weights->SetCentAxis(*(realCent->GetXaxis()));
97  // Loop over defined centrality bins
98  for (Int_t i = 1; i <= realCent->GetNbinsX(); i++) {
99  Double_t c1 = realCent->GetXaxis()->GetBinLowEdge(i);
100  Double_t c2 = realCent->GetXaxis()->GetBinUpEdge (i);
101 
102  ProcessBin(i, dimen, scaleToTail, c1, c2, realTop, simTop, weights);
103  }
104 
105  weights->Write();
106  out->Write();
107  Printf("Output writtten to %s", outName.Data());
108 
109  return weights;
110  }
116  void Normalize(TH3* h)
117  {
118  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
119  for (Int_t j = 1; j <= h->GetNbinsZ(); j++) {
120  // Project out delta distribution for this bin
121  TH1* tmp = h->ProjectionY(Form("%s_%02d_%02d",
122  h->GetName(), i, j),
123  i, i, j, j);
124  Double_t intg = tmp->Integral();
125  tmp->SetDirectory(0);
126  // Scale to average
127  if (intg > 1e-3) tmp->Scale(tmp->GetNbinsX()*1./intg);
128  else tmp->Reset();
129  for (Int_t k = 1; k <= h->GetNbinsY(); k++) {
130  h->SetBinContent(i, k, j, tmp->GetBinContent(k));
131  h->SetBinError (i, k, j, tmp->GetBinError (k));
132  }
133  delete tmp;
134  }
135  }
136  }
146  TH1* Weights1(Container* realMeas,
147  Container* simMeas,
148  Bool_t scaleToTail)
149  {
150  TH1* realDelta = GetH1(realMeas, "delta");
151  TH1* simDelta = GetH1(simMeas, "delta");
152  if (scaleToTail) {
153  Double_t realTail = GetD(realMeas, "deltaTail");
154  Double_t realTailE = GetD(realMeas, "deltaTailError");
155  Double_t simTail = GetD(simMeas, "deltaTail");
156  Double_t simTailE = GetD(simMeas, "deltaTailError");
157  Double_t scaleE, scale = RatioE(realTail, realTailE,
158  simTail, simTailE,
159  scaleE);
160  Scale(simDelta, scale, scaleE);
161  }
162 
163  TH1* ratio = static_cast<TH1*>(realDelta->Clone("tmp"));
164  ratio->Divide(simDelta);
165  ratio->Smooth(3);
166 
167  return ratio;
168  }
178  TH2* Weights2(Container* realMeas,
179  Container* simMeas,
180  Bool_t scaleToTail)
181  {
182  TH2* realDelta = GetH2(realMeas, "etaDelta");
183  TH2* simDelta = GetH2(simMeas, "etaDelta");
184  if (scaleToTail) {
185  TH1* realTail = GetH1(realMeas, "etaDeltaTail");
186  TH1* simTail = GetH1(simMeas, "etaDeltaTail");
187  TH1* scale = static_cast<TH1*>(realTail->Clone("scale"));
188  scale->Divide(simTail);
189  scale->SetDirectory(0);
190  Scale(simDelta, scale);
191  delete scale;
192  }
193 
194  TH2* ratio = static_cast<TH2*>(realDelta->Clone("tmp"));
195  ratio->Divide(simDelta);
196  ratio->Smooth(1);
197 
198  return ratio;
199  }
209  TH3* Weights3(Container* realMeas,
210  Container* simMeas,
211  Bool_t scaleToTail)
212  {
213  TH3* realDelta = GetH3(realMeas, "etaDeltaIPz");
214  TH3* simDelta = GetH3(simMeas, "etaDeltaIPz");
215  if (scaleToTail) {
216  TH2* realTail = GetH2(realMeas, "etaIPzDeltaTail");
217  TH2* simTail = GetH2(simMeas, "etaIPzDeltaTail");
218  TH2* scale = static_cast<TH2*>(realTail->Clone("scale"));
219  scale->Divide(simTail);
220  scale->SetDirectory(0);
221  ScaleDelta(simDelta, scale);
222  delete scale;
223  }
224 
225  TH3* ratio = static_cast<TH3*>(realDelta->Clone("tmp"));
226  ratio->Divide(simDelta);
227 
228  return ratio;
229  }
245  Int_t dimen,
246  Bool_t scaleToTail,
247  Double_t c1,
248  Double_t c2,
249  Container* realTop,
250  Container* simTop,
251  AliTrackletDeltaWeights* weights)
252  {
253  Printf(" %5.2f-%5.2f%%", c1, c2);
254  TString name = CentName(c1, c2);
255  Color_t color = CentColor(bin);
256  Container* realBin = GetC(realTop, name);
257  Container* simBin = GetC(simTop, name);
258  if (!realBin || !simBin) return false;
259  Container* realMeas = GetC(realBin, "measured");
260  Container* simMeas = GetC(simBin, "measured");
261  if (!realMeas || !simMeas) return false;
262  TH1* ratio = 0;
263  switch (dimen) {
264  case 2:
265  ratio = Weights2(realMeas, simMeas, scaleToTail);
266  break;
267  case 3:
268  ratio = Weights3(realMeas, simMeas, scaleToTail);
269  break;
270  default:
271  ratio = Weights1(realMeas, simMeas, scaleToTail);
272  break;
273  }
274  if (!ratio) return false;
275 
276  ratio->SetName(name);
277  ratio->SetDirectory(0);
278  ratio->SetMarkerColor(color);
279  ratio->SetLineColor(color);
280  ratio->SetFillColor(color);
281  // Normalize(ratio);
282 
283  weights->SetHisto(bin, ratio);
284  return true;
285  }
286 };
287 
288 //
289 // EOF
290 //
291 
Int_t color[]
double Double_t
Definition: External.C:58
Definition: External.C:244
TH2 * Weights2(Container *realMeas, Container *simMeas, Bool_t scaleToTail)
Bool_t SetHisto(Int_t bin, TH1 *h)
static TFile * OpenFile(const char *filename)
TSystem * gSystem
static TH1 * GetH1(Container *parent, const char *name, Bool_t verb=true)
static const char * CentName(Double_t c1, Double_t c2)
static TH3 * ScaleDelta(TH3 *h, TH2 *scale)
TH1 * Weights1(Container *realMeas, Container *simMeas, Bool_t scaleToTail)
static Double_t GetD(Container *parent, const char *name, Double_t def=-1, Bool_t verb=true)
static Bool_t CheckConsistency(const TH1 *h1, const TH1 *h2)
static TH1 * Scale(TH1 *h, Double_t x, Double_t xe)
Utilities for midrapidity analysis.
void SetCentAxis(const TAxis &axis)
int Int_t
Definition: External.C:63
Encode simulation weights for 2nd pass.
static Color_t CentColor(const TAxis &axis, Double_t c1, Double_t c2)
static TH3 * GetH3(Container *parent, const char *name, Bool_t verb=true)
static Double_t RatioE(Double_t n, Double_t en, Double_t d, Double_t ed, Double_t &er)
static TH2 * GetH2(Container *parent, const char *name, Bool_t verb=true)
Definition: External.C:220
AliTrackletDeltaWeights * Run(const char *realFileName, const char *simFileName, Int_t dimen=2, Bool_t scaleToTail=true)
TH3 * Weights3(Container *realMeas, Container *simMeas, Bool_t scaleToTail)
static Container * GetC(Container *parent, const char *name, Bool_t verb=true)
bool Bool_t
Definition: External.C:53
Bool_t ProcessBin(Int_t bin, Int_t dimen, Bool_t scaleToTail, Double_t c1, Double_t c2, Container *realTop, Container *simTop, AliTrackletDeltaWeights *weights)
void Normalize(TH3 *h)
Definition: External.C:196