AliPhysics  63e47e1 (63e47e1)
AliAnaWeights.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // Root
17 #include <TSystem.h>
18 #include <TFile.h>
19 #include <TTree.h>
20 #include <TKey.h>
21 #include <TH1F.h>
22 #include <TList.h>
23 #include <TProfile.h>
24 
25 // AliRoot
26 #include "AliAnalysisManager.h"
27 #include "AliLog.h"
28 #include "AliGenPythiaEventHeader.h"
29 
30 #include "AliAnaWeights.h"
31 
33 ClassImp(AliAnaWeights) ;
35 
36 //______________________________________________________________
38 //______________________________________________________________
40 : TObject(), fDebug(0),
41  fhCentralityWeight(0),
42  fCentrality(0),
43  fUseCentralityWeight(0),
44  fDoMCParticlePtWeights(1),
45  fEtaFunction(0), fPi0Function(0),
46  fMCWeight(1.),
47  fCurrFileName(0),
48  fCheckMCCrossSection(kFALSE),
49  fJustFillCrossSecHist(0),
50  fhXsec(0),
51  fhTrials(0),
52  fPyEventHeader(0),
53  fCheckPythiaEventHeader(0)
54 {
55 }
56 
57 //______________________________________________________________
59 //______________________________________________________________
61 {
62  if ( fhCentralityWeight ) delete fhCentralityWeight ;
63 
64  if ( fEtaFunction ) delete fEtaFunction ;
65 
66  if ( fPi0Function ) delete fPi0Function ;
67 }
68 
69 //____________________________________________________
72 //____________________________________________________
74 {
75  TList * outputContainer = new TList() ;
76  outputContainer->SetName("MCWeightHistogram") ;
77 
78  if(!fCheckMCCrossSection) return outputContainer;
79 
80  outputContainer->SetOwner(kFALSE);
81 
82  fhXsec = new TH1F("hXsec","xsec from pyxsec.root",1,0,1);
83  fhXsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
84  outputContainer->Add(fhXsec);
85 
86  fhTrials = new TH1F("hTrials","trials root file",1,0,1);
87  fhTrials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
88  outputContainer->Add(fhTrials);
89 
90  return outputContainer ;
91 }
92 
93 //__________________________________
98 //_________________________________
100 {
101  Double_t weight = 1.;
102 
103  if ( fCheckMCCrossSection )
104  {
106 
107  AliDebug(1,Form("MC pT-hard weight: %e",temp));
108 
109  weight*=temp;
110  }
111 
112  if ( fUseCentralityWeight )
113  {
114  Double_t temp = fhCentralityWeight->GetBinContent(fhCentralityWeight->FindBin(fCentrality));
115 
116  AliDebug(1,Form("Centrality %2.1f, weight: %2.2f",fCentrality,temp));
117 
118  weight*=temp;
119  }
120 
121  AliDebug(1,Form("Event weight %e",weight));
122 
123  return weight ;
124 }
125 
126 //_____________________________________________
127 //
139 //_____________________________________________
141 {
142  Double_t weight = 1.;
143 
144  if ( !fDoMCParticlePtWeights ) return weight ;
145 
146  if (pdg == 111 && fPi0Function &&
147  (genName.Contains("Pi0") || genName.Contains("pi0") || genName.Contains("PARAM") || genName.Contains("BOX")))
148  weight = fPi0Function->Eval(pt);
149  else if (pdg == 221 && fEtaFunction &&
150  (genName.Contains("Eta") || genName.Contains("eta") || genName.Contains("PARAM") || genName.Contains("BOX")))
151  weight = fEtaFunction->Eval(pt);
152 
153  AliDebug(1,Form("MC particle pdg %d, pt %2.2f, generator %s with index %d: weight %f",pdg,pt,genName.Data(),igen, weight));
154 
155  return weight ;
156 }
157 
158 //_____________________________________________
159 //
168 //_____________________________________________
170 {
171  Float_t xsection = 0;
172  Float_t trials = 1;
173  Float_t avgTrials = 0;
174 
175  if ( !fhXsec || !fhTrials ) return 1;
176 
177  // -----------------------------------------------
178  // Check if info is already in Pythia event header
179  // Do not apply the weight per event, too much number of trial variation
180  // -----------------------------------------------
182  {
183  fMCWeight = 1 ;
184 
185  AliDebug(fDebug,Form("Pythia event header: xs %2.2e, trial %d", fPyEventHeader->GetXsection(),fPyEventHeader->Trials()));
186 
187  fhXsec ->Fill("<#sigma>" ,fPyEventHeader->GetXsection());
188  fhTrials->Fill("#sum{ntrials}",fPyEventHeader->Trials());
189 
190  return 1 ;
191  }
192 
193  // -----------------------------------------------
194  // Get cross section from corresponding files
195  // -----------------------------------------------
196 
197  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
198  if ( !tree ) return 1;
199 
200  TFile *curfile = tree->GetCurrentFile();
201 
202  if ( !curfile ) return 1;
203 
204  // Check if file not accessed previously, if so
205  // return the previously calculated weight
206  if(fCurrFileName == curfile->GetName()) return fMCWeight;
207 
208  fCurrFileName = TString(curfile->GetName());
209 
210  Bool_t ok = GetPythiaInfoFromFile(fCurrFileName,xsection,trials);
211 
212  if ( !ok )
213  {
214  AliWarning("Parameters from file not recovered properly");
215  return 1;
216  }
217 
218  fhXsec->Fill("<#sigma>",xsection);
219 
220  // average number of trials
221  Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
222 
223  if(trials >= nEntries && nEntries > 0.) avgTrials = trials/nEntries;
224 
225  fhTrials->Fill("#sum{ntrials}",avgTrials);
226 
227  AliInfo(Form("xs %2.2e, trial %e, avg trials %2.2e, events per file %e",
228  xsection,trials,avgTrials,nEntries));
229 
230  AliDebug(1,Form("Reading File %s",curfile->GetName()));
231 
232  if(avgTrials > 0.)
233  {
235  {
236  fMCWeight = xsection / avgTrials ;
237 
238  AliInfo(Form("MC Weight: %e",fMCWeight));
239  }
240  else fMCWeight = 1; // do not add weight to histograms
241  }
242  else
243  {
244  AliWarning(Form("Average number of trials is NULL!! Set weight to 1: xs : %e, trials %e, entries %e",
245  xsection,trials,nEntries));
246 
247  fMCWeight = 1;
248  }
249 
250  return fMCWeight ;
251 }
252 
253 //_______________________________________________________________________________________
260 //_______________________________________________________________________________________
262 {
263  xsec = 0;
264  trials = 1;
265 
266  if(file.Contains("root_archive.zip#"))
267  {
268  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
269  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
270  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
271  file.Replace(pos+1,pos2-pos1,"");
272  }
273  else
274  {
275  // not an archive take the basename....
276  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
277  }
278 
279  //Printf("%s",file.Data());
280 
281  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
282  if(!fxsec)
283  {
284  // next trial fetch the histgram file
285  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
286  if(!fxsec)
287  {
288  // not a severe condition but inciate that we have no information
289  return kFALSE;
290  }
291  else
292  {
293  // find the tlist we want to be independtent of the name so use the Tkey
294  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
295  if(!key)
296  {
297  fxsec->Close();
298  return kFALSE;
299  }
300 
301  TList *list = dynamic_cast<TList*>(key->ReadObj());
302  if(!list)
303  {
304  fxsec->Close();
305  return kFALSE;
306  }
307 
308  xsec = ((TProfile*)list->FindObject("h1Xsec")) ->GetBinContent(1);
309  trials = ((TH1F*) list->FindObject("h1Trials"))->GetBinContent(1);
310  fxsec->Close();
311  }
312  } // no tree pyxsec.root
313  else
314  {
315  TTree *xtree = (TTree*)fxsec->Get("Xsection");
316  if(!xtree)
317  {
318  fxsec->Close();
319  return kFALSE;
320  }
321 
322  UInt_t ntrials = 0;
323  Double_t xsection = 0;
324  xtree->SetBranchAddress("xsection",&xsection);
325  xtree->SetBranchAddress("ntrials",&ntrials);
326  xtree->GetEntry(0);
327  trials = ntrials;
328  xsec = xsection;
329  fxsec->Close();
330  }
331 
332  return kTRUE;
333 }
334 
335 //_______________________________________________________________________________________
341 //_______________________________________________________________________________________
343 {
344  if ( fhCentralityWeight ) delete fhCentralityWeight ;
345 
346  fhCentralityWeight = new TH1F("hCentralityWeights","Centrality weights",nbins,minCen,maxCen);
347 
348  for(Int_t ibin = 0; ibin < fhCentralityWeight->GetNbinsX(); ibin++)
349  fhCentralityWeight->SetBinContent(ibin,1.) ;
350 }
351 
352 //_________________________________________________
354 //_________________________________________________
356 {
358 
359  return fhCentralityWeight ;
360 }
361 
Int_t pdg
TList * GetCreateOutputHistograms()
double Double_t
Definition: External.C:58
TSystem * gSystem
Calculate the weight to the event to be applied when filling histograms.
Definition: AliAnaWeights.h:34
TH1F * fhCentralityWeight
Container of centrality weights.
Double_t GetParticlePtWeight(Float_t pt, Int_t pdg, TString genName, Int_t igen) const
Int_t fDebug
Debug level.
Float_t fCentrality
Container of centrality percentile.
Bool_t fJustFillCrossSecHist
Do not provide a weight, just fill cross section histograms.
virtual ~AliAnaWeights()
Destructor.
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
void InitCentralityWeightsHistogram(Int_t nbins=100, Int_t minCen=0, Int_t maxCen=100)
TF1 * fPi0Function
! pi0 spectrum parametrization
TH1F * GetCentralityWeightsHistogram()
AliAnaWeights()
Constructor.
TH1F * fhXsec
! Cross section in PYTHIA.
TF1 * fEtaFunction
! eta spectrum parametrization
Double_t fMCWeight
pT-hard bin MC weight. It is used only internally.
TFile * file
TList with histograms for a given trigger.
AliGenPythiaEventHeader * fPyEventHeader
! Pythia event header, needed to retrieve cross section, only in recent MC
TString fCurrFileName
Current file path name.
Bool_t fDoMCParticlePtWeights
activate the generation of a pT weight depending on MC particle pdg and generator ...
virtual Double_t GetWeight()
const Int_t nbins
bool Bool_t
Definition: External.C:53
Bool_t fUseCentralityWeight
Return the centratlity weight.
Bool_t fCheckPythiaEventHeader
Get cross section from pythia event header.
Bool_t fCheckMCCrossSection
Retrieve from the pyxsec.root file the cross section, only if requested.
virtual Double_t GetPythiaCrossSection()
static Bool_t GetPythiaInfoFromFile(TString currFile, Float_t &xsec, Float_t &trials)
TH1F * fhTrials
! Number of event trials in PYTHIA.