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