switchable response matrix
[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){
700 if(fDataType==kESD)
701 return CalculateCentrality(dynamic_cast<AliESDEvent*>(ev));
702 else if(fDataType==kAOD)
703 return CalculateCentrality(dynamic_cast<AliAODEvent*>(ev));
704 else
705 return 5;
706}
707
708//________________________________________________________________________
709Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliESDEvent *esd){
710
1ea145ef 711
d756027f 712 Float_t cent = 999;
1ea145ef 713
714 if(esd){
715 if(esd->GetCentrality()){
716 cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
717 if(fDebug>3) cout << "centrality: " << cent << endl;
718 }
d756027f 719 }
720
721 if(cent>80)return 4;
722 if(cent>50)return 3;
723 if(cent>30)return 2;
724 if(cent>10)return 1;
725 return 0;
726
727}
728
729//________________________________________________________________________
730Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliAODEvent *aod){
731
1ea145ef 732 if(!aod)return 4;
d756027f 733 Float_t cent = aod->GetHeader()->GetCentrality();
734 cout << "centrality: " << cent << endl;
735 if(cent>80)return 4;
736 if(cent>50)return 3;
737 if(cent>30)return 2;
738 if(cent>10)return 1;
739 return 0;
740
741}
742
743//________________________________________________________________________
744void AliPWG4HighPtTrackQA::UserExec(Option_t *) {
745 // Main loop
746 // Called for each event
747 AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserExec \n"));
748
749 fEvent = InputEvent();
750 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
751
752 // All events without selection
753 fNEventAll->Fill(0.);
754
755 if(!SelectEvent()) {
756 // Post output data
757 PostData(1, fHistList);
758 return;
759 }
760
761
762 //Need to keep track of selected events
763 fNEventSel->Fill(0.);
764
765 fVariables = new TArrayF(fNVariables);
766
767 if(fDataType==kESD) DoAnalysisESD();
768 if(fDataType==kAOD) DoAnalysisAOD();
769
770 //Delete old fVariables
771 if(fVariables) delete fVariables;
772
773 // Post output data
774 PostData(1, fHistList);
775
776}
777
778//________________________________________________________________________
779void AliPWG4HighPtTrackQA::DoAnalysisESD() {
780
781 if(!fESD) {
782 PostData(1, fHistList);
783 return;
784 }
785
786 // ---- Get MC Header information (for MC productions in pThard bins) ----
787 Double_t ptHard = 0.;
788 Double_t nTrials = 1; // trials for MC trigger weight for real data
789
790 AliMCEventHandler *eventHandlerMC = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
791 if (eventHandlerMC) {
792
793 if(eventHandlerMC->MCEvent()){
794 AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(eventHandlerMC->MCEvent());
795 if(pythiaGenHeader){
796 nTrials = pythiaGenHeader->Trials();
797 ptHard = pythiaGenHeader->GetPtHard();
798
799 fh1PtHard->Fill(ptHard);
800 fh1PtHardTrials->Fill(ptHard,nTrials);
801
802 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
803 }
804 }
805 }
806
807 Int_t nTracks = fESD->GetNumberOfTracks();
808 AliDebug(2,Form("nTracks ESD%d", nTracks));
809
810 /*
811 Variables to be put in fVariables
812 0: pt
813 1: phi
814 2: eta
815 3: dca2D
816 4: dcaZ
817 5: nClustersTPC
818 6: nPointITS
819 7: chi2C
820 8: nSigmaToVertex
821 9: relUncertainty1Pt
822 10: chi2PerClusterTPC
823 11: #crossed rows
824 12: (#crossed rows)/(#findable clusters)
2b553e6f 825 13: SigmaY2
826 14: SigmaZ2
827 15: SigmaSnp2
828 16: SigmaTgl2
829 17: Sigma1Pt2
d756027f 830 */
831
832 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
833 fh1NTracksAll->Fill(0.);
834
835 //Get track for analysis
2b553e6f 836 AliESDtrack *track = 0x0;
d756027f 837 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
838 if(!esdtrack) {
839 fh1NTracksReject->Fill("noESDtrack",1);
840 continue;
841 }
842
aa3ba8d2 843 if(fTrackType==4) {
844 FillSystematicCutHist(esdtrack);
845 if (!(fTrackCuts->AcceptTrack(esdtrack))) {
846 fh1NTracksReject->Fill("trackCuts",1);
847 continue;
848 }
849 }
850
d756027f 851 if(fTrackType==1)
852 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
aa3ba8d2 853 else if(fTrackType==2 || fTrackType==4) {
d756027f 854 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
855 if(!track) {
856 fh1NTracksReject->Fill("noTPConly",1);
d756027f 857 continue;
858 }
859 AliExternalTrackParam exParam;
860 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
861 if( !relate ) {
862 fh1NTracksReject->Fill("relate",1);
863 delete track;
864 continue;
865 }
866 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
867 }
868 else
869 track = esdtrack;
870
871 if(!track) {
afc75e88 872 // if(fTrackType==1 || fTrackType==2 || fTrackType==4) delete track;
d756027f 873 continue;
874 }
875
aa3ba8d2 876 if(fTrackType==2 || fTrackType==4) {
2b553e6f 877 //Cut on chi2 of constrained fit
aa3ba8d2 878 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
2b553e6f 879 fh1NTracksReject->Fill("chi2",1);
880 delete track;
881 continue;
882 }
883 }
884
aa3ba8d2 885 if(fTrackType!=4)
886 FillSystematicCutHist(track);
887
d756027f 888 fPtAll->Fill(track->Pt());
889
aa3ba8d2 890 if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4) {
d756027f 891 fh1NTracksReject->Fill("trackCuts",1);
892 if(fTrackType==1 || fTrackType==2) delete track;
893 continue;
894 }
aa3ba8d2 895
896 //Cut out laser tracks
d756027f 897 if(track->GetTPCsignal()<10) { //Cut on laser tracks
898 fh1NTracksReject->Fill("laser",1);
aa3ba8d2 899 if(fTrackType==1 || fTrackType==2 || fTrackType==4) delete track;
d756027f 900 continue;
901 }
902
903 fh1NTracksSel->Fill(0.);
904
905 fVariables->Reset(0.);
906
2b553e6f 907 fVariables->SetAt(track->Pt(),0);
908 fVariables->SetAt(track->Phi(),1);
909 fVariables->SetAt(track->Eta(),2);
d756027f 910
911 Float_t dca2D = 0.;
912 Float_t dcaz = 0.;
913 if(fTrackType==0) { //Global
914 track->GetImpactParameters(dca2D,dcaz);
915 }
aa3ba8d2 916 else if(fTrackType==1 || fTrackType==2 || fTrackType==4) { //TPConly
d756027f 917 track->GetImpactParametersTPC(dca2D,dcaz);
918 }
919 fVariables->SetAt(dca2D,3);
aa3ba8d2 920 fVariables->SetAt(dcaz,5);
d756027f 921
922 fVariables->SetAt((float)track->GetTPCNcls(),5);
923
924 Int_t nPointITS = 0;
925 UChar_t itsMap = track->GetITSClusterMap();
926 for (Int_t i=0; i < 6; i++) {
927 if (itsMap & (1 << i))
928 nPointITS ++;
929 }
930 fVariables->SetAt((float)nPointITS,6);
2b553e6f 931 Float_t chi2C = (float)track->GetConstrainedChi2();
aa3ba8d2 932 if(fTrackType==1 || fTrackType==2 || fTrackType==4)
2b553e6f 933 chi2C = (float)track->GetConstrainedChi2TPC();
934 fVariables->SetAt(chi2C,7);
d756027f 935 fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
936
937 fVariables->SetAt(TMath::Sqrt(track->GetSigma1Pt2())*fVariables->At(0),9);
938
939 if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
940
941 //cout << "#crossed rows (1): " << track->GetTPCClusterInfo(1) << endl;
942 //cout << "#crossed rows (2): " << track->GetTPCClusterInfo(2) << endl;
943 fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
944 Float_t crossedRowsTPCNClsF = track->GetTPCClusterInfo(2,0);
945 //if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
946 fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
2b553e6f 947 fVariables->SetAt(track->GetSigmaY2(),13);
948 fVariables->SetAt(track->GetSigmaZ2(),14);
949 fVariables->SetAt(track->GetSigmaSnp2(),15);
950 fVariables->SetAt(track->GetSigmaTgl2(),16);
951 fVariables->SetAt(track->GetSigma1Pt2(),17);
d756027f 952
953 FillHistograms();
954
955 // int mult = fTrackCuts->CountAcceptedTracks(fESD);
956
aa3ba8d2 957 if(fTrackType==1 || fTrackType==2 || fTrackType==4) delete track;
d756027f 958
959 }//track loop
960
961}
962
963//________________________________________________________________________
964void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
965
966 AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1ea145ef 967 if(!aod)return;
d756027f 968 for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
969
970 AliAODTrack *aodtrack = aod->GetTrack(iTrack);
971 if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
972
973 fVariables->Reset(0.);
974
975 fVariables->SetAt(aodtrack->Pt(),0);
976 fVariables->SetAt(aodtrack->Phi(),1);
977 fVariables->SetAt(aodtrack->Eta(),2);
978
979 Double_t dca[2] = {1e6,1e6};
980 Double_t covar[3] = {1e6,1e6,1e6};
981 if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
982 fVariables->SetAt(dca[0],3);
983 fVariables->SetAt(dca[1],4);
984 }
985
986 fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
987 fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
aa3ba8d2 988 fVariables->SetAt(aodtrack->Chi2perNDF(),7);
d756027f 989 fVariables->SetAt(0.,8);
990 fVariables->SetAt(0.,9);
991 fVariables->SetAt(aodtrack->Chi2perNDF(),10);
992 fVariables->SetAt(GetTPCClusterInfo(aodtrack,2),11);
993 Float_t crossedRowsTPCNClsF = 0.;
994 if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
995 fVariables->SetAt(crossedRowsTPCNClsF,12);
996
997 fPtAll->Fill(fVariables->At(0));
998
999 FillHistograms();
1000
1001 }
1002
1003}
1004
1005//________________________________________________________________________
1006void AliPWG4HighPtTrackQA::FillHistograms() {
1007
1008 fPtSel->Fill(fVariables->At(0));
1009 fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
1010 fPtEta->Fill(fVariables->At(0),fVariables->At(2));
1011 fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
1012 fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
1013 fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
1014 fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
1015 if(fDataType==kESD) {
1016 fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
1017 fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
1018 fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(9));
2b553e6f 1019 fPtUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)/fVariables->At(9));
1020 fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1021 fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1022 fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1023 fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1024 fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
aa3ba8d2 1025
1026 fProfPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1027 fProfPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1028 fProfPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1029 fProfPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1030 fProfPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1031 fProfPtSigma1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1032 fProfPtPtSigma1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
d756027f 1033 }
1034 fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
1035 fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
1036 fPtNCrossedRowsNClusF->Fill(fVariables->At(0),fVariables->At(12));
1037 fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
1038}
1039
1040//________________________________________________________________________
1041Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1042 //
1043 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1044 // This is to called in Notify and should provide the path to the AOD/ESD file
1045 // Copied from AliAnalysisTaskJetSpectrum2
1046 //
1047
1048 TString file(currFile);
1049 fXsec = 0;
1050 fTrials = 1;
1051
1052 if(file.Contains("root_archive.zip#")){
1053 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1054 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1055 file.Replace(pos+1,20,"");
1056 }
1057 else {
1058 // not an archive take the basename....
1059 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1060 }
1061 // Printf("%s",file.Data());
1062
1063
1064 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
1065 if(!fxsec){
1066 // next trial fetch the histgram file
1067 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1068 if(!fxsec){
1069 // not a severe condition but inciate that we have no information
1070 return kFALSE;
1071 }
1072 else{
1073 // find the tlist we want to be independtent of the name so use the Tkey
1074 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1075 if(!key){
1076 fxsec->Close();
1077 return kFALSE;
1078 }
1079 TList *list = dynamic_cast<TList*>(key->ReadObj());
1080 if(!list){
1081 fxsec->Close();
1082 return kFALSE;
1083 }
1084 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1085 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1086 fxsec->Close();
1087 }
1088 } // no tree pyxsec.root
1089 else {
1090 TTree *xtree = (TTree*)fxsec->Get("Xsection");
1091 if(!xtree){
1092 fxsec->Close();
1093 return kFALSE;
1094 }
1095 UInt_t ntrials = 0;
1096 Double_t xsection = 0;
1097 xtree->SetBranchAddress("xsection",&xsection);
1098 xtree->SetBranchAddress("ntrials",&ntrials);
1099 xtree->GetEntry(0);
1100 fTrials = ntrials;
1101 fXsec = xsection;
1102 fxsec->Close();
1103 }
1104 return kTRUE;
1105}
1106//________________________________________________________________________
1107Bool_t AliPWG4HighPtTrackQA::Notify()
1108{
1109 //
1110 // Implemented Notify() to read the cross sections
1111 // and number of trials from pyxsec.root
1112 // Copied from AliAnalysisTaskJetSpectrum2
1113 //
1114
1115 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1116 Float_t xsection = 0;
1117 Float_t ftrials = 1;
1118
1119 fAvgTrials = 1;
1120 if(tree){
1121 TFile *curfile = tree->GetCurrentFile();
1122 if (!curfile) {
1123 Error("Notify","No current file");
1124 return kFALSE;
1125 }
1126 if(!fh1Xsec||!fh1Trials){
1127 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1128 return kFALSE;
1129 }
1130 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1131 fh1Xsec->Fill("<#sigma>",xsection);
1132 // construct a poor man average trials
1133 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1134 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1135 }
1136 return kTRUE;
1137}
1138
1139//________________________________________________________________________
1140AliGenPythiaEventHeader* AliPWG4HighPtTrackQA::GetPythiaEventHeader(AliMCEvent *mcEvent){
1141
1142 if(!mcEvent)return 0;
1143 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1144 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1145 if(!pythiaGenHeader){
1146 // cocktail ??
1147 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1148
1149 if (!genCocktailHeader) {
1150 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1151 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1152 return 0;
1153 }
1154 TList* headerList = genCocktailHeader->GetHeaders();
1155 for (Int_t i=0; i<headerList->GetEntries(); i++) {
1156 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1157 if (pythiaGenHeader)
1158 break;
1159 }
1160 if(!pythiaGenHeader){
1161 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1162 return 0;
1163 }
1164 }
1165 return pythiaGenHeader;
1166
1167}
1168
1169//_______________________________________________________________________
1170Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1171{
1172 //MV: copied from AliESDtrack since method is not available in AliAODTrack
1173
1174 //
1175 // TPC cluster information
1176 // type 0: get fraction of found/findable clusters with neighbourhood definition
1177 // 1: findable clusters with neighbourhood definition
1178 // 2: found clusters
1179 //
1180 // definition of findable clusters:
1181 // a cluster is defined as findable if there is another cluster
1182 // within +- nNeighbours pad rows. The idea is to overcome threshold
1183 // effects with a very simple algorithm.
1184 //
1185
1186 TBits fTPCClusterMap = tr->GetTPCClusterMap();
1187 if (type==2) return fTPCClusterMap.CountBits();
1188
1189 Int_t found=0;
1190 Int_t findable=0;
1191 Int_t last=-nNeighbours;
1192
1193 for (Int_t i=row0; i<row1; ++i){
1194 //look to current row
1195 if (fTPCClusterMap[i]) {
1196 last=i;
1197 ++found;
1198 ++findable;
1199 continue;
1200 }
1201 //look to nNeighbours before
1202 if ((i-last)<=nNeighbours) {
1203 ++findable;
1204 continue;
1205 }
1206 //look to nNeighbours after
1207 for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1208 if (fTPCClusterMap[j]){
1209 ++findable;
1210 break;
1211 }
1212 }
1213 }
1214 if (type==1) return findable;
1215
1216 if (type==0){
1217 Float_t fraction=0;
1218 if (findable>0)
1219 fraction=(Float_t)found/(Float_t)findable;
1220 else
1221 fraction=0;
1222 return fraction;
1223 }
1224 return 0; // undefined type - default value
1225}
1226
aa3ba8d2 1227void AliPWG4HighPtTrackQA::FillSystematicCutHist(AliESDtrack *track) {
1228
1229 fSystTrackCuts->Fill("noCut",1);
1230 if(TMath::Abs(track->Eta())>0.9)
1231 fSystTrackCuts->Fill("eta",1);
1232 if(track->Pt()<0.15 || track->Pt()>1e10)
1233 fSystTrackCuts->Fill("0.15<pT<1e10",1);
1234 if(track->GetKinkIndex(0)>0)
1235 fSystTrackCuts->Fill("kink",1);
1236 if(track->GetTPCclusters(0)<70)
1237 fSystTrackCuts->Fill("NClusterTPC",1);
1238 if(track->GetTPCclusters(0)>0.) {
1239 if(track->GetTPCchi2()/Float_t(track->GetTPCclusters(0))>4.)
1240 fSystTrackCuts->Fill("Chi2PerNClusTPC",1);
1241 }
1242
1243 Float_t dcaToVertexXY = 0.;
1244 Float_t dcaToVertexZ = 0.;
1245 track->GetImpactParameters(dcaToVertexXY,dcaToVertexZ);
1246 Float_t fCutMaxDCAToVertexXY = 2.4;
1247 Float_t fCutMaxDCAToVertexZ = 3.2;
1248 Float_t dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1249 if(dcaToVertex>1.)
1250 fSystTrackCuts->Fill("DCA2D",1);
1251
1252}
d756027f 1253
1254//________________________________________________________________________
1255void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1256{
1257 // The Terminate() function is the last function to be called during
1258 // a query. It always runs on the client, it can be used to present
1259 // the results graphically or save the results to file.
1260
1261}
1262
1263#endif