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