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