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