]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGJE/AliPWG4HighPtQAMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGJE / AliPWG4HighPtQAMC.cxx
... / ...
CommitLineData
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//-----------------------------------------------------------------------
17// This class compares the global reconstruction with the MC information
18// Momentum resolution is stored as function of track cuts and pt.
19// Output: Histograms for different set of cuts
20//-----------------------------------------------------------------------
21// Author : Marta Verweij - UU
22//-----------------------------------------------------------------------
23
24#ifndef ALIPWG4HighPtQAMC_CXX
25#define ALIPWG4HighPtQAMC_CXX
26
27#include "AliPWG4HighPtQAMC.h"
28
29#include "TH1.h"
30#include "TH2.h"
31#include "TH3.h"
32#include "TProfile.h"
33#include "TList.h"
34#include "TFile.h"
35#include "TChain.h"
36#include "TH3F.h"
37#include "TKey.h"
38#include "TSystem.h"
39
40#include "AliAnalysisTask.h"
41#include "AliAnalysisManager.h"
42#include "AliESDInputHandler.h"
43#include "AliMCEvent.h"
44#include "AliMCEventHandler.h"
45#include "AliStack.h"
46#include "AliESDtrack.h"
47#include "AliESDtrackCuts.h"
48#include "AliExternalTrackParam.h"
49#include "AliLog.h"
50#include "AliGenPythiaEventHeader.h"
51#include "AliGenCocktailEventHeader.h"
52//#include "AliAnalysisHelperJetTasks.h"
53
54using namespace std; //required for resolving the 'cout' symbol
55
56ClassImp(AliPWG4HighPtQAMC)
57
58AliPWG4HighPtQAMC::AliPWG4HighPtQAMC()
59: AliAnalysisTask("AliPWG4HighPtQAMC", ""),
60 fESD(0),
61 fMC(0),
62 fStack(0),
63 fVtx(0x0),
64 fTrackCuts(0),
65 fTrackCutsReject(),
66 fTrackType(0),
67 fSigmaConstrainedMax(1e6),
68 fPtMax(100.),
69 fAvgTrials(1),
70 fNEventAll(0),
71 fNEventSel(0),
72 fNEventReject(0),
73 fh1Xsec(0),
74 fh1Trials(0),
75 fh1PtHard(0),
76 fh1PtHardTrials(0),
77 fPtAll(0),
78 fPtSel(0),
79 fPtSelFakes(0),
80 fNPointTPCFakes(0),
81 fPtSelLargeLabel(0),
82 fMultRec(0),
83 fNPointTPCMultRec(0),
84 fDeltaPtMultRec(0),
85 fPtAllvsPtMC(0),
86 fPtAllminPtMCvsPtMC(0),
87 fPtAllminPtMCvsPtAll(0),
88 fPtAllvsPtMCvsMult(0),
89 fPtAllminPtMCvsPtAllvsMult(0),
90 fPtAllminPtMCvsPtAllNPointTPC(0),
91 fPtAllminPtMCvsPtAllNPointTPCIter1(0),
92 fPtAllminPtMCvsPtAllChi2TPC(0),
93 fPtAllminPtMCvsPtAllChi2TPCIter1(0),
94 fPtAllminPtMCvsPtAllDCAR(0),
95 fPtAllminPtMCvsPtAllDCAZ(0),
96 fPtAllminPtMCvsPtAllPhi(0),
97 fPtAllminPtMCvsPtAllNPointITS(0),
98 fPtAllminPtMCvsPtAllNSigmaToVertex(0),
99 fPtAllminPtMCvsPtAllChi2C(0),
100 fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
101 fPtAllMC(0),
102 fPtSelMC(0),
103 fHistList(0)
104{
105
106 fPtBinEdges[0][0] = 10.;
107 fPtBinEdges[0][1] = 1.;
108 fPtBinEdges[1][0] = 20.;
109 fPtBinEdges[1][1] = 2.;
110 fPtBinEdges[2][0] = 100.;
111 fPtBinEdges[2][1] = 10.;
112
113}
114//________________________________________________________________________
115AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
116 AliAnalysisTask(name,""),
117 fESD(0),
118 fMC(0),
119 fStack(0),
120 fVtx(0x0),
121 fTrackCuts(),
122 fTrackCutsReject(),
123 fTrackType(0),
124 fSigmaConstrainedMax(1e6),
125 fPtMax(100.),
126 fAvgTrials(1),
127 fNEventAll(0),
128 fNEventSel(0),
129 fNEventReject(0),
130 fh1Xsec(0),
131 fh1Trials(0),
132 fh1PtHard(0),
133 fh1PtHardTrials(0),
134 fPtAll(0),
135 fPtSel(0),
136 fPtSelFakes(0),
137 fNPointTPCFakes(0),
138 fPtSelLargeLabel(0),
139 fMultRec(0),
140 fNPointTPCMultRec(0),
141 fDeltaPtMultRec(0),
142 fPtAllvsPtMC(0),
143 fPtAllminPtMCvsPtMC(0),
144 fPtAllminPtMCvsPtAll(0),
145 fPtAllvsPtMCvsMult(0),
146 fPtAllminPtMCvsPtAllvsMult(0),
147 fPtAllminPtMCvsPtAllNPointTPC(0),
148 fPtAllminPtMCvsPtAllNPointTPCIter1(0),
149 fPtAllminPtMCvsPtAllChi2TPC(0),
150 fPtAllminPtMCvsPtAllChi2TPCIter1(0),
151 fPtAllminPtMCvsPtAllDCAR(0),
152 fPtAllminPtMCvsPtAllDCAZ(0),
153 fPtAllminPtMCvsPtAllPhi(0),
154 fPtAllminPtMCvsPtAllNPointITS(0),
155 fPtAllminPtMCvsPtAllNSigmaToVertex(0),
156 fPtAllminPtMCvsPtAllChi2C(0),
157 fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
158 fPtAllMC(0),
159 fPtSelMC(0),
160 fHistList(0)
161{
162 //
163 // Constructor. Initialization of Inputs and Outputs
164 //
165 AliDebug(2,Form("AliPWG4HighPtQAMC Calling Constructor"));
166
167 fPtBinEdges[0][0] = 10.;
168 fPtBinEdges[0][1] = 1.;
169 fPtBinEdges[1][0] = 20.;
170 fPtBinEdges[1][1] = 2.;
171 fPtBinEdges[2][0] = 100.;
172 fPtBinEdges[2][1] = 10.;
173
174 // Input slot #0 works with a TChain ESD
175 DefineInput(0, TChain::Class());
176 // Output slot #0, #1 write into a TList
177 DefineOutput(0, TList::Class());
178}
179
180//________________________________________________________________________
181void AliPWG4HighPtQAMC::SetPtBinEdges(Int_t region, Double_t ptmax, Double_t ptBinWidth) {
182
183 if(region<3) {
184 fPtBinEdges[region][0] = ptmax;
185 fPtBinEdges[region][1] = ptBinWidth;
186 }
187 else {
188 AliError("Only 3 regions alowed. Use region 0/1/2\n");
189 return;
190 }
191
192}
193
194//________________________________________________________________________
195void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
196{
197 // Connect ESD and MC event handler here
198 // Called once
199 AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
200 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
201 if (!tree) {
202 AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
203 return;
204 }
205
206 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
207
208 if (!esdH) {
209 AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
210 return;
211 } else
212 fESD = esdH->GetEvent();
213
214 AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
215 if (!eventHandler) {
216 AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
217 }
218 else
219 fMC = eventHandler->MCEvent();
220
221}
222
223
224//________________________________________________________________________
225void AliPWG4HighPtQAMC::CreateOutputObjects() {
226 //Create output objects
227 AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
228
229 Bool_t oldStatus = TH1::AddDirectoryStatus();
230 TH1::AddDirectory(kFALSE);
231
232 OpenFile(0);
233 fHistList = new TList();
234 fHistList->SetOwner(kTRUE);
235
236 Float_t fgkPtMin = 0.;
237
238 //fPtBinEdges[region][0] = ptmax of region ; fPtBinEdges[region][1] = binWidth of region
239 const Float_t ptmin1 = fgkPtMin;
240 const Float_t ptmax1 = fPtBinEdges[0][0];
241 const Float_t ptmin2 = ptmax1 ;
242 const Float_t ptmax2 = fPtBinEdges[1][0];
243 const Float_t ptmin3 = ptmax2 ;
244 const Float_t ptmax3 = fPtBinEdges[2][0];//fgkPtMax;
245 const Int_t nbin11 = (int)((ptmax1-ptmin1)/fPtBinEdges[0][1]);
246 const Int_t nbin12 = (int)((ptmax2-ptmin2)/fPtBinEdges[1][1])+nbin11;
247 const Int_t nbin13 = (int)((ptmax3-ptmin3)/fPtBinEdges[2][1])+nbin12;
248 Int_t fgkNPtBins=nbin13;
249 //Create array with low edges of each bin
250 Double_t *binsPt=new Double_t[fgkNPtBins+1];
251 for(Int_t i=0; i<=fgkNPtBins; i++) {
252 if(i<=nbin11) binsPt[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
253 if(i<=nbin12 && i>nbin11) binsPt[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
254 if(i<=nbin13 && i>nbin12) binsPt[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
255 }
256
257 Int_t fgkNPhiBins = 18*6;
258 Float_t kMinPhi = 0.;
259 Float_t kMaxPhi = 2.*TMath::Pi();
260 Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
261 for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
262
263 Int_t fgkResPtBins=80;
264 Float_t kMinResPt = -1.;
265 Float_t kMaxResPt = 1.;
266 Double_t *binsResPt = new Double_t[fgkResPtBins+1];
267 for(Int_t i=0; i<=fgkResPtBins; i++) binsResPt[i]=(Double_t)kMinResPt + (kMaxResPt-kMinResPt)/fgkResPtBins*(Double_t)i ;
268
269 Int_t fgkMultBins=50;
270 Float_t kMinMult = 0.;
271 Float_t kMaxMult = 3000.;
272 Double_t *binsMult = new Double_t[fgkMultBins+1];
273 for(Int_t i=0; i<=fgkMultBins; i++) binsMult[i]=(Double_t)kMinMult + (kMaxMult-kMinMult)/fgkMultBins*(Double_t)i ;
274
275 Int_t fgkNEtaBins=20;
276 Float_t fgkEtaMin = -1.;
277 Float_t fgkEtaMax = 1.;
278 Double_t *binsEta=new Double_t[fgkNEtaBins+1];
279 for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
280
281 Int_t fgkNNClustersTPCBins=80;
282 Float_t fgkNClustersTPCMin = 0.5;
283 Float_t fgkNClustersTPCMax = 160.5;
284 Double_t *binsNClustersTPC=new Double_t[fgkNNClustersTPCBins+1];
285 for(Int_t i=0; i<=fgkNNClustersTPCBins; i++) binsNClustersTPC[i]=(Double_t)fgkNClustersTPCMin + (fgkNClustersTPCMax-fgkNClustersTPCMin)/fgkNNClustersTPCBins*(Double_t)i ;
286
287 Int_t fgkNDCA2DBins=80;
288 Float_t fgkDCA2DMin = -0.2;
289 Float_t fgkDCA2DMax = 0.2;
290 if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
291 fgkDCA2DMin = -2.;
292 fgkDCA2DMax = 2.;
293 }
294 Double_t *binsDCA2D=new Double_t[fgkNDCA2DBins+1];
295 for(Int_t i=0; i<=fgkNDCA2DBins; i++) binsDCA2D[i]=(Double_t)fgkDCA2DMin + (fgkDCA2DMax-fgkDCA2DMin)/fgkNDCA2DBins*(Double_t)i ;
296
297 Int_t fgkNDCAZBins=80;
298 Float_t fgkDCAZMin = -2.;
299 Float_t fgkDCAZMax = 2.;
300 if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
301 fgkDCAZMin = -5.;
302 fgkDCAZMax = 5.;
303 }
304 Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
305 for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
306
307 Int_t fgkNNPointITSBins=9;
308 Float_t fgkNPointITSMin = -0.5;
309 Float_t fgkNPointITSMax = 8.5;
310 Double_t *binsNPointITS=new Double_t[fgkNNPointITSBins+1];
311 for(Int_t i=0; i<=fgkNNPointITSBins; i++) binsNPointITS[i]=(Double_t)fgkNPointITSMin + (fgkNPointITSMax-fgkNPointITSMin)/fgkNNPointITSBins*(Double_t)i ;
312
313 Int_t fgkNNSigmaToVertexBins=20;
314 Float_t fgkNSigmaToVertexMin = 0.;
315 Float_t fgkNSigmaToVertexMax = 8.;
316 Double_t *binsNSigmaToVertex=new Double_t[fgkNNSigmaToVertexBins+1];
317 for(Int_t i=0; i<=fgkNNSigmaToVertexBins; i++) binsNSigmaToVertex[i]=(Double_t)fgkNSigmaToVertexMin + (fgkNSigmaToVertexMax-fgkNSigmaToVertexMin)/fgkNNSigmaToVertexBins*(Double_t)i ;
318
319 Int_t fgkNChi2CBins=20;
320 Float_t fgkChi2CMin = 0.;
321 Float_t fgkChi2CMax = 100.;
322 Double_t *binsChi2C=new Double_t[fgkNChi2CBins+1];
323 for(Int_t i=0; i<=fgkNChi2CBins; i++) binsChi2C[i]=(Double_t)fgkChi2CMin + (fgkChi2CMax-fgkChi2CMin)/fgkNChi2CBins*(Double_t)i ;
324
325 Int_t fgkNRel1PtUncertaintyBins=50;
326 Float_t fgkRel1PtUncertaintyMin = 0.;
327 Float_t fgkRel1PtUncertaintyMax = 1.;
328 Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
329 for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
330
331 Float_t fgkChi2PerClusMin = 0.;
332 Float_t fgkChi2PerClusMax = 4.;
333 Int_t fgkNChi2PerClusBins = (int)(fgkChi2PerClusMax*10.);
334 Double_t *binsChi2PerClus=new Double_t[fgkNChi2PerClusBins+1];
335 for(Int_t i=0; i<=fgkNChi2PerClusBins; i++) binsChi2PerClus[i]=(Double_t)fgkChi2PerClusMin + (fgkChi2PerClusMax-fgkChi2PerClusMin)/fgkNChi2PerClusBins*(Double_t)i ;
336
337
338 fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
339 fHistList->Add(fNEventAll);
340 fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
341 fHistList->Add(fNEventSel);
342 fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
343 //Set labels
344 fNEventReject->Fill("noESD",0);
345 fNEventReject->Fill("Trigger",0);
346 fNEventReject->Fill("noMCEvent",0);
347 fNEventReject->Fill("noStack",0);
348 fNEventReject->Fill("NTracks<2",0);
349 fNEventReject->Fill("noVTX",0);
350 fNEventReject->Fill("VtxStatus",0);
351 fNEventReject->Fill("NCont<2",0);
352 fNEventReject->Fill("ZVTX>10",0);
353 fHistList->Add(fNEventReject);
354
355 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
356 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
357 fHistList->Add(fh1Xsec);
358
359 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
360 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
361 fHistList->Add(fh1Trials);
362
363 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
364 fHistList->Add(fh1PtHard);
365 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
366 fHistList->Add(fh1PtHardTrials);
367
368 fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
369 fHistList->Add(fPtAll);
370 fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
371 fHistList->Add(fPtSel);
372
373 fPtSelFakes = new TH1F("fPtSelFakes","PtSelFakes",fgkNPtBins, binsPt);
374 fHistList->Add(fPtSelFakes);
375
376 fNPointTPCFakes = new TH1F("fNPointTPCFakes","fNPointTPCFakes",fgkNNClustersTPCBins, binsNClustersTPC);
377 fHistList->Add(fNPointTPCFakes);
378
379 fPtSelLargeLabel = new TH1F("fPtSelLargeLabel","PtSelLargeLabel",fgkNPtBins, binsPt);
380 fHistList->Add(fPtSelLargeLabel);
381
382 fMultRec = new TH1F("fMultRec","Multiple reconstruction of tracks",fgkMultBins, binsMult);
383 fHistList->Add(fMultRec);
384
385 fNPointTPCMultRec = new TH1F("fNPointTPCMultRec","Multiple reconstruction of tracks",fgkNNClustersTPCBins, binsNClustersTPC);
386 fHistList->Add(fNPointTPCMultRec);
387
388 fDeltaPtMultRec = new TH2F("fDeltaPtMultRec","Delta pT vs pT for multiple reconstructed tracks",50,0.,50.,40,-20.,20.);
389 fHistList->Add(fDeltaPtMultRec);
390
391 fPtAllvsPtMC = new TH2F("fPtAllvsPtMC","fPtAllvsPtMC;p_{T,MC};p_{T,rec}",fgkNPtBins, binsPt,fgkNPtBins, binsPt);
392 fHistList->Add(fPtAllvsPtMC);
393
394 fPtAllminPtMCvsPtMC = new TH2F("fPtAllminPtMCvsPtMC","PtAllminPtMCvsPtMC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt);
395 fPtAllminPtMCvsPtMC->SetXTitle("p_{t}^{MC}");
396 fPtAllminPtMCvsPtMC->SetYTitle("(1/p_{t}^{rec}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
397 fHistList->Add(fPtAllminPtMCvsPtMC);
398
399 fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, binsPt,fgkResPtBins,binsResPt);
400 fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{rec}");
401 fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{rec}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
402 fHistList->Add(fPtAllminPtMCvsPtAll);
403
404 fPtAllvsPtMCvsMult = new TH3F("fPtAllvsPtMCvsMult","fPtAllvsPtMCvsMult;p_{T,MC};p_{T,rec};N_{tracks}",fgkNPtBins, binsPt,fgkNPtBins, binsPt,fgkMultBins,binsMult);
405 fHistList->Add(fPtAllvsPtMCvsMult);
406
407 fPtAllminPtMCvsPtAllvsMult = new TH3F("fPtAllminPtMCvsPtAllvsMult","fPtAllminPtMCvsPtAllvsMult",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkMultBins,binsMult);
408 fPtAllminPtMCvsPtAllvsMult->SetXTitle("p_{t}^{MC}");
409 fPtAllminPtMCvsPtAllvsMult->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
410 fPtAllminPtMCvsPtAllvsMult->SetZTitle("N_{tracks}");
411 fHistList->Add(fPtAllminPtMCvsPtAllvsMult);
412
413 fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNClustersTPCBins, binsNClustersTPC);
414 fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
415 fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
416 fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{cls,TPC}");
417 fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
418
419 fPtAllminPtMCvsPtAllNPointTPCIter1 = new TH3F("fPtAllminPtMCvsPtAllNPointTPCIter1","PtAllminPtMCvsPtAllNPointTPCIter1",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNClustersTPCBins, binsNClustersTPC);
420 fPtAllminPtMCvsPtAllNPointTPCIter1->SetXTitle("p_{t}^{MC}");
421 fPtAllminPtMCvsPtAllNPointTPCIter1->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
422 fPtAllminPtMCvsPtAllNPointTPCIter1->SetZTitle("N_{clsIter1,TPC}");
423 fHistList->Add(fPtAllminPtMCvsPtAllNPointTPCIter1);
424
425 fPtAllminPtMCvsPtAllChi2TPC = new TH3F("fPtAllminPtMCvsPtAllChi2TPC","PtAllminPtMCvsPtAllChi2TPC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2PerClusBins, binsChi2PerClus);
426 fPtAllminPtMCvsPtAllChi2TPC->SetXTitle("p_{t}^{MC}");
427 fPtAllminPtMCvsPtAllChi2TPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
428 fPtAllminPtMCvsPtAllChi2TPC->SetZTitle("#chi^{2} TPC");
429 fHistList->Add(fPtAllminPtMCvsPtAllChi2TPC);
430
431 fPtAllminPtMCvsPtAllChi2TPCIter1 = new TH3F("fPtAllminPtMCvsPtAllChi2TPCIter1","PtAllminPtMCvsPtAllChi2TPCIter1",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2PerClusBins, binsChi2PerClus);
432 fPtAllminPtMCvsPtAllChi2TPCIter1->SetXTitle("p_{t}^{MC}");
433 fPtAllminPtMCvsPtAllChi2TPCIter1->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
434 fPtAllminPtMCvsPtAllChi2TPCIter1->SetZTitle("#chi^{2} TPC Iter1");
435 fHistList->Add(fPtAllminPtMCvsPtAllChi2TPCIter1);
436
437 fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNDCA2DBins,binsDCA2D);
438 fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
439 fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
440 fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
441 fHistList->Add(fPtAllminPtMCvsPtAllDCAR);
442
443 fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNDCAZBins,binsDCAZ);
444 fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
445 fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
446 fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
447 fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);
448
449 fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNPhiBins,binsPhi);
450 fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
451 fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
452 fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
453 fHistList->Add(fPtAllminPtMCvsPtAllPhi);
454
455 fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNPointITSBins,binsNPointITS);
456 fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
457 fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
458 fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
459 fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
460
461 fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
462 fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
463 fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
464 fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
465 fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);
466
467 fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2CBins,binsChi2C);
468 fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
469 fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
470 fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
471 fHistList->Add(fPtAllminPtMCvsPtAllChi2C);
472
473 fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
474 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
475 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
476 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
477 fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);
478
479
480 fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, binsPt);
481 fHistList->Add(fPtAllMC);
482 fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, binsPt);
483 fHistList->Add(fPtSelMC);
484
485 TH1::AddDirectory(oldStatus);
486
487 PostData(0, fHistList);
488
489 if(binsPt) delete [] binsPt;
490 if(binsResPt) delete [] binsResPt;
491 if(binsPhi) delete [] binsPhi;
492 if(binsNClustersTPC) delete [] binsNClustersTPC;
493 if(binsDCA2D) delete [] binsDCA2D;
494 if(binsDCAZ) delete [] binsDCAZ;
495 if(binsNPointITS) delete [] binsNPointITS;
496 if(binsNSigmaToVertex) delete [] binsNSigmaToVertex;
497 if(binsChi2C) delete [] binsChi2C;
498 if(binsEta) delete [] binsEta;
499 if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
500 if(binsChi2PerClus) delete [] binsChi2PerClus;
501 if(binsMult) delete [] binsMult;
502
503}
504
505//________________________________________________________________________
506Bool_t AliPWG4HighPtQAMC::SelectEvent() {
507 //
508 // Decide if event should be selected for analysis
509 //
510
511 // Checks following requirements:
512 // - fESD available
513 // - trigger info from AliPhysicsSelection
514 // - MCevent available
515 // - number of reconstructed tracks > 1
516 // - primary vertex reconstructed
517 // - z-vertex < 10 cm
518
519 Bool_t selectEvent = kTRUE;
520
521 //fESD object available?
522 if (!fESD) {
523 AliDebug(2,Form("ERROR: fInputEvent not available\n"));
524 fNEventReject->Fill("noESD",1);
525 selectEvent = kFALSE;
526 return selectEvent;
527 }
528
529 //Trigger
530 UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
531 if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
532 AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
533 fNEventReject->Fill("Trigger",1);
534 selectEvent = kFALSE;
535 return selectEvent;
536 }
537
538 //MCEvent available?
539 //if yes: get stack
540 if(fMC) {
541 AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
542 fStack = fMC->Stack(); //Particles Stack
543 if(fStack) {
544 AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
545 }
546 else {
547 AliDebug(2,Form("ERROR: Could not retrieve MC eventHandler"));
548 fNEventReject->Fill("noStack",1);
549 selectEvent = kFALSE;
550 return selectEvent;
551 }
552 } else {
553 AliDebug(2,Form("ERROR: Could not retrieve stack"));
554 fNEventReject->Fill("noMCEvent",1);
555 selectEvent = kFALSE;
556 return selectEvent;
557 }
558
559 //Check if number of reconstructed tracks is larger than 1
560 if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2) {
561 fNEventReject->Fill("NTracks<2",1);
562 selectEvent = kFALSE;
563 return selectEvent;
564 }
565
566 //Check if vertex is reconstructed
567 if(fTrackType==1) fVtx = fESD->GetPrimaryVertexTPC();
568 else fVtx = fESD->GetPrimaryVertexSPD();
569
570 if(!fVtx) {
571 fNEventReject->Fill("noVTX",1);
572 selectEvent = kFALSE;
573 return selectEvent;
574 }
575
576 if(!fVtx->GetStatus()) {
577 fNEventReject->Fill("VtxStatus",1);
578 selectEvent = kFALSE;
579 return selectEvent;
580 }
581
582 // Need vertex cut
583 // TString vtxName(fVtx->GetName());
584 if(fVtx->GetNContributors()<2) {
585 fNEventReject->Fill("NCont<2",1);
586 selectEvent = kFALSE;
587 return selectEvent;
588 }
589
590 //Check if z-vertex < 10 cm
591 double primVtx[3];
592 fVtx->GetXYZ(primVtx);
593 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
594 fNEventReject->Fill("ZVTX>10",1);
595 selectEvent = kFALSE;
596 return selectEvent;
597 }
598
599 AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",fVtx->GetTitle(), fVtx->GetStatus(), fVtx->GetNContributors()));
600
601 return selectEvent;
602
603}
604
605//________________________________________________________________________
606void AliPWG4HighPtQAMC::Exec(Option_t *) {
607 // Main loop
608 // Called for each event
609 AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
610 // All events without selection
611 fNEventAll->Fill(0.);
612
613 if(!SelectEvent()) {
614 // Post output data
615 PostData(0, fHistList);
616 return;
617 }
618
619 // ---- Get MC Header information (for MC productions in pThard bins) ----
620 Double_t ptHard = 0.;
621 Double_t nTrials = 1; // trials for MC trigger weight for real data
622
623 if(fMC){
624 AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(fMC);
625 if(pythiaGenHeader){
626 nTrials = pythiaGenHeader->Trials();
627 ptHard = pythiaGenHeader->GetPtHard();
628
629 fh1PtHard->Fill(ptHard);
630 fh1PtHardTrials->Fill(ptHard,nTrials);
631
632 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
633 }
634 }
635
636 //Need to keep track of selected events
637 fNEventSel->Fill(0.);
638
639 Int_t nTracks = fESD->GetNumberOfTracks();
640 AliDebug(2,Form("nTracks ESD%d", nTracks));
641
642 int nMCtracks = fStack->GetNtrack();
643
644 Float_t pt = 0.;
645 Float_t ptMC = 0.;
646 Float_t phi = 0.;
647 Float_t dca2D = 0.;
648 Float_t dcaZ = 0.;
649 Int_t nPointITS = 0;
650 Float_t chi2C = 0.;
651 Float_t nSigmaToVertex = 0.;
652 Float_t relUncertainty1Pt = 0.;
653
654 int mult = fTrackCuts->CountAcceptedTracks(fESD);
655
656 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
657
658 AliESDtrack *track = 0;
659 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
660 if(!esdtrack) continue;
661
662 if(fTrackType==1) {
663 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
664 if(!track) continue;
665 }
666 else if(fTrackType==2) {
667 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
668 if(!track) continue;
669
670 AliExternalTrackParam exParam;
671 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
672 if( !relate ) {
673 delete track;
674 continue;
675 }
676 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
677 }
678 else if(fTrackType==7) {
679 //use global constrained track
680 track = new AliESDtrack(*esdtrack);
681 // track = esdtrack;
682 // track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
683 }
684 else
685 track = esdtrack;
686
687 if(!track) continue;
688
689 if(fTrackType==2) {
690 //Cut on chi2 of constrained fit
691 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
692 delete track;
693 continue;
694 }
695 }
696
697 Int_t label = TMath::Abs(track->GetLabel());
698 if(label>=nMCtracks) {
699 if (fTrackCuts->AcceptTrack(track)) { fPtSelLargeLabel->Fill(pt); }
700 if(fTrackType==1 || fTrackType==2 || fTrackType==7) delete track;
701 continue;
702 }
703
704 TParticle *particle = fStack->Particle(label) ;
705 if(!particle) {
706 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
707 if(track) delete track;
708 }
709 continue;
710 }
711
712 ptMC = particle->Pt();
713
714 if(fTrackType==1 || fTrackType==2)
715 track->GetImpactParametersTPC(dca2D,dcaZ); //TPConly
716 else
717 track->GetImpactParameters(dca2D,dcaZ); //Global
718
719 UChar_t itsMap = track->GetITSClusterMap();
720 for (Int_t i=0; i < 6; i++) {
721 if (itsMap & (1 << i))
722 nPointITS ++;
723 }
724
725 // fPtAll->Fill(pt);
726 fPtAllMC->Fill(ptMC);
727
728 if (fTrackCuts->AcceptTrack(track)) {
729
730 if(fTrackType==7) {
731 if(fTrackCutsReject ) {
732 if(fTrackCutsReject->AcceptTrack(track) ) {
733 if(track) delete track;
734 continue;
735 }
736 }
737
738 if(esdtrack->GetConstrainedParam())
739 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
740 }
741
742 pt = track->Pt();
743 phi = track->Phi();
744
745 nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
746 chi2C = track->GetConstrainedChi2();
747 relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
748
749 fPtSel->Fill(pt);
750 if(track->GetLabel()<0) {
751 fPtSelFakes->Fill(pt);
752 fNPointTPCFakes->Fill(track->GetTPCNcls());
753 }
754 fPtSelMC->Fill(ptMC);
755 fPtAllvsPtMC->Fill(ptMC,pt);
756 fPtAllminPtMCvsPtMC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
757 fPtAllminPtMCvsPtAll->Fill(pt,(1./pt-1./ptMC)/(1./ptMC) );
758 fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
759 fPtAllminPtMCvsPtAllNPointTPCIter1->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNclsIter1());
760 if(track->GetTPCNcls()>0.)
761 fPtAllminPtMCvsPtAllChi2TPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCchi2()/track->GetTPCNcls());
762 if(track->GetTPCNclsIter1()>0.)
763 fPtAllminPtMCvsPtAllChi2TPCIter1->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCchi2Iter1()/track->GetTPCNclsIter1());
764 fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
765 fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
766 fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
767 fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
768 fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
769 fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
770 fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
771
772 fPtAllvsPtMCvsMult->Fill(ptMC,pt,mult);
773 fPtAllminPtMCvsPtAllvsMult->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC), mult);
774
775 //Check if track is reconstructed multiple times
776 /*
777 int multCounter = 1;
778 for (Int_t iTrack2 = iTrack+1; iTrack2 < nTracks; iTrack2++) {
779 // AliESDtrack *track2 = GetTrackForAnalysis(iTrack2);
780 AliESDtrack *track2;
781 AliESDtrack *esdtrack2 = fESD->GetTrack(iTrack2);
782 if(!esdtrack2) continue;
783
784 if(fTrackType==1)
785 track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
786 else if(fTrackType==2) {
787 track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
788 if(!track2) {
789 continue;
790 }
791 AliExternalTrackParam exParam2;
792 Bool_t relate = track2->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam2);
793 if( !relate ) {
794 delete track2;
795 continue;
796 }
797 track2->Set(exParam2.GetX(),exParam2.GetAlpha(),exParam2.GetParameter(),exParam2.GetCovariance());
798 }
799 else
800 track2 = esdtrack2;
801 if(!track2) {
802 continue;
803 }
804
805 if (fTrackCuts->AcceptTrack(track2)) {
806 Int_t label2 = TMath::Abs(track2->GetLabel());
807 if(label==label2) {
808 fNPointTPCMultRec->Fill(track->GetTPCNcls());
809 fNPointTPCMultRec->Fill(track2->GetTPCNcls());
810 fDeltaPtMultRec->Fill(track->Pt(),track->Pt()-track2->Pt());
811 multCounter++;
812 }
813 }
814 if(fTrackType==1 || fTrackType==2) delete track2;
815 }//track2 loop
816 if(multCounter>1) fMultRec->Fill(multCounter);
817 */
818
819 }//fTrackCuts selection
820
821
822 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
823 if(track) delete track;
824 }
825
826 }//ESD track loop
827
828 // Post output data
829 PostData(0, fHistList);
830
831}
832//________________________________________________________________________
833Bool_t AliPWG4HighPtQAMC::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
834 //
835 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
836 // This is to called in Notify and should provide the path to the AOD/ESD file
837 // Copied from AliAnalysisTaskJetSpectrum2
838 //
839
840 TString file(currFile);
841 fXsec = 0;
842 fTrials = 1;
843
844 if(file.Contains("root_archive.zip#")){
845 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
846 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
847 file.Replace(pos+1,20,"");
848 }
849 else {
850 // not an archive take the basename....
851 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
852 }
853 // Printf("%s",file.Data());
854
855
856 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
857 if(!fxsec){
858 // next trial fetch the histgram file
859 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
860 if(!fxsec){
861 // not a severe condition but inciate that we have no information
862 return kFALSE;
863 }
864 else{
865 // find the tlist we want to be independtent of the name so use the Tkey
866 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
867 if(!key){
868 fxsec->Close();
869 return kFALSE;
870 }
871 TList *list = dynamic_cast<TList*>(key->ReadObj());
872 if(!list){
873 fxsec->Close();
874 return kFALSE;
875 }
876 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
877 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
878 fxsec->Close();
879 }
880 } // no tree pyxsec.root
881 else {
882 TTree *xtree = (TTree*)fxsec->Get("Xsection");
883 if(!xtree){
884 fxsec->Close();
885 return kFALSE;
886 }
887 UInt_t ntrials = 0;
888 Double_t xsection = 0;
889 xtree->SetBranchAddress("xsection",&xsection);
890 xtree->SetBranchAddress("ntrials",&ntrials);
891 xtree->GetEntry(0);
892 fTrials = ntrials;
893 fXsec = xsection;
894 fxsec->Close();
895 }
896 return kTRUE;
897}
898//________________________________________________________________________
899Bool_t AliPWG4HighPtQAMC::Notify()
900{
901 //
902 // Implemented Notify() to read the cross sections
903 // and number of trials from pyxsec.root
904 // Copied from AliAnalysisTaskJetSpectrum2
905 //
906
907 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
908 Float_t xsection = 0;
909 Float_t ftrials = 1;
910
911 fAvgTrials = 1;
912 if(tree){
913 TFile *curfile = tree->GetCurrentFile();
914 if (!curfile) {
915 Error("Notify","No current file");
916 return kFALSE;
917 }
918 if(!fh1Xsec||!fh1Trials){
919 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
920 return kFALSE;
921 }
922 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
923 fh1Xsec->Fill("<#sigma>",xsection);
924 // construct a poor man average trials
925 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
926 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
927 }
928 return kTRUE;
929}
930
931//________________________________________________________________________
932AliGenPythiaEventHeader* AliPWG4HighPtQAMC::GetPythiaEventHeader(AliMCEvent *mcEvent){
933
934 if(!mcEvent)return 0;
935 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
936 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
937 if(!pythiaGenHeader){
938 // cocktail ??
939 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
940
941 if (!genCocktailHeader) {
942 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
943 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
944 return 0;
945 }
946 TList* headerList = genCocktailHeader->GetHeaders();
947 for (Int_t i=0; i<headerList->GetEntries(); i++) {
948 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
949 if (pythiaGenHeader)
950 break;
951 }
952 if(!pythiaGenHeader){
953 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
954 return 0;
955 }
956 }
957 return pythiaGenHeader;
958
959}
960
961//________________________________________________________________________
962void AliPWG4HighPtQAMC::Terminate(Option_t *)
963{
964 // The Terminate() function is the last function to be called during
965 // a query. It always runs on the client, it can be used to present
966 // the results graphically or save the results to file.
967
968}
969
970#endif