Bug fix - Coverity warnings 16767 and 16768
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliPWG4HighPtTrackQA.cxx
CommitLineData
d756027f 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 stores QA variables as function of pT for different type
18// of tracks and track selection criteria
19// Output: Histograms for different set of cuts
20//-----------------------------------------------------------------------
21// Author : Marta Verweij - UU
22//-----------------------------------------------------------------------
23
24#ifndef ALIPWG4HIGHPTTRACKQA_CXX
25#define ALIPWG4HIGHPTTRACKQA_CXX
26
27#include "AliPWG4HighPtTrackQA.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#include "TBits.h"
40
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 "AliCentrality.h"
53#include "AliAODVertex.h"
54#include "AliAODEvent.h"
55//#include "AliAnalysisHelperJetTasks.h"
56
57using namespace std; //required for resolving the 'cout' symbol
58
59ClassImp(AliPWG4HighPtTrackQA)
60
61AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA()
62: AliAnalysisTaskSE(),
63 fDataType(kESD),
64 fEvent(0x0),
65 fESD(0x0),
66 fVtx(0x0),
67 fTrackCuts(0),
68 fTrackType(0),
69 fFilterMask(0),
2b553e6f 70 fSigmaConstrainedMax(5.),
d756027f 71 fPtMax(100.),
72 fIsPbPb(0),
73 fCentClass(10),
2b553e6f 74 fNVariables(18),
d756027f 75 fVariables(0x0),
76 fAvgTrials(1),
77 fNEventAll(0),
78 fNEventSel(0),
79 fNEventReject(0),
80 fh1Centrality(0x0),
81 fh1Xsec(0),
82 fh1Trials(0),
83 fh1PtHard(0),
84 fh1PtHardTrials(0),
85 fh1NTracksAll(0x0),
86 fh1NTracksReject(0x0),
87 fh1NTracksSel(0x0),
88 fPtAll(0),
89 fPtSel(0),
90 fPtPhi(0x0),
91 fPtEta(0x0),
92 fPtDCA2D(0x0),
93 fPtDCAZ(0x0),
94 fPtNClustersTPC(0x0),
95 fPtNPointITS(0x0),
96 fPtChi2C(0x0),
97 fPtNSigmaToVertex(0x0),
98 fPtRelUncertainty1Pt(0x0),
2b553e6f 99 fPtUncertainty1Pt(0x0),
d756027f 100 fPtChi2PerClusterTPC(0x0),
101 fPtNCrossedRows(0x0),
102 fPtNCrossedRowsNClusF(0x0),
103 fPtNCrRNCrRNClusF(0x0),
2b553e6f 104 fPtSigmaY2(0x0),
105 fPtSigmaZ2(0x0),
106 fPtSigmaSnp2(0x0),
107 fPtSigmaTgl2(0x0),
108 fPtSigma1Pt2(0x0),
d756027f 109 fHistList(0)
110{
2b553e6f 111 SetNVariables(18);
d756027f 112}
113//________________________________________________________________________
114AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA(const char *name):
115 AliAnalysisTaskSE(name),
116 fDataType(kESD),
117 fEvent(0x0),
118 fESD(0x0),
119 fVtx(0x0),
120 fTrackCuts(),
121 fTrackType(0),
122 fFilterMask(0),
2b553e6f 123 fSigmaConstrainedMax(5.),
d756027f 124 fPtMax(100.),
125 fIsPbPb(0),
126 fCentClass(10),
2b553e6f 127 fNVariables(18),
d756027f 128 fVariables(0x0),
129 fAvgTrials(1),
130 fNEventAll(0),
131 fNEventSel(0),
132 fNEventReject(0),
133 fh1Centrality(0x0),
134 fh1Xsec(0),
135 fh1Trials(0),
136 fh1PtHard(0),
137 fh1PtHardTrials(0),
138 fh1NTracksAll(0x0),
139 fh1NTracksReject(0x0),
140 fh1NTracksSel(0x0),
141 fPtAll(0),
142 fPtSel(0),
143 fPtPhi(0x0),
144 fPtEta(0x0),
145 fPtDCA2D(0x0),
146 fPtDCAZ(0x0),
147 fPtNClustersTPC(0x0),
148 fPtNPointITS(0x0),
149 fPtChi2C(0x0),
150 fPtNSigmaToVertex(0x0),
151 fPtRelUncertainty1Pt(0x0),
2b553e6f 152 fPtUncertainty1Pt(0x0),
d756027f 153 fPtChi2PerClusterTPC(0x0),
154 fPtNCrossedRows(0x0),
155 fPtNCrossedRowsNClusF(0x0),
156 fPtNCrRNCrRNClusF(0x0),
2b553e6f 157 fPtSigmaY2(0x0),
158 fPtSigmaZ2(0x0),
159 fPtSigmaSnp2(0x0),
160 fPtSigmaTgl2(0x0),
161 fPtSigma1Pt2(0x0),
d756027f 162 fHistList(0)
163{
164 //
165 // Constructor. Initialization of Inputs and Outputs
166 //
167 AliDebug(2,Form("AliPWG4HighPtTrackQA Calling Constructor"));
168
2b553e6f 169 SetNVariables(18);
d756027f 170
171 // Input slot #0 works with a TChain ESD
172 DefineInput(0, TChain::Class());
173 // Output slot #1 write into a TList
174 DefineOutput(1, TList::Class());
175}
176
177//________________________________________________________________________
178void AliPWG4HighPtTrackQA::UserCreateOutputObjects() {
179 //Create output objects
180 AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserCreateOutputObjects \n"));
181
182 Bool_t oldStatus = TH1::AddDirectoryStatus();
183 TH1::AddDirectory(kFALSE);
184
185 OpenFile(1);
186 fHistList = new TList();
187 fHistList->SetOwner(kTRUE);
188
189 Float_t fgkPtMin = 0.;
190 Float_t fgkPtMax = fPtMax;
191
192 Float_t ptBinEdges[2][2];
193 ptBinEdges[0][0] = 10.;
194 ptBinEdges[0][1] = 1.;
195 ptBinEdges[1][0] = 20.;
196 ptBinEdges[1][1] = 2.;
197 Float_t binWidth3 = 5.;
198 if(fPtMax>100.) {
199 ptBinEdges[0][0] = 100.;
200 ptBinEdges[0][1] = 5.;
201 ptBinEdges[1][0] = 300.;
202 ptBinEdges[1][1] = 10.;
203 binWidth3 = 20.;
204 }
205
206 const Float_t ptmin1 = fgkPtMin;
207 const Float_t ptmax1 = ptBinEdges[0][0];
208 const Float_t ptmin2 = ptmax1 ;
209 const Float_t ptmax2 = ptBinEdges[1][0];
210 const Float_t ptmin3 = ptmax2 ;
211 const Float_t ptmax3 = fgkPtMax;
212 const Int_t nbin11 = (int)((ptmax1-ptmin1)/ptBinEdges[0][1]);
213 const Int_t nbin12 = (int)((ptmax2-ptmin2)/ptBinEdges[1][1])+nbin11;
214 const Int_t nbin13 = (int)((ptmax3-ptmin3)/binWidth3)+nbin12;
215 Int_t fgkNPtBins=nbin13;
216 //Create array with low edges of each bin
217 Double_t *binsPt=new Double_t[fgkNPtBins+1];
218 for(Int_t i=0; i<=fgkNPtBins; i++) {
219 if(i<=nbin11) binsPt[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
220 if(i<=nbin12 && i>nbin11) binsPt[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
221 if(i<=nbin13 && i>nbin12) binsPt[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
222 }
223
224 Int_t fgkNPhiBins = 18*6;
225 Float_t kMinPhi = 0.;
226 Float_t kMaxPhi = 2.*TMath::Pi();
227 Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
228 for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
229
230 Int_t fgkNEtaBins=20;
231 Float_t fgkEtaMin = -1.;
232 Float_t fgkEtaMax = 1.;
233 Double_t *binsEta=new Double_t[fgkNEtaBins+1];
234 for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
235
236 Int_t fgkNNClustersTPCBins=80;
237 Float_t fgkNClustersTPCMin = 0.5;
238 Float_t fgkNClustersTPCMax = 160.5;
239 Double_t *binsNClustersTPC=new Double_t[fgkNNClustersTPCBins+1];
240 for(Int_t i=0; i<=fgkNNClustersTPCBins; i++) binsNClustersTPC[i]=(Double_t)fgkNClustersTPCMin + (fgkNClustersTPCMax-fgkNClustersTPCMin)/fgkNNClustersTPCBins*(Double_t)i ;
241
242 Int_t fgkNDCA2DBins=80;
243 Float_t fgkDCA2DMin = -0.2;
244 Float_t fgkDCA2DMax = 0.2;
2b553e6f 245 if(fTrackType==1 || fTrackType==2) {
246 fgkDCA2DMin = -2.;
247 fgkDCA2DMax = 2.;
248 }
d756027f 249 Double_t *binsDCA2D=new Double_t[fgkNDCA2DBins+1];
250 for(Int_t i=0; i<=fgkNDCA2DBins; i++) binsDCA2D[i]=(Double_t)fgkDCA2DMin + (fgkDCA2DMax-fgkDCA2DMin)/fgkNDCA2DBins*(Double_t)i ;
251
252 Int_t fgkNDCAZBins=80;
253 Float_t fgkDCAZMin = -2.;
254 Float_t fgkDCAZMax = 2.;
2b553e6f 255 if(fTrackType==1 || fTrackType==2) {
256 fgkDCAZMin = -5.;
257 fgkDCAZMax = 5.;
258 }
d756027f 259 Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
260 for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
261
262 Int_t fgkNNPointITSBins=9;
263 Float_t fgkNPointITSMin = -0.5;
264 Float_t fgkNPointITSMax = 8.5;
265 Double_t *binsNPointITS=new Double_t[fgkNNPointITSBins+1];
266 for(Int_t i=0; i<=fgkNNPointITSBins; i++) binsNPointITS[i]=(Double_t)fgkNPointITSMin + (fgkNPointITSMax-fgkNPointITSMin)/fgkNNPointITSBins*(Double_t)i ;
267
268 Int_t fgkNNSigmaToVertexBins=40;
269 Float_t fgkNSigmaToVertexMin = 0.;
270 Float_t fgkNSigmaToVertexMax = 8.;
271 Double_t *binsNSigmaToVertex=new Double_t[fgkNNSigmaToVertexBins+1];
272 for(Int_t i=0; i<=fgkNNSigmaToVertexBins; i++) binsNSigmaToVertex[i]=(Double_t)fgkNSigmaToVertexMin + (fgkNSigmaToVertexMax-fgkNSigmaToVertexMin)/fgkNNSigmaToVertexBins*(Double_t)i ;
273
274 Int_t fgkNChi2CBins=20;
275 Float_t fgkChi2CMin = 0.;
2b553e6f 276 Float_t fgkChi2CMax = 100.;
d756027f 277 Double_t *binsChi2C=new Double_t[fgkNChi2CBins+1];
278 for(Int_t i=0; i<=fgkNChi2CBins; i++) binsChi2C[i]=(Double_t)fgkChi2CMin + (fgkChi2CMax-fgkChi2CMin)/fgkNChi2CBins*(Double_t)i ;
279
280 Int_t fgkNRel1PtUncertaintyBins=30;
281 Float_t fgkRel1PtUncertaintyMin = 0.;
282 Float_t fgkRel1PtUncertaintyMax = 0.3;
2b553e6f 283 if(fTrackType==1 || fTrackType==2) fgkRel1PtUncertaintyMax = 0.5;
d756027f 284 Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
285 for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
286
2b553e6f 287 Int_t fgkNUncertainty1PtBins = 30;
288 Float_t fgkUncertainty1PtMin = 0.;
289 Float_t fgkUncertainty1PtMax = 0.1;
290 if(fTrackType==1 || fTrackType==2) fgkUncertainty1PtMax = 0.2;
291 Double_t *binsUncertainty1Pt=new Double_t[fgkNUncertainty1PtBins+1];
292 for(Int_t i=0; i<=fgkNUncertainty1PtBins; i++) binsUncertainty1Pt[i]=(Double_t)fgkUncertainty1PtMin + (fgkUncertainty1PtMax-fgkUncertainty1PtMin)/fgkNUncertainty1PtBins*(Double_t)i ;
293
d756027f 294 Float_t fgkChi2PerClusMin = 0.;
295 Float_t fgkChi2PerClusMax = 4.;
296 Int_t fgkNChi2PerClusBins = (int)(fgkChi2PerClusMax*10.);
297 Double_t *binsChi2PerClus=new Double_t[fgkNChi2PerClusBins+1];
298 for(Int_t i=0; i<=fgkNChi2PerClusBins; i++) binsChi2PerClus[i]=(Double_t)fgkChi2PerClusMin + (fgkChi2PerClusMax-fgkChi2PerClusMin)/fgkNChi2PerClusBins*(Double_t)i ;
299
300 Int_t fgkNCrossedRowsNClusFBins = 50;
301 Float_t fgkNCrossedRowsNClusFMin = 0.;
302 Float_t fgkNCrossedRowsNClusFMax = 1.;
303 Double_t *binsNCrossedRowsNClusF=new Double_t[fgkNCrossedRowsNClusFBins+1];
304 for(Int_t i=0; i<=fgkNCrossedRowsNClusFBins; i++) binsNCrossedRowsNClusF[i]=(Double_t)fgkNCrossedRowsNClusFMin + (fgkNCrossedRowsNClusFMax-fgkNCrossedRowsNClusFMin)/fgkNCrossedRowsNClusFBins*(Double_t)i ;
305
2b553e6f 306 Int_t fgkN1PtBins = 50;
307 Float_t fgk1PtMin = 0.;
308 Float_t fgk1PtMax = 6.;
309 Double_t *bins1Pt=new Double_t[fgkN1PtBins+1];
310 for(Int_t i=0; i<=fgkN1PtBins; i++) bins1Pt[i]=(Double_t)fgk1PtMin + (fgk1PtMax-fgk1PtMin)/fgkN1PtBins*(Double_t)i ;
311
312 Int_t fgkNSigmaY2Bins = 50;
313 Float_t fgkSigmaY2Min = 0.;
314 Float_t fgkSigmaY2Max = 2.;
315 Double_t *binsSigmaY2=new Double_t[fgkNSigmaY2Bins+1];
316 for(Int_t i=0; i<=fgkNSigmaY2Bins; i++) binsSigmaY2[i]=(Double_t)fgkSigmaY2Min + (fgkSigmaY2Max-fgkSigmaY2Min)/fgkNSigmaY2Bins*(Double_t)i ;
317
318 Int_t fgkNSigmaZ2Bins = 50;
319 Float_t fgkSigmaZ2Min = 0.;
320 Float_t fgkSigmaZ2Max = 2.;
321 Double_t *binsSigmaZ2=new Double_t[fgkNSigmaZ2Bins+1];
322 for(Int_t i=0; i<=fgkNSigmaZ2Bins; i++) binsSigmaZ2[i]=(Double_t)fgkSigmaZ2Min + (fgkSigmaZ2Max-fgkSigmaZ2Min)/fgkNSigmaZ2Bins*(Double_t)i ;
323
324 Int_t fgkNSigmaSnp2Bins = 50;
325 Float_t fgkSigmaSnp2Min = 0.;
326 Float_t fgkSigmaSnp2Max = 2.;
327 Double_t *binsSigmaSnp2=new Double_t[fgkNSigmaSnp2Bins+1];
328 for(Int_t i=0; i<=fgkNSigmaSnp2Bins; i++) binsSigmaSnp2[i]=(Double_t)fgkSigmaSnp2Min + (fgkSigmaSnp2Max-fgkSigmaSnp2Min)/fgkNSigmaSnp2Bins*(Double_t)i ;
329
330 Int_t fgkNSigmaTgl2Bins = 50;
331 Float_t fgkSigmaTgl2Min = 0.;
332 Float_t fgkSigmaTgl2Max = 2.;
333 Double_t *binsSigmaTgl2=new Double_t[fgkNSigmaTgl2Bins+1];
334 for(Int_t i=0; i<=fgkNSigmaTgl2Bins; i++) binsSigmaTgl2[i]=(Double_t)fgkSigmaTgl2Min + (fgkSigmaTgl2Max-fgkSigmaTgl2Min)/fgkNSigmaTgl2Bins*(Double_t)i ;
335
336 Int_t fgkNSigma1Pt2Bins = 50;
337 Float_t fgkSigma1Pt2Min = 0.;
338 Float_t fgkSigma1Pt2Max = 2.;
339 Double_t *binsSigma1Pt2=new Double_t[fgkNSigma1Pt2Bins+1];
340 for(Int_t i=0; i<=fgkNSigma1Pt2Bins; i++) binsSigma1Pt2[i]=(Double_t)fgkSigma1Pt2Min + (fgkSigma1Pt2Max-fgkSigma1Pt2Min)/fgkNSigma1Pt2Bins*(Double_t)i ;
341
342
d756027f 343 fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
344 fHistList->Add(fNEventAll);
345 fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
346 fHistList->Add(fNEventSel);
347 fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
348 //Set labels
349 fNEventReject->Fill("noESD",0);
350 fNEventReject->Fill("Trigger",0);
351 fNEventReject->Fill("NTracks<2",0);
352 fNEventReject->Fill("noVTX",0);
353 fNEventReject->Fill("VtxStatus",0);
354 fNEventReject->Fill("NCont<2",0);
355 fNEventReject->Fill("ZVTX>10",0);
356 fNEventReject->Fill("cent",0);
357 fNEventReject->Fill("cent>90",0);
358 fHistList->Add(fNEventReject);
359
360 fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
361 fHistList->Add(fh1Centrality);
362
363 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
364 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
365 fHistList->Add(fh1Xsec);
366
367 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
368 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
369 fHistList->Add(fh1Trials);
370
371 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
372 fHistList->Add(fh1PtHard);
373 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
374 fHistList->Add(fh1PtHardTrials);
375
376 fh1NTracksAll = new TH1F("fh1NTracksAll","fh1NTracksAll",1,-0.5,0.5);
377 fHistList->Add(fh1NTracksAll);
378
379 fh1NTracksReject = new TH1F("fh1NTracksReject","fh1NTracksReject",1,-0.5,0.5);
380 fh1NTracksReject->Fill("noESDtrack",0);
381 fh1NTracksReject->Fill("noTPCInner",0);
382 fh1NTracksReject->Fill("FillTPC",0);
383 fh1NTracksReject->Fill("noTPConly",0);
384 fh1NTracksReject->Fill("relate",0);
385 fh1NTracksReject->Fill("trackCuts",0);
386 fh1NTracksReject->Fill("laser",0);
2b553e6f 387 fh1NTracksReject->Fill("chi2",0);
d756027f 388 fHistList->Add(fh1NTracksReject);
389
390 fh1NTracksSel = new TH1F("fh1NTracksSel","fh1NTracksSel",1,-0.5,0.5);
391 fHistList->Add(fh1NTracksSel);
392
393
394 fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
395 fHistList->Add(fPtAll);
396 fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
397 fHistList->Add(fPtSel);
398
399 fPtPhi = new TH2F("fPtPhi","fPtPhi",fgkNPtBins,binsPt,fgkNPhiBins,binsPhi);
400 fHistList->Add(fPtPhi);
401
402 fPtEta = new TH2F("fPtEta","fPtEta",fgkNPtBins,binsPt,fgkNEtaBins,binsEta);
403 fHistList->Add(fPtEta);
404
405 fPtDCA2D = new TH2F("fPtDCA2D","fPtDCA2D",fgkNPtBins,binsPt,fgkNDCA2DBins,binsDCA2D);
406 fHistList->Add(fPtDCA2D);
407
408 fPtDCAZ = new TH2F("fPtDCAZ","fPtDCAZ",fgkNPtBins,binsPt,fgkNDCAZBins,binsDCAZ);
409 fHistList->Add(fPtDCAZ);
410
411 fPtNClustersTPC = new TH2F("fPtNClustersTPC","fPtNClustersTPC",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
412 fHistList->Add(fPtNClustersTPC);
413
414 fPtNPointITS = new TH2F("fPtNPointITS","fPtNPointITS",fgkNPtBins,binsPt,fgkNNPointITSBins,binsNPointITS);
415 fHistList->Add(fPtNPointITS);
416
417 fPtChi2C = new TH2F("fPtChi2C","fPtChi2C",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
418 fHistList->Add(fPtChi2C);
419
420 fPtNSigmaToVertex = new TH2F("fPtNSigmaToVertex","fPtNSigmaToVertex",fgkNPtBins,binsPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
421 fHistList->Add(fPtNSigmaToVertex);
422
423 fPtRelUncertainty1Pt = new TH2F("fPtRelUncertainty1Pt","fPtRelUncertainty1Pt",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
424 fHistList->Add(fPtRelUncertainty1Pt);
2b553e6f 425
426 fPtUncertainty1Pt = new TH2F("fPtUncertainty1Pt","fPtUncertainty1Pt",fgkNPtBins,binsPt,fgkNUncertainty1PtBins,binsUncertainty1Pt);
427 fHistList->Add(fPtUncertainty1Pt);
d756027f 428
429 fPtChi2PerClusterTPC = new TH2F("fPtChi2PerClusterTPC","fPtChi2PerClusterTPC",fgkNPtBins,binsPt,fgkNChi2PerClusBins,binsChi2PerClus);
430 fHistList->Add(fPtChi2PerClusterTPC);
431
432 fPtNCrossedRows = new TH2F("fPtNCrossedRows","fPtNCrossedRows",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
433 fHistList->Add(fPtNCrossedRows);
434
435 fPtNCrossedRowsNClusF = new TH2F("fPtNCrossedRowsNClusF","fPtNCrossedRowsNClusF",fgkNPtBins,binsPt,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
436 fHistList->Add(fPtNCrossedRowsNClusF);
437
438 fPtNCrRNCrRNClusF = new TH3F("fPtNCrRNCrRNClusF","fPtNCrRNCrRNClusF",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
439 fHistList->Add(fPtNCrRNCrRNClusF);
2b553e6f 440
441 fPtSigmaY2 = new TH2F("fPtSigmaY2","fPtSigmaY2",fgkN1PtBins,bins1Pt,fgkNSigmaZ2Bins,binsSigmaY2);
442 fHistList->Add(fPtSigmaY2);
443
444 fPtSigmaZ2 = new TH2F("fPtSigmaZ2","fPtSigmaZ2",fgkN1PtBins,bins1Pt,fgkNSigmaZ2Bins,binsSigmaZ2);
445 fHistList->Add(fPtSigmaZ2);
446
447 fPtSigmaSnp2 = new TH2F("fPtSigmaSnp2","fPtSigmaSnp2",fgkN1PtBins,bins1Pt,fgkNSigmaZ2Bins,binsSigmaSnp2);
448 fHistList->Add(fPtSigmaSnp2);
449
450 fPtSigmaTgl2 = new TH2F("fPtSigmaTgl2","fPtSigmaTgl2",fgkN1PtBins,bins1Pt,fgkNSigmaZ2Bins,binsSigmaTgl2);
451 fHistList->Add(fPtSigmaTgl2);
452
453 fPtSigma1Pt2 = new TH2F("fPtSigma1Pt2","fPtSigma1Pt2",fgkN1PtBins,bins1Pt,fgkNSigmaZ2Bins,binsSigma1Pt2);
454 fHistList->Add(fPtSigma1Pt2);
455
d756027f 456 TH1::AddDirectory(oldStatus);
457
458 PostData(1, fHistList);
459
460 if(binsPhi) delete [] binsPhi;
461 if(binsPt) delete [] binsPt;
462 if(binsNClustersTPC) delete [] binsNClustersTPC;
463 if(binsDCA2D) delete [] binsDCA2D;
464 if(binsDCAZ) delete [] binsDCAZ;
465 if(binsNPointITS) delete [] binsNPointITS;
466 if(binsNSigmaToVertex) delete [] binsNSigmaToVertex;
467 if(binsChi2C) delete [] binsChi2C;
468 if(binsEta) delete [] binsEta;
469 if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
2b553e6f 470 if(binsUncertainty1Pt) delete [] binsUncertainty1Pt;
d756027f 471 if(binsChi2PerClus) delete [] binsChi2PerClus;
472 if(binsChi2PerClus) delete [] binsNCrossedRowsNClusF;
2b553e6f 473 if(bins1Pt) delete [] bins1Pt;
474 if(binsSigmaY2) delete [] binsSigmaY2;
475 if(binsSigmaZ2) delete [] binsSigmaZ2;
476 if(binsSigmaSnp2) delete [] binsSigmaSnp2;
477 if(binsSigmaTgl2) delete [] binsSigmaTgl2;
478 if(binsSigma1Pt2) delete [] binsSigma1Pt2;
d756027f 479}
480
481//________________________________________________________________________
482Bool_t AliPWG4HighPtTrackQA::SelectEvent() {
483 //
484 // Decide if event should be selected for analysis
485 //
486
487 // Checks following requirements:
488 // - fEvent available
489 // - trigger info from AliPhysicsSelection
490 // - MCevent available
491 // - number of reconstructed tracks > 1
492 // - primary vertex reconstructed
493 // - z-vertex < 10 cm
494 // - centrality in case of PbPb
495
496 Bool_t selectEvent = kTRUE;
497
498 //fEvent object available?
499 if (!fEvent) {
500 AliDebug(2,Form("ERROR: fInputEvent not available\n"));
501 fNEventReject->Fill("noAliVEvent",1);
502 selectEvent = kFALSE;
503 return selectEvent;
504 }
505
506 //Check if number of reconstructed tracks is larger than 1
507 if(!fEvent->GetNumberOfTracks() || fEvent->GetNumberOfTracks()<2) {
508 fNEventReject->Fill("NTracks<2",1);
509 selectEvent = kFALSE;
510 return selectEvent;
511 }
512
513 //Check if vertex is reconstructed
1ea145ef 514 if(fDataType==kESD&&dynamic_cast<AliESDEvent*>(fEvent)) {
515 fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexSPD();
d756027f 516
517 if(!fVtx) {
518 fNEventReject->Fill("noVTX",1);
519 selectEvent = kFALSE;
520 return selectEvent;
521 }
522
523 if(!fVtx->GetStatus()) {
524 fNEventReject->Fill("VtxStatus",1);
525 selectEvent = kFALSE;
526 return selectEvent;
527 }
528
529 // Need vertex cut
530 if(fVtx->GetNContributors()<2) {
531 fNEventReject->Fill("NCont<2",1);
532 selectEvent = kFALSE;
533 return selectEvent;
534 }
535
536 //Check if z-vertex < 10 cm
537 double primVtx[3];
538 fVtx->GetXYZ(primVtx);
539 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
540 fNEventReject->Fill("ZVTX>10",1);
541 selectEvent = kFALSE;
542 return selectEvent;
543 }
544 }
1ea145ef 545 else if(fDataType==kAOD&&dynamic_cast<AliAODEvent*>(fEvent)) {
546 const AliAODVertex *vtx = ((AliAODEvent*)fEvent)->GetPrimaryVertexSPD();
d756027f 547 if(!vtx) {
548 fNEventReject->Fill("noVTX",1);
549 selectEvent = kFALSE;
550 return selectEvent;
551 }
552
553 // Need vertex cut
554 if(vtx->GetNContributors()<2) {
555 fNEventReject->Fill("NCont<2",1);
556 selectEvent = kFALSE;
557 return selectEvent;
558 }
559
560 //Check if z-vertex < 10 cm
561 double primVtx[3];
562 vtx->GetXYZ(primVtx);
563 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
564 fNEventReject->Fill("ZVTX>10",1);
565 selectEvent = kFALSE;
566 return selectEvent;
567 }
568
569 }
570
571 //Centrality selection should only be done in case of PbPb
572 if(IsPbPb()) {
573 Float_t cent = 0.;
574 if(fCentClass!=CalculateCentrality(fEvent) && fCentClass!=10) {
575 fNEventReject->Fill("cent",1);
576 selectEvent = kFALSE;
577 return selectEvent;
578 }
579 else {
580 if(fDataType==kESD) {
581 if(dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()) {
582 cent = dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()->GetCentralityPercentile("V0M");
583 }
584 }
585 else if(fDataType==kAOD) {
586 if(dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality())
587 cent = dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality();
588 }
589 if(cent>90.) {
590 fNEventReject->Fill("cent>90",1);
591 selectEvent = kFALSE;
592 return selectEvent;
593 }
594 fh1Centrality->Fill(cent);
595 }
596 }
597
598 return selectEvent;
599
600}
601
602//________________________________________________________________________
603Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliVEvent *ev){
604 if(fDataType==kESD)
605 return CalculateCentrality(dynamic_cast<AliESDEvent*>(ev));
606 else if(fDataType==kAOD)
607 return CalculateCentrality(dynamic_cast<AliAODEvent*>(ev));
608 else
609 return 5;
610}
611
612//________________________________________________________________________
613Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliESDEvent *esd){
614
1ea145ef 615
d756027f 616 Float_t cent = 999;
1ea145ef 617
618 if(esd){
619 if(esd->GetCentrality()){
620 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
621 if(fDebug>3) cout << "centrality: " << cent << endl;
622 }
d756027f 623 }
624
625 if(cent>80)return 4;
626 if(cent>50)return 3;
627 if(cent>30)return 2;
628 if(cent>10)return 1;
629 return 0;
630
631}
632
633//________________________________________________________________________
634Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliAODEvent *aod){
635
1ea145ef 636 if(!aod)return 4;
d756027f 637 Float_t cent = aod->GetHeader()->GetCentrality();
638 cout << "centrality: " << cent << endl;
639 if(cent>80)return 4;
640 if(cent>50)return 3;
641 if(cent>30)return 2;
642 if(cent>10)return 1;
643 return 0;
644
645}
646
647//________________________________________________________________________
648void AliPWG4HighPtTrackQA::UserExec(Option_t *) {
649 // Main loop
650 // Called for each event
651 AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserExec \n"));
652
653 fEvent = InputEvent();
654 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
655
656 // All events without selection
657 fNEventAll->Fill(0.);
658
659 if(!SelectEvent()) {
660 // Post output data
661 PostData(1, fHistList);
662 return;
663 }
664
665
666 //Need to keep track of selected events
667 fNEventSel->Fill(0.);
668
669 fVariables = new TArrayF(fNVariables);
670
671 if(fDataType==kESD) DoAnalysisESD();
672 if(fDataType==kAOD) DoAnalysisAOD();
673
674 //Delete old fVariables
675 if(fVariables) delete fVariables;
676
677 // Post output data
678 PostData(1, fHistList);
679
680}
681
682//________________________________________________________________________
683void AliPWG4HighPtTrackQA::DoAnalysisESD() {
684
685 if(!fESD) {
686 PostData(1, fHistList);
687 return;
688 }
689
690 // ---- Get MC Header information (for MC productions in pThard bins) ----
691 Double_t ptHard = 0.;
692 Double_t nTrials = 1; // trials for MC trigger weight for real data
693
694 AliMCEventHandler *eventHandlerMC = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
695 if (eventHandlerMC) {
696
697 if(eventHandlerMC->MCEvent()){
698 AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(eventHandlerMC->MCEvent());
699 if(pythiaGenHeader){
700 nTrials = pythiaGenHeader->Trials();
701 ptHard = pythiaGenHeader->GetPtHard();
702
703 fh1PtHard->Fill(ptHard);
704 fh1PtHardTrials->Fill(ptHard,nTrials);
705
706 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
707 }
708 }
709 }
710
711 Int_t nTracks = fESD->GetNumberOfTracks();
712 AliDebug(2,Form("nTracks ESD%d", nTracks));
713
714 /*
715 Variables to be put in fVariables
716 0: pt
717 1: phi
718 2: eta
719 3: dca2D
720 4: dcaZ
721 5: nClustersTPC
722 6: nPointITS
723 7: chi2C
724 8: nSigmaToVertex
725 9: relUncertainty1Pt
726 10: chi2PerClusterTPC
727 11: #crossed rows
728 12: (#crossed rows)/(#findable clusters)
2b553e6f 729 13: SigmaY2
730 14: SigmaZ2
731 15: SigmaSnp2
732 16: SigmaTgl2
733 17: Sigma1Pt2
d756027f 734 */
735
736 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
737 fh1NTracksAll->Fill(0.);
738
739 //Get track for analysis
2b553e6f 740 AliESDtrack *track = 0x0;
d756027f 741 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
742 if(!esdtrack) {
743 fh1NTracksReject->Fill("noESDtrack",1);
744 continue;
745 }
746
747 if(fTrackType==1)
748 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
749 else if(fTrackType==2) {
750 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
751 if(!track) {
752 fh1NTracksReject->Fill("noTPConly",1);
d756027f 753 continue;
754 }
755 AliExternalTrackParam exParam;
756 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
757 if( !relate ) {
758 fh1NTracksReject->Fill("relate",1);
759 delete track;
760 continue;
761 }
762 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
763 }
764 else
765 track = esdtrack;
766
767 if(!track) {
d756027f 768 continue;
769 }
770
2b553e6f 771 if(fTrackType==2) {
772 //Cut on chi2 of constrained fit
773 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
774 fh1NTracksReject->Fill("chi2",1);
775 delete track;
776 continue;
777 }
778 }
779
d756027f 780 fPtAll->Fill(track->Pt());
781
782 if (!(fTrackCuts->AcceptTrack(track))) {
783 fh1NTracksReject->Fill("trackCuts",1);
784 if(fTrackType==1 || fTrackType==2) delete track;
785 continue;
786 }
787 if(track->GetTPCsignal()<10) { //Cut on laser tracks
788 fh1NTracksReject->Fill("laser",1);
789 if(fTrackType==1 || fTrackType==2) delete track;
790 continue;
791 }
792
793 fh1NTracksSel->Fill(0.);
794
795 fVariables->Reset(0.);
796
2b553e6f 797 fVariables->SetAt(track->Pt(),0);
798 fVariables->SetAt(track->Phi(),1);
799 fVariables->SetAt(track->Eta(),2);
d756027f 800
801 Float_t dca2D = 0.;
802 Float_t dcaz = 0.;
803 if(fTrackType==0) { //Global
804 track->GetImpactParameters(dca2D,dcaz);
805 }
806 else if(fTrackType==1 || fTrackType==2) { //TPConly
807 track->GetImpactParametersTPC(dca2D,dcaz);
808 }
809 fVariables->SetAt(dca2D,3);
2b553e6f 810 fVariables->SetAt(dcaz,4);
d756027f 811
812 fVariables->SetAt((float)track->GetTPCNcls(),5);
813
814 Int_t nPointITS = 0;
815 UChar_t itsMap = track->GetITSClusterMap();
816 for (Int_t i=0; i < 6; i++) {
817 if (itsMap & (1 << i))
818 nPointITS ++;
819 }
820 fVariables->SetAt((float)nPointITS,6);
2b553e6f 821 Float_t chi2C = (float)track->GetConstrainedChi2();
822 if(fTrackType==1 || fTrackType==2)
823 chi2C = (float)track->GetConstrainedChi2TPC();
824 fVariables->SetAt(chi2C,7);
d756027f 825 fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
826
827 fVariables->SetAt(TMath::Sqrt(track->GetSigma1Pt2())*fVariables->At(0),9);
828
829 if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
830
831 //cout << "#crossed rows (1): " << track->GetTPCClusterInfo(1) << endl;
832 //cout << "#crossed rows (2): " << track->GetTPCClusterInfo(2) << endl;
833 fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
834 Float_t crossedRowsTPCNClsF = track->GetTPCClusterInfo(2,0);
835 //if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
836 fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
2b553e6f 837 fVariables->SetAt(track->GetSigmaY2(),13);
838 fVariables->SetAt(track->GetSigmaZ2(),14);
839 fVariables->SetAt(track->GetSigmaSnp2(),15);
840 fVariables->SetAt(track->GetSigmaTgl2(),16);
841 fVariables->SetAt(track->GetSigma1Pt2(),17);
d756027f 842
843 FillHistograms();
844
845 // int mult = fTrackCuts->CountAcceptedTracks(fESD);
846
847 if(fTrackType==1 || fTrackType==2) delete track;
848
849 }//track loop
850
851}
852
853//________________________________________________________________________
854void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
855
856 AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1ea145ef 857 if(!aod)return;
d756027f 858 for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
859
860 AliAODTrack *aodtrack = aod->GetTrack(iTrack);
861 if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
862
863 fVariables->Reset(0.);
864
865 fVariables->SetAt(aodtrack->Pt(),0);
866 fVariables->SetAt(aodtrack->Phi(),1);
867 fVariables->SetAt(aodtrack->Eta(),2);
868
869 Double_t dca[2] = {1e6,1e6};
870 Double_t covar[3] = {1e6,1e6,1e6};
871 if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
872 fVariables->SetAt(dca[0],3);
873 fVariables->SetAt(dca[1],4);
874 }
875
876 fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
877 fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
878 fVariables->SetAt(0.,7);
879 fVariables->SetAt(0.,8);
880 fVariables->SetAt(0.,9);
881 fVariables->SetAt(aodtrack->Chi2perNDF(),10);
882 fVariables->SetAt(GetTPCClusterInfo(aodtrack,2),11);
883 Float_t crossedRowsTPCNClsF = 0.;
884 if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
885 fVariables->SetAt(crossedRowsTPCNClsF,12);
886
887 fPtAll->Fill(fVariables->At(0));
888
889 FillHistograms();
890
891 }
892
893}
894
895//________________________________________________________________________
896void AliPWG4HighPtTrackQA::FillHistograms() {
897
898 fPtSel->Fill(fVariables->At(0));
899 fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
900 fPtEta->Fill(fVariables->At(0),fVariables->At(2));
901 fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
902 fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
903 fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
904 fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
905 if(fDataType==kESD) {
906 fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
907 fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
908 fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(9));
2b553e6f 909 fPtUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)/fVariables->At(9));
910 fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
911 fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
912 fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
913 fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
914 fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
d756027f 915 }
916 fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
917 fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
918 fPtNCrossedRowsNClusF->Fill(fVariables->At(0),fVariables->At(12));
919 fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
920}
921
922//________________________________________________________________________
923Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
924 //
925 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
926 // This is to called in Notify and should provide the path to the AOD/ESD file
927 // Copied from AliAnalysisTaskJetSpectrum2
928 //
929
930 TString file(currFile);
931 fXsec = 0;
932 fTrials = 1;
933
934 if(file.Contains("root_archive.zip#")){
935 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
936 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
937 file.Replace(pos+1,20,"");
938 }
939 else {
940 // not an archive take the basename....
941 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
942 }
943 // Printf("%s",file.Data());
944
945
946 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
947 if(!fxsec){
948 // next trial fetch the histgram file
949 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
950 if(!fxsec){
951 // not a severe condition but inciate that we have no information
952 return kFALSE;
953 }
954 else{
955 // find the tlist we want to be independtent of the name so use the Tkey
956 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
957 if(!key){
958 fxsec->Close();
959 return kFALSE;
960 }
961 TList *list = dynamic_cast<TList*>(key->ReadObj());
962 if(!list){
963 fxsec->Close();
964 return kFALSE;
965 }
966 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
967 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
968 fxsec->Close();
969 }
970 } // no tree pyxsec.root
971 else {
972 TTree *xtree = (TTree*)fxsec->Get("Xsection");
973 if(!xtree){
974 fxsec->Close();
975 return kFALSE;
976 }
977 UInt_t ntrials = 0;
978 Double_t xsection = 0;
979 xtree->SetBranchAddress("xsection",&xsection);
980 xtree->SetBranchAddress("ntrials",&ntrials);
981 xtree->GetEntry(0);
982 fTrials = ntrials;
983 fXsec = xsection;
984 fxsec->Close();
985 }
986 return kTRUE;
987}
988//________________________________________________________________________
989Bool_t AliPWG4HighPtTrackQA::Notify()
990{
991 //
992 // Implemented Notify() to read the cross sections
993 // and number of trials from pyxsec.root
994 // Copied from AliAnalysisTaskJetSpectrum2
995 //
996
997 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
998 Float_t xsection = 0;
999 Float_t ftrials = 1;
1000
1001 fAvgTrials = 1;
1002 if(tree){
1003 TFile *curfile = tree->GetCurrentFile();
1004 if (!curfile) {
1005 Error("Notify","No current file");
1006 return kFALSE;
1007 }
1008 if(!fh1Xsec||!fh1Trials){
1009 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1010 return kFALSE;
1011 }
1012 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1013 fh1Xsec->Fill("<#sigma>",xsection);
1014 // construct a poor man average trials
1015 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1016 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1017 }
1018 return kTRUE;
1019}
1020
1021//________________________________________________________________________
1022AliGenPythiaEventHeader* AliPWG4HighPtTrackQA::GetPythiaEventHeader(AliMCEvent *mcEvent){
1023
1024 if(!mcEvent)return 0;
1025 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1026 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1027 if(!pythiaGenHeader){
1028 // cocktail ??
1029 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1030
1031 if (!genCocktailHeader) {
1032 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1033 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1034 return 0;
1035 }
1036 TList* headerList = genCocktailHeader->GetHeaders();
1037 for (Int_t i=0; i<headerList->GetEntries(); i++) {
1038 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1039 if (pythiaGenHeader)
1040 break;
1041 }
1042 if(!pythiaGenHeader){
1043 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1044 return 0;
1045 }
1046 }
1047 return pythiaGenHeader;
1048
1049}
1050
1051//_______________________________________________________________________
1052Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1053{
1054 //MV: copied from AliESDtrack since method is not available in AliAODTrack
1055
1056 //
1057 // TPC cluster information
1058 // type 0: get fraction of found/findable clusters with neighbourhood definition
1059 // 1: findable clusters with neighbourhood definition
1060 // 2: found clusters
1061 //
1062 // definition of findable clusters:
1063 // a cluster is defined as findable if there is another cluster
1064 // within +- nNeighbours pad rows. The idea is to overcome threshold
1065 // effects with a very simple algorithm.
1066 //
1067
1068 TBits fTPCClusterMap = tr->GetTPCClusterMap();
1069 if (type==2) return fTPCClusterMap.CountBits();
1070
1071 Int_t found=0;
1072 Int_t findable=0;
1073 Int_t last=-nNeighbours;
1074
1075 for (Int_t i=row0; i<row1; ++i){
1076 //look to current row
1077 if (fTPCClusterMap[i]) {
1078 last=i;
1079 ++found;
1080 ++findable;
1081 continue;
1082 }
1083 //look to nNeighbours before
1084 if ((i-last)<=nNeighbours) {
1085 ++findable;
1086 continue;
1087 }
1088 //look to nNeighbours after
1089 for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1090 if (fTPCClusterMap[j]){
1091 ++findable;
1092 break;
1093 }
1094 }
1095 }
1096 if (type==1) return findable;
1097
1098 if (type==0){
1099 Float_t fraction=0;
1100 if (findable>0)
1101 fraction=(Float_t)found/(Float_t)findable;
1102 else
1103 fraction=0;
1104 return fraction;
1105 }
1106 return 0; // undefined type - default value
1107}
1108
1109
1110//________________________________________________________________________
1111void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1112{
1113 // The Terminate() function is the last function to be called during
1114 // a query. It always runs on the client, it can be used to present
1115 // the results graphically or save the results to file.
1116
1117}
1118
1119#endif