AliPhysics  9fe175b (9fe175b)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
35 
36 //______________________________________________________________
38 //______________________________________________________________
40 : TObject(), fDebug(0),
41  fhCentralityWeight(0),
42  fCentrality(0),
43  fUseCentralityWeight(0),
44  fCurrFileName(0),
45  fCheckMCCrossSection(kFALSE),
46  fJustFillCrossSecHist(0),
47  fhXsec(0),
48  fhTrials(0),
49  fPyEventHeader(0),
50  fCheckPythiaEventHeader(0)
51 {
52 }
53 
54 //____________________________________________________
57 //____________________________________________________
59 {
60  TList * outputContainer = new TList() ;
61  outputContainer->SetName("MCWeightHistogram") ;
62 
63  if(!fCheckMCCrossSection) return outputContainer;
64 
65  outputContainer->SetOwner(kFALSE);
66 
67  fhXsec = new TH1F("hXsec","xsec from pyxsec.root",1,0,1);
68  fhXsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
69  outputContainer->Add(fhXsec);
70 
71  fhTrials = new TH1F("hTrials","trials root file",1,0,1);
72  fhTrials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
73  outputContainer->Add(fhTrials);
74 
75  return outputContainer ;
76 }
77 
78 //__________________________________
83 //_________________________________
85 {
86  Double_t weight = 1.;
87 
89  {
90  Double_t temp = GetPythiaCrossSection() ;
91 
92  AliDebug(1,Form("MC pT-hard weight: %e",temp));
93 
94  weight*=temp;
95  }
96 
98  {
99  Double_t temp = fhCentralityWeight->GetBinContent(fhCentralityWeight->FindBin(fCentrality));
100 
101  AliDebug(1,Form("Centrality %2.1f, weight: %2.2f",fCentrality,temp));
102 
103  weight*=temp;
104  }
105 
106  AliDebug(1,Form("Event weight %e",weight));
107 
108  return weight ;
109 }
110 
111 //_____________________________________________
112 //
121 //_____________________________________________
123 {
124  Float_t xsection = 0;
125  Float_t trials = 1;
126  Float_t avgTrials = 0;
127 
128  if ( !fhXsec || !fhTrials ) return 1;
129 
130  // -----------------------------------------------
131  // Check if info is already in Pythia event header
132  // Do not apply the weight per event, too much number of trial variation
133  // -----------------------------------------------
135  {
136  fMCWeight = 1 ;
137 
138  AliDebug(fDebug,Form("Pythia event header: xs %2.2e, trial %d", fPyEventHeader->GetXsection(),fPyEventHeader->Trials()));
139 
140  fhXsec ->Fill("<#sigma>" ,fPyEventHeader->GetXsection());
141  fhTrials->Fill("#sum{ntrials}",fPyEventHeader->Trials());
142 
143  return 1 ;
144  }
145 
146  // -----------------------------------------------
147  // Get cross section from corresponding files
148  // -----------------------------------------------
149 
150  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
151  if ( !tree ) return 1;
152 
153  TFile *curfile = tree->GetCurrentFile();
154 
155  if ( !curfile ) return 1;
156 
157  // Check if file not accessed previously, if so
158  // return the previously calculated weight
159  if(fCurrFileName == curfile->GetName()) return fMCWeight;
160 
161  fCurrFileName = TString(curfile->GetName());
162 
163  Bool_t ok = GetPythiaInfoFromFile(fCurrFileName,xsection,trials);
164 
165  if ( !ok )
166  {
167  AliWarning("Parameters from file not recovered properly");
168  return 1;
169  }
170 
171  fhXsec->Fill("<#sigma>",xsection);
172 
173  // average number of trials
174  Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
175 
176  if(trials >= nEntries && nEntries > 0.) avgTrials = trials/nEntries;
177 
178  fhTrials->Fill("#sum{ntrials}",avgTrials);
179 
180  AliInfo(Form("xs %2.2e, trial %e, avg trials %2.2e, events per file %e",
181  xsection,trials,avgTrials,nEntries));
182 
183  AliDebug(1,Form("Reading File %s",curfile->GetName()));
184 
185  if(avgTrials > 0.)
186  {
188  {
189  fMCWeight = xsection / avgTrials ;
190 
191  AliInfo(Form("MC Weight: %e",fMCWeight));
192  }
193  else fMCWeight = 1; // do not add weight to histograms
194  }
195  else
196  {
197  AliWarning(Form("Average number of trials is NULL!! Set weight to 1: xs : %e, trials %e, entries %e",
198  xsection,trials,nEntries));
199 
200  fMCWeight = 1;
201  }
202 
203  return fMCWeight ;
204 }
205 
206 //_______________________________________________________________________________________
213 //_______________________________________________________________________________________
214 Bool_t AliAnaWeights::GetPythiaInfoFromFile(TString file,Float_t & xsec,Float_t & trials)
215 {
216  xsec = 0;
217  trials = 1;
218 
219  if(file.Contains("root_archive.zip#"))
220  {
221  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
222  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
223  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
224  file.Replace(pos+1,pos2-pos1,"");
225  }
226  else
227  {
228  // not an archive take the basename....
229  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
230  }
231 
232  //Printf("%s",file.Data());
233 
234  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
235  if(!fxsec)
236  {
237  // next trial fetch the histgram file
238  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
239  if(!fxsec)
240  {
241  // not a severe condition but inciate that we have no information
242  return kFALSE;
243  }
244  else
245  {
246  // find the tlist we want to be independtent of the name so use the Tkey
247  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
248  if(!key)
249  {
250  fxsec->Close();
251  return kFALSE;
252  }
253 
254  TList *list = dynamic_cast<TList*>(key->ReadObj());
255  if(!list)
256  {
257  fxsec->Close();
258  return kFALSE;
259  }
260 
261  xsec = ((TProfile*)list->FindObject("h1Xsec")) ->GetBinContent(1);
262  trials = ((TH1F*) list->FindObject("h1Trials"))->GetBinContent(1);
263  fxsec->Close();
264  }
265  } // no tree pyxsec.root
266  else
267  {
268  TTree *xtree = (TTree*)fxsec->Get("Xsection");
269  if(!xtree)
270  {
271  fxsec->Close();
272  return kFALSE;
273  }
274 
275  UInt_t ntrials = 0;
276  Double_t xsection = 0;
277  xtree->SetBranchAddress("xsection",&xsection);
278  xtree->SetBranchAddress("ntrials",&ntrials);
279  xtree->GetEntry(0);
280  trials = ntrials;
281  xsec = xsection;
282  fxsec->Close();
283  }
284 
285  return kTRUE;
286 }
287 
288 //_______________________________________________________________________________________
294 //_______________________________________________________________________________________
295 void AliAnaWeights::InitCentralityWeightsHistogram(Int_t nbins, Int_t minCen, Int_t maxCen)
296 {
297  if ( fhCentralityWeight ) delete fhCentralityWeight ;
298 
299  fhCentralityWeight = new TH1F("hCentralityWeights","Centrality weights",nbins,minCen,maxCen);
300 
301  for(Int_t ibin = 0; ibin < fhCentralityWeight->GetNbinsX(); ibin++)
302  fhCentralityWeight->SetBinContent(ibin,1.) ;
303 }
304 
305 //_________________________________________________
307 //_________________________________________________
309 {
311 
312  return fhCentralityWeight ;
313 }
314 
TList * GetCreateOutputHistograms()
TSystem * gSystem
Calculate the weight to the event to be applied when filling histograms.
Definition: AliAnaWeights.h:32
TList * list
TH1F * fhCentralityWeight
Container of centrality weights.
Int_t fDebug
Debug level.
Definition: AliAnaWeights.h:98
Float_t fCentrality
Container of centrality percentile.
Bool_t fJustFillCrossSecHist
Do not provide a weight, just fill cross section histograms.
void InitCentralityWeightsHistogram(Int_t nbins=100, Int_t minCen=0, Int_t maxCen=100)
TH1F * GetCentralityWeightsHistogram()
AliAnaWeights()
Constructor.
TH1F * fhXsec
! Cross section in PYTHIA.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Double_t fMCWeight
pT-hard bin MC weight. It is used only internally.
TFile * file
AliGenPythiaEventHeader * fPyEventHeader
! Pythia event header, needed to retrieve cross section, only in recent MC
TString fCurrFileName
Current file path name.
virtual Double_t GetWeight()
const Int_t nbins
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.