AliPhysics  608b256 (608b256)
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  weight*=temp;
109  }
110 
111  if ( fUseCentralityWeight )
112  {
113  Double_t temp = fhCentralityWeight->GetBinContent(fhCentralityWeight->FindBin(fCentrality));
114 
115  AliDebug(1,Form("Centrality %2.1f, weight: %2.2f",fCentrality,temp));
116 
117  weight*=temp;
118  }
119 
120  AliDebug(1,Form("Event weight %e",weight));
121 
122  return weight ;
123 }
124 
125 //_____________________________________________
126 //
138 //_____________________________________________
140 {
141  Double_t weight = 1.;
142 
143  if ( !fDoMCParticlePtWeights ) return weight ;
144 
145  if (pdg == 111 && fPi0Function &&
146  (genName.Contains("Pi0") || genName.Contains("pi0") || genName.Contains("PARAM") || genName.Contains("BOX")))
147  weight = fPi0Function->Eval(pt);
148  else if (pdg == 221 && fEtaFunction &&
149  (genName.Contains("Eta") || genName.Contains("eta") || genName.Contains("PARAM") || genName.Contains("BOX")))
150  weight = fEtaFunction->Eval(pt);
151 
152  AliDebug(1,Form("MC particle pdg %d, pt %2.2f, generator %s with index %d: weight %f",pdg,pt,genName.Data(),igen, weight));
153 
154  return weight ;
155 }
156 
157 //_____________________________________________
158 //
167 //_____________________________________________
169 {
170  Float_t xsection = 0;
171  Float_t trials = 1;
172  Float_t avgTrials = 0;
173 
174  if ( !fhXsec || !fhTrials ) return 1;
175 
176  // -----------------------------------------------
177  // Check if info is already in Pythia event header
178  // Do not apply the weight per event, too much number of trial variation
179  // -----------------------------------------------
181  {
182  AliDebug(fDebug,Form("Pythia event header: xs %2.2e, trial %d",
183  fPyEventHeader->GetXsection(),fPyEventHeader->Trials()));
184 
185  fhXsec ->Fill("<#sigma>" ,fPyEventHeader->GetXsection());
186  fhTrials->Fill("#sum{ntrials}",fPyEventHeader->Trials());
187 
188  if ( !fJustFillCrossSecHist )
189  {
190  fMCWeight = fPyEventHeader->GetXsection() / fPyEventHeader->Trials() ;
191  AliDebug(1,Form("MC Weight: %e",fMCWeight));
192  }
193  else fMCWeight = 1;
194 
195  return fMCWeight ;
196  }
197 
198  // -----------------------------------------------
199  // Get cross section from corresponding files
200  // -----------------------------------------------
201 
202  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
203  if ( !tree ) return 1;
204 
205  TFile *curfile = tree->GetCurrentFile();
206 
207  if ( !curfile ) return 1;
208 
209  // Check if file not accessed previously, if so
210  // return the previously calculated weight
211  if(fCurrFileName == curfile->GetName()) return fMCWeight;
212 
213  fCurrFileName = TString(curfile->GetName());
214 
215  Bool_t ok = GetPythiaInfoFromFile(fCurrFileName,xsection,trials);
216 
217  if ( !ok )
218  {
219  AliWarning("Parameters from file not recovered properly");
220  return 1;
221  }
222 
223  fhXsec->Fill("<#sigma>",xsection);
224 
225  // average number of trials
226  Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
227 
228  if(trials >= nEntries && nEntries > 0.) avgTrials = trials/nEntries;
229 
230  fhTrials->Fill("#sum{ntrials}",avgTrials);
231 
232  AliInfo(Form("xs %2.2e, trial %e, avg trials %2.2e, events per file %e",
233  xsection,trials,avgTrials,nEntries));
234 
235  AliDebug(1,Form("Reading File %s",curfile->GetName()));
236 
237  if(avgTrials > 0.)
238  {
240  {
241  fMCWeight = xsection / avgTrials ;
242 
243  AliInfo(Form("MC Weight: %e",fMCWeight));
244  }
245  else fMCWeight = 1; // do not add weight to histograms
246  }
247  else
248  {
249  AliWarning(Form("Average number of trials is NULL!! Set weight to 1: xs : %e, trials %e, entries %e",
250  xsection,trials,nEntries));
251 
252  fMCWeight = 1;
253  }
254 
255  return fMCWeight ;
256 }
257 
258 //_______________________________________________________________________________________
265 //_______________________________________________________________________________________
267 {
268  xsec = 0;
269  trials = 1;
270 
271  if(file.Contains("root_archive.zip#"))
272  {
273  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
274  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
275  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
276  file.Replace(pos+1,pos2-pos1,"");
277  }
278  else
279  {
280  // not an archive take the basename....
281  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
282  }
283 
284  //Printf("%s",file.Data());
285 
286  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
287  if(!fxsec)
288  {
289  // next trial fetch the histgram file
290  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
291  if(!fxsec)
292  {
293  // not a severe condition but inciate that we have no information
294  return kFALSE;
295  }
296  else
297  {
298  // find the tlist we want to be independtent of the name so use the Tkey
299  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
300  if(!key)
301  {
302  fxsec->Close();
303  return kFALSE;
304  }
305 
306  TList *list = dynamic_cast<TList*>(key->ReadObj());
307  if(!list)
308  {
309  fxsec->Close();
310  return kFALSE;
311  }
312 
313  xsec = ((TProfile*)list->FindObject("h1Xsec")) ->GetBinContent(1);
314  trials = ((TH1F*) list->FindObject("h1Trials"))->GetBinContent(1);
315  fxsec->Close();
316  }
317  } // no tree pyxsec.root
318  else
319  {
320  TTree *xtree = (TTree*)fxsec->Get("Xsection");
321  if(!xtree)
322  {
323  fxsec->Close();
324  return kFALSE;
325  }
326 
327  UInt_t ntrials = 0;
328  Double_t xsection = 0;
329  xtree->SetBranchAddress("xsection",&xsection);
330  xtree->SetBranchAddress("ntrials",&ntrials);
331  xtree->GetEntry(0);
332  trials = ntrials;
333  xsec = xsection;
334  fxsec->Close();
335  }
336 
337  return kTRUE;
338 }
339 
340 //_______________________________________________________________________________________
346 //_______________________________________________________________________________________
348 {
349  if ( fhCentralityWeight ) delete fhCentralityWeight ;
350 
351  fhCentralityWeight = new TH1F("hCentralityWeights","Centrality weights",nbins,minCen,maxCen);
352 
353  for(Int_t ibin = 0; ibin < fhCentralityWeight->GetNbinsX(); ibin++)
354  fhCentralityWeight->SetBinContent(ibin,1.) ;
355 }
356 
357 //_________________________________________________
359 //_________________________________________________
361 {
363 
364  return fhCentralityWeight ;
365 }
366 
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.