AliPhysics  cdeda5a (cdeda5a)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UnfoldMult.C
Go to the documentation of this file.
1 
11 #include <iostream>
12 #include <TList.h>
13 #include <TFile.h>
14 #include <TObject.h>
15 #include <TDirectory.h>
16 #include "TH1D.h"
17 #include "TH2D.h"
18 #include <TROOT.h>
19 #include <TSystem.h>
20 #include "RooUnfoldBayes.h"
21 #include "RooUnfoldResponse.h"
22 #include <TMath.h>
38 TObject* GetObject(TDirectory* d, const char* name)
39 {
40  TObject* o = d->Get(name);
41  if (!o) {
42  Error("GetCollection", "%s not found in %s", name, d->GetName());
43  return 0;
44  }
45  return o;
46 }
57 TObject* GetObject(TCollection* d, const char* name)
58 {
59  TObject* o = d->FindObject(name);
60  if (!o) {
61  Error("GetCollection", "%s not found in %s", name, d->GetName());
62  d->ls();
63  return 0;
64  }
65  return o;
66 }
77 TCollection* GetCollection(TDirectory* d, const char* name)
78 {
79  return static_cast<TCollection*>(GetObject(d, name));
80 }
91 TCollection* GetCollection(TCollection* d, const char* name)
92 {
93  return static_cast<TCollection*>(GetObject(d, name));
94 }
105 TH1* GetH1(TDirectory* d, const char* name)
106 {
107  return static_cast<TH1*>(GetObject(d, name));
108 }
119 TH1* GetH1(TCollection* d, const char* name)
120 {
121  return static_cast<TH1*>(GetObject(d, name));
122 }
131 TH2* GetH2(TDirectory* d, const char* name)
132 {
133  return static_cast<TH2*>(GetObject(d, name));
134 }
145 TH2* GetH2(TCollection* d, const char* name)
146 {
147  return static_cast<TH2*>(GetObject(d, name));
148 }
159 TFile* OpenFile(const char* filename, Bool_t readNotWrite)
160 {
161  TFile* ret = TFile::Open(filename, readNotWrite ? "READ" : "RECREATE");
162  if (!ret) {
163  Error("OpenFile", "Failed to open the file \"%s\"", filename);
164  return 0;
165  }
166  return ret;
167 }
173 typedef TCollection Dir;
174 
187 void GetHistos(Double_t lowLim,
188  Double_t highLim,
189  Dir* resp,
190  Dir* data,
191  TDirectory* out);
202 void ProcessUnfold(TH2* response, TH1* measured, TH1* mcHist, TH1* esdHist);
213 TH1* GetTriggerBias(TH1* mcHist, TH1* esdHist);
221 void DoCorrection(TH1* unfoldBefore, TH1* triggerBias, TH1* mcHist);
222 
223 
233 void UnfoldMult(const char* respFile="forward_response.root",
234  const char* dataFile="forward_multiplicy.root",
235  const char* outFile="unfolded.root")
236 {
237  TFile* output = OpenFile(outFile,false);
238  TFile* resp = OpenFile(respFile, true);
239  TFile* data = OpenFile(dataFile, true);
240  if (!output || !resp || !data) return;
241 
242  Dir* topResp = GetCollection(resp, "ResponseMatricesSums");
243  Dir* topData = GetCollection(data, "MultSums");
244  if (!topData || !topResp) return;
245 
246  Double_t limits[] = {5.1, 3.4, 3.0, 2.4, 2.0, 1.5, 1.0, 0.5, 0. };
247  Double_t* limit = limits;
248 
249  while ((*limit) > 0.1) {
250  if ((*limit) >5.) GetHistos(-3.4, (*limit), topResp, topData, output);
251  else GetHistos(-(*limit), (*limit), topResp, topData, output);
252  limit++;
253  output->cd();
254  }
255  output->Close();
256  resp->Close();
257  data->Close();
258 }
259 
260 //____________________________________________________________________
261 void GetHistos(Double_t lowLim,
262  Double_t highLim,
263  Dir* topResp,
264  Dir* topData,
265  TDirectory* out)
266 {
267  // Format the name of the bin
268  TString sLow( TString::Format("%+03d",int(10*lowLim)));
269  TString sHigh(TString::Format("%+03d",int(10*highLim)));
270  sLow .ReplaceAll("+", "plus");
271  sLow .ReplaceAll("-", "minus");
272  sHigh.ReplaceAll("+", "plus");
273  sHigh.ReplaceAll("-", "minus");
274  TString name = TString::Format("%s_%s", sLow.Data(), sHigh.Data());
275  TDirectory* dir = out->mkdir(name);
276  dir->cd();
277 
278  // Get our collections
279  Dir* listResp = GetCollection(topResp, name);
280  Dir* listData = GetCollection(topData,name);
281  if(!listResp || !listData) return;
282 
283  // Get the MC histograms
284  TH2* response = GetH2(listResp, "responseMatrix");
285  TH1* mcHist = GetH1(listResp, "fMCNSD");
286  TH1* esdHist = GetH1(listResp, "fESDNSD");
287  if (!response || !mcHist || !esdHist) return;
288 
289  // Get the data histogram
290  TH1* measured = GetH1(listData, "mult");
291  if(!measured) return;
292 
293  // Now do the unfolding
294  ProcessUnfold(response, measured, mcHist, esdHist);
295 }
296 
297 
298 
299 //____________________________________________________________________
300 void ProcessUnfold(TH2* response, TH1* measured, TH1* mcHist, TH1* esdHist)
301 {
302  Int_t nX = response->GetNbinsX();
303  Int_t nY = response->GetNbinsY();
304  TH1* projY = static_cast<TH1*>(response->ProjectionY("projY",1,nX,""));
305  TH1* projX = static_cast<TH1*>(response->ProjectionX("projX",1,nY,""));
306  projX->SetDirectory(0);
307  projY->SetDirectory(0);
308 
309  TH2* responseTranspose = static_cast<TH2*>(response->Clone("response"));
310  for(Int_t i = 1; i <= nX; i++) {
311  for(Int_t j = 1; j <= nY; j++) {
312  responseTranspose->SetBinContent(j,i, response->GetBinContent(i,j));
313  responseTranspose->SetBinError(j,i, response->GetBinError(i,j));
314  }
315  }
316 
317  RooUnfoldResponse* responseObject = new RooUnfoldResponse(projY, projX,
318  responseTranspose,
319  "", "");
320  RooUnfold* unfold = new RooUnfoldBayes(responseObject,
321  measured, 5);
322  TH1* unfolded = static_cast<TH1*>(unfold->Hreco());
323  TH1* triggerBias = GetTriggerBias(mcHist, esdHist);
324  DoCorrection(unfolded, TriggerBias, mcHist);
325 
326  Double_t scale_unf = 1/unfolded->Integral();
327  unfolded->Scale(scale_unf);
328  unfolded->Write("Unfolded");
329 
330  Double_t scale_raw = 1/measured->Integral();
331  measured->Scale(scale_raw);
332  measured->Write("Raw");
333  triggerBias->Write("TriggerBiasCorr");
334 }
335 
336 
337 
338 //____________________________________________________________________
339 TH1* GetTriggerBias(TH1* mcHist, TH1* esdHist)
340 {
341  mcHist->Sumw2();
342  esdHist->Sumw2();
343 
344  TH1* corr = new TH1D("corr", "corr", 40, -0.5, 39.5);
345  for (Int_t j=1; j<corr->GetNbinsX(); j++) {
346  Double_t errSqr = 0.;
347  Double_t cESD = esdHist->GetBinContent(j);
348  Double_t cMS = mcHist->GetBinContent(j);
349  Double_t c = (cMC == 0 ? 1 :
350  cESD == 0 ? 1/cMC : cESD/cMC);
351  corr->SetBinContent(j, c);
352  corr->SetBinError(j, c*TMath::Sqrt(errSqr));
353  }
354  return corr;
355 }
356 
357 //____________________________________________________________________
358 void DoCorrection(TH1* unfoldBefore, TH1* triggerBias, TH1*)
359 {
360  for (Int_t k=1; k<=35; k++) {
361  Double_t tb = triggerBias->GetBinContent(k);
362  Double_t ub = unfoldBefore->GetBinContent(k);
363  if (tb > 1e-5 && ub > 0) {
364  unfoldBefore->SetBinContent(k, ub / tb);
365  Double_t eub = UnfoldBefore->GetBinError(k);
366  Double_t etb = TriggerBias->GetBinError(k);
367  unfoldBefore->SetBinError(k, TMath::Sqrt(TMath::Power(eub/ub)+
368  TMath::Power(etb/tb))*ub/tb);
369  }
370  }
371 }
372 //
373 // EOF
374 //
375 
376 
377 
const char * filename
Definition: TestFCM.C:1
double Double_t
Definition: External.C:58
TCanvas * c
Definition: TestFitELoss.C:172
void UnfoldMult(const char *respFile="forward_response.root", const char *dataFile="forward_multiplicy.root", const char *outFile="unfolded.root")
Definition: UnfoldMult.C:233
TFile * OpenFile(const char *filename, Bool_t readNotWrite)
Definition: UnfoldMult.C:159
void GetHistos(Double_t lowLim, Double_t highLim, Dir *resp, Dir *data, TDirectory *out)
Definition: UnfoldMult.C:261
TH1 * GetTriggerBias(TH1 *mcHist, TH1 *esdHist)
Definition: UnfoldMult.C:339
int Int_t
Definition: External.C:63
Definition: External.C:212
TCollection Dir
Definition: UnfoldMult.C:173
TH1 * GetH1(TDirectory *d, const char *name)
Definition: UnfoldMult.C:105
void ProcessUnfold(TH2 *response, TH1 *measured, TH1 *mcHist, TH1 *esdHist)
Definition: UnfoldMult.C:300
TH2 * GetH2(TDirectory *d, const char *name)
Definition: UnfoldMult.C:131
TCollection * GetCollection(TDirectory *d, const char *name)
Definition: UnfoldMult.C:77
Definition: External.C:220
void DoCorrection(TH1 *unfoldBefore, TH1 *triggerBias, TH1 *mcHist)
Definition: UnfoldMult.C:358
TObject * GetObject(TDirectory *d, const char *name)
Definition: UnfoldMult.C:38
bool Bool_t
Definition: External.C:53
Definition: External.C:196