Updated QA task for constraint global tracks (M. Verweij)
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliPWG4HighPtTrackQA.cxx
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
57 using namespace std; //required for resolving the 'cout' symbol
58
59 ClassImp(AliPWG4HighPtTrackQA)
60
61 AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA()
62 : AliAnalysisTaskSE(), 
63   fDataType(kESD),
64   fEvent(0x0),
65   fESD(0x0),
66   fVtx(0x0),
67   fTrackCuts(0x0), 
68   fTrackCutsITSLoose(0x0), 
69   fTrackCutsTPConly(0x0), 
70   fTrackType(0),
71   fFilterMask(0),
72   fSigmaConstrainedMax(-1.),
73   fPtMax(100.),
74   fIsPbPb(0),
75   fCentClass(10),
76   fNVariables(18),
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),
101   fPtUncertainty1Pt(0x0),
102   fPtChi2PerClusterTPC(0x0),
103   fPtNCrossedRows(0x0),
104   fPtNCrossedRowsNClusF(0x0),
105   fPtNCrRNCrRNClusF(0x0),
106   fPtSigmaY2(0x0),
107   fPtSigmaZ2(0x0),
108   fPtSigmaSnp2(0x0),
109   fPtSigmaTgl2(0x0),
110   fPtSigma1Pt2(0x0),
111   fProfPtSigmaY2(0x0),
112   fProfPtSigmaZ2(0x0),
113   fProfPtSigmaSnp2(0x0),
114   fProfPtSigmaTgl2(0x0),
115   fProfPtSigma1Pt2(0x0),
116   fProfPtSigma1Pt(0x0),
117   fProfPtPtSigma1Pt(0x0),
118   fSystTrackCuts(0x0),
119   fHistList(0)
120 {
121   SetNVariables(18);
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
130 }
131 //________________________________________________________________________
132 AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA(const char *name): 
133   AliAnalysisTaskSE(name), 
134   fDataType(kESD),
135   fEvent(0x0),
136   fESD(0x0),
137   fVtx(0x0),
138   fTrackCuts(0x0),
139   fTrackCutsITSLoose(0x0), 
140   fTrackCutsTPConly(0x0), 
141   fTrackType(0),
142   fFilterMask(0),
143   fSigmaConstrainedMax(-1.),
144   fPtMax(100.),
145   fIsPbPb(0),
146   fCentClass(10),
147   fNVariables(18),
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),
172   fPtUncertainty1Pt(0x0),
173   fPtChi2PerClusterTPC(0x0),
174   fPtNCrossedRows(0x0),
175   fPtNCrossedRowsNClusF(0x0),
176   fPtNCrRNCrRNClusF(0x0),
177   fPtSigmaY2(0x0),
178   fPtSigmaZ2(0x0),
179   fPtSigmaSnp2(0x0),
180   fPtSigmaTgl2(0x0),
181   fPtSigma1Pt2(0x0),
182   fProfPtSigmaY2(0x0),
183   fProfPtSigmaZ2(0x0),
184   fProfPtSigmaSnp2(0x0),
185   fProfPtSigmaTgl2(0x0),
186   fProfPtSigma1Pt2(0x0),
187   fProfPtSigma1Pt(0x0),
188   fProfPtPtSigma1Pt(0x0),
189   fSystTrackCuts(0x0),
190   fHistList(0)
191 {
192   //
193   // Constructor. Initialization of Inputs and Outputs
194   //
195   AliDebug(2,Form("AliPWG4HighPtTrackQA Calling Constructor"));
196
197   SetNVariables(18);
198
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
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
212 //________________________________________________________________________
213 void 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
226 //________________________________________________________________________
227 void 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
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
256   const Float_t ptmin1 =  fgkPtMin;
257   const Float_t ptmax1 =  fPtBinEdges[0][0];
258   const Float_t ptmin2 =  ptmax1 ;
259   const Float_t ptmax2 =  fPtBinEdges[1][0];
260   const Float_t ptmin3 =  ptmax2 ;
261   const Float_t ptmax3 =  fgkPtMax;
262   const Int_t nbin11 = (int)((ptmax1-ptmin1)/fPtBinEdges[0][1]);
263   const Int_t nbin12 = (int)((ptmax2-ptmin2)/fPtBinEdges[1][1])+nbin11;
264   const Int_t nbin13 = (int)((ptmax3-ptmin3)/fPtBinEdges[2][1])+nbin12;
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;
295   if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
296     fgkDCA2DMin = -2.;
297     fgkDCA2DMax = 2.;
298   }
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.;
305   if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
306     fgkDCAZMin = -5.;
307     fgkDCAZMax = 5.;
308   }
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.;
326   Float_t fgkChi2CMax = 100.;
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;
333   if(fTrackType==1 || fTrackType==2 || fTrackType==4) 
334     fgkRel1PtUncertaintyMax = 0.5; 
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
338   Int_t fgkNUncertainty1PtBins = 30;
339   Float_t fgkUncertainty1PtMin = 0.;
340   Float_t fgkUncertainty1PtMax = 0.1;
341   if(fTrackType==1 || fTrackType==2 || fTrackType==4) 
342     fgkUncertainty1PtMax = 0.2; 
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
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
358   Float_t fgk1PtMin = 0.;
359   Float_t fgk1PtMax = 6.;
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;
366   Double_t *bins1Pt=new Double_t[fgkN1PtBins+1];
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   }
374
375   Int_t fgkNSigmaY2Bins = 50;
376   Float_t fgkSigmaY2Min = 0.;
377   Float_t fgkSigmaY2Max = 1.;
378   if(fTrackType==1) fgkSigmaY2Max = 4.;
379   if(fTrackType==2 || fTrackType==4) fgkSigmaY2Max = 0.1;
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.;
385   Float_t fgkSigmaZ2Max = 0.4;
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.;
391   Float_t fgkSigmaSnp2Max = 0.05;
392   if(fTrackType==1) fgkSigmaSnp2Max = 0.2;
393   if(fTrackType==2 || fTrackType==4) fgkSigmaSnp2Max = 0.1;
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.;
399   Float_t fgkSigmaTgl2Max = 0.1;
400   if(fTrackType==1) fgkSigmaTgl2Max = 0.2;
401   if(fTrackType==2 || fTrackType==4) fgkSigmaTgl2Max = 0.1;
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.;
407   Float_t fgkSigma1Pt2Max = 1.;
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
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);
456   fh1NTracksReject->Fill("chi2",0);
457   fHistList->Add(fh1NTracksReject);
458
459   fh1NTracksSel = new TH1F("fh1NTracksSel","fh1NTracksSel",1,-0.5,0.5);
460   fHistList->Add(fh1NTracksSel);
461
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
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);
504
505   fPtUncertainty1Pt = new TH2F("fPtUncertainty1Pt","fPtUncertainty1Pt",fgkNPtBins,binsPt,fgkNUncertainty1PtBins,binsUncertainty1Pt);
506   fHistList->Add(fPtUncertainty1Pt);
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);
519
520   fPtSigmaY2 = new TH2F("fPtSigmaY2","fPtSigmaY2",fgkN1PtBins,bins1Pt,fgkNSigmaY2Bins,binsSigmaY2);
521   fHistList->Add(fPtSigmaY2);
522
523   fPtSigmaZ2 = new TH2F("fPtSigmaZ2","fPtSigmaZ2",fgkN1PtBins,bins1Pt,fgkNSigmaZ2Bins,binsSigmaZ2);
524   fHistList->Add(fPtSigmaZ2);
525
526   fPtSigmaSnp2 = new TH2F("fPtSigmaSnp2","fPtSigmaSnp2",fgkN1PtBins,bins1Pt,fgkNSigmaSnp2Bins,binsSigmaSnp2);
527   fHistList->Add(fPtSigmaSnp2);
528
529   fPtSigmaTgl2 = new TH2F("fPtSigmaTgl2","fPtSigmaTgl2",fgkN1PtBins,bins1Pt,fgkNSigmaTgl2Bins,binsSigmaTgl2);
530   fHistList->Add(fPtSigmaTgl2);
531
532   fPtSigma1Pt2 = new TH2F("fPtSigma1Pt2","fPtSigma1Pt2",fgkN1PtBins,bins1Pt,fgkNSigma1Pt2Bins,binsSigma1Pt2);
533   fHistList->Add(fPtSigma1Pt2);
534
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
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;
570   if(binsUncertainty1Pt)    delete [] binsUncertainty1Pt;
571   if(binsChi2PerClus)       delete [] binsChi2PerClus;
572   if(binsChi2PerClus)       delete [] binsNCrossedRowsNClusF;
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;
579 }
580
581 //________________________________________________________________________
582 Bool_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
614   if(fDataType==kESD&&dynamic_cast<AliESDEvent*>(fEvent)) {
615     fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexSPD();
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   }
645   else if(fDataType==kAOD&&dynamic_cast<AliAODEvent*>(fEvent)) {
646     const AliAODVertex *vtx = ((AliAODEvent*)fEvent)->GetPrimaryVertexSPD();
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 //________________________________________________________________________
703 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliVEvent *ev){
704   //
705   // Get centrality from ESD or AOD
706   //
707
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 //________________________________________________________________________
717 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliESDEvent *esd){
718   //
719   // Get centrality from ESD
720   //
721
722   Float_t cent = -1;
723
724   if(esd){
725     if(esd->GetCentrality()){
726       cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
727       if(fDebug>3) printf("centrality: %f\n",cent);
728     }
729   }
730
731   return GetCentralityClass(cent);
732
733 }
734
735 //________________________________________________________________________
736 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliAODEvent *aod){
737   //
738   // Get centrality from AOD
739   //
740
741   if(!aod) return 5;
742   Float_t cent = aod->GetHeader()->GetCentrality();
743   if(fDebug>3) printf("centrality: %f\n",cent);
744
745   return GetCentralityClass(cent);
746
747 }
748
749 //________________________________________________________________________
750 Int_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;
760   return 0;
761
762 }
763
764 //________________________________________________________________________
765 void 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 //________________________________________________________________________
800 void 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)
846     13: SigmaY2
847     14: SigmaZ2
848     15: SigmaSnp2
849     16: SigmaTgl2
850     17: Sigma1Pt2
851   */
852
853   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
854     fh1NTracksAll->Fill(0.);
855
856     //Get track for analysis
857     AliESDtrack *track = 0x0;
858     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
859     if(!esdtrack) {
860       fh1NTracksReject->Fill("noESDtrack",1);
861       continue;
862     }
863
864     if(fTrackType==4) {
865       FillSystematicCutHist(esdtrack);
866       if (!(fTrackCuts->AcceptTrack(esdtrack))) {
867         fh1NTracksReject->Fill("trackCuts",1);
868         continue;
869       }
870     }
871
872     if(fTrackType==1)
873       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
874     else if(fTrackType==2 || fTrackType==4) {
875       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
876       if(!track) {
877         fh1NTracksReject->Fill("noTPConly",1);
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);
884         delete track;
885         continue;
886       }
887       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
888     }
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());
899             AliExternalTrackParam exParam;
900             Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
901             if( !relate ) {
902               fh1NTracksReject->Fill("relate",1);
903               delete track;
904               continue;
905             }
906             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
907           }
908           else if(fTrackType==6) {
909             //use global constrained track
910             track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
911
912           }
913         }
914       }
915     }
916     else
917       track = esdtrack;
918     
919     if(!track) {
920       //      if(fTrackType==1 || fTrackType==2 || fTrackType==4) delete track;
921       continue;
922     }
923
924     if(fTrackType==2 || fTrackType==4) {
925       //Cut on chi2 of constrained fit
926       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
927         fh1NTracksReject->Fill("chi2",1);
928         delete track;
929         continue;
930       }
931     }
932
933     if(fTrackType!=4)
934       FillSystematicCutHist(track);
935
936     fPtAll->Fill(track->Pt());
937
938     if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
939       fh1NTracksReject->Fill("trackCuts",1);
940       if(fTrackType==1 || fTrackType==2) delete track;
941       continue;
942     }
943
944     fh1NTracksSel->Fill(0.);
945
946     fVariables->Reset(0.);
947       
948     fVariables->SetAt(track->Pt(),0); 
949     fVariables->SetAt(track->Phi(),1); 
950     fVariables->SetAt(track->Eta(),2); 
951
952     Float_t dca2D = 0.;
953     Float_t dcaz  = 0.;
954     if(fTrackType==0) {       //Global
955       track->GetImpactParameters(dca2D,dcaz);
956     }
957     else if(fTrackType==1 || fTrackType==2 || fTrackType==4) {  //TPConly
958       track->GetImpactParametersTPC(dca2D,dcaz);
959     }
960     fVariables->SetAt(dca2D,3);
961     fVariables->SetAt(dcaz,5);
962
963     fVariables->SetAt((float)track->GetTPCNcls(),5);
964
965     Int_t nPointITS = 0;
966     UChar_t itsMap = track->GetITSClusterMap();
967     for (Int_t i=0; i < 6; i++) {
968       if (itsMap & (1 << i))
969         nPointITS ++;
970     }
971     fVariables->SetAt((float)nPointITS,6);
972     Float_t chi2C = (float)track->GetConstrainedChi2();
973     if(fTrackType==1 || fTrackType==2 || fTrackType==4)
974       chi2C = (float)track->GetConstrainedChi2TPC();
975     fVariables->SetAt(chi2C,7);
976     fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
977   
978     fVariables->SetAt(TMath::Sqrt(track->GetSigma1Pt2())*fVariables->At(0),9);
979   
980     if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
981     
982     fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
983     Float_t crossedRowsTPCNClsF = 1.;//track->GetTPCClusterInfo(2,0);
984     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
985     fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
986     fVariables->SetAt(track->GetSigmaY2(),13);
987     fVariables->SetAt(track->GetSigmaZ2(),14);
988     fVariables->SetAt(track->GetSigmaSnp2(),15);
989     fVariables->SetAt(track->GetSigmaTgl2(),16);
990     fVariables->SetAt(track->GetSigma1Pt2(),17);
991     
992     FillHistograms();
993   
994     //      int mult = fTrackCuts->CountAcceptedTracks(fESD);
995
996     if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5) delete track;
997     
998   }//track loop
999
1000 }
1001
1002 //________________________________________________________________________
1003 void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
1004
1005   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1006   if(!aod)return;
1007   for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
1008
1009     AliAODTrack *aodtrack = aod->GetTrack(iTrack);
1010     if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
1011
1012     fVariables->Reset(0.);
1013
1014     fVariables->SetAt(aodtrack->Pt(),0);
1015     fVariables->SetAt(aodtrack->Phi(),1);
1016     fVariables->SetAt(aodtrack->Eta(),2);
1017
1018     Double_t dca[2] = {1e6,1e6};
1019     Double_t covar[3] = {1e6,1e6,1e6};
1020     if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
1021       fVariables->SetAt(dca[0],3);
1022       fVariables->SetAt(dca[1],4);
1023     }
1024
1025     fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
1026     fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
1027     fVariables->SetAt(aodtrack->Chi2perNDF(),7);
1028     fVariables->SetAt(0.,8);
1029     fVariables->SetAt(0.,9);
1030     fVariables->SetAt(aodtrack->Chi2perNDF(),10);
1031     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2),11);
1032     Float_t crossedRowsTPCNClsF = 0.;
1033     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
1034     fVariables->SetAt(crossedRowsTPCNClsF,12);
1035     
1036     fPtAll->Fill(fVariables->At(0));
1037
1038     FillHistograms();
1039
1040   }
1041
1042 }
1043
1044 //________________________________________________________________________
1045 void AliPWG4HighPtTrackQA::FillHistograms() {
1046
1047   fPtSel->Fill(fVariables->At(0));
1048   fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
1049   fPtEta->Fill(fVariables->At(0),fVariables->At(2));
1050   fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
1051   fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
1052   fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
1053   fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
1054   if(fDataType==kESD) {
1055     fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
1056     fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
1057     fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(9));
1058     fPtUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)/fVariables->At(9));
1059     fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1060     fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1061     fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1062     fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1063     fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1064
1065     fProfPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1066     fProfPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1067     fProfPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1068     fProfPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1069     fProfPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1070     fProfPtSigma1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1071     fProfPtPtSigma1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1072   }
1073   fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
1074   fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
1075   fPtNCrossedRowsNClusF->Fill(fVariables->At(0),fVariables->At(12));
1076   fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
1077 }
1078
1079 //________________________________________________________________________
1080 Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1081   //
1082   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1083   // This is to called in Notify and should provide the path to the AOD/ESD file
1084   // Copied from AliAnalysisTaskJetSpectrum2
1085   //
1086
1087   TString file(currFile);  
1088   fXsec = 0;
1089   fTrials = 1;
1090
1091   if(file.Contains("root_archive.zip#")){
1092     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1093     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1094     file.Replace(pos+1,20,"");
1095   }
1096   else {
1097     // not an archive take the basename....
1098     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1099   }
1100   //  Printf("%s",file.Data());
1101    
1102
1103   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
1104   if(!fxsec){
1105     // next trial fetch the histgram file
1106     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1107     if(!fxsec){
1108         // not a severe condition but inciate that we have no information
1109       return kFALSE;
1110     }
1111     else{
1112       // find the tlist we want to be independtent of the name so use the Tkey
1113       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
1114       if(!key){
1115         fxsec->Close();
1116         return kFALSE;
1117       }
1118       TList *list = dynamic_cast<TList*>(key->ReadObj());
1119       if(!list){
1120         fxsec->Close();
1121         return kFALSE;
1122       }
1123       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1124       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1125       fxsec->Close();
1126     }
1127   } // no tree pyxsec.root
1128   else {
1129     TTree *xtree = (TTree*)fxsec->Get("Xsection");
1130     if(!xtree){
1131       fxsec->Close();
1132       return kFALSE;
1133     }
1134     UInt_t   ntrials  = 0;
1135     Double_t  xsection  = 0;
1136     xtree->SetBranchAddress("xsection",&xsection);
1137     xtree->SetBranchAddress("ntrials",&ntrials);
1138     xtree->GetEntry(0);
1139     fTrials = ntrials;
1140     fXsec = xsection;
1141     fxsec->Close();
1142   }
1143   return kTRUE;
1144 }
1145 //________________________________________________________________________
1146 Bool_t AliPWG4HighPtTrackQA::Notify()
1147 {
1148   //
1149   // Implemented Notify() to read the cross sections
1150   // and number of trials from pyxsec.root
1151   // Copied from AliAnalysisTaskJetSpectrum2
1152   // 
1153
1154   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1155   Float_t xsection = 0;
1156   Float_t ftrials  = 1;
1157
1158   fAvgTrials = 1;
1159   if(tree){
1160     TFile *curfile = tree->GetCurrentFile();
1161     if (!curfile) {
1162       Error("Notify","No current file");
1163       return kFALSE;
1164     }
1165     if(!fh1Xsec||!fh1Trials){
1166       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1167       return kFALSE;
1168     }
1169     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1170     fh1Xsec->Fill("<#sigma>",xsection);
1171     // construct a poor man average trials 
1172     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1173     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1174   }  
1175   return kTRUE;
1176 }
1177
1178 //________________________________________________________________________
1179 AliGenPythiaEventHeader*  AliPWG4HighPtTrackQA::GetPythiaEventHeader(AliMCEvent *mcEvent){
1180   
1181   if(!mcEvent)return 0;
1182   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1183   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1184   if(!pythiaGenHeader){
1185     // cocktail ??
1186     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1187     
1188     if (!genCocktailHeader) {
1189       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1190       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1191       return 0;
1192     }
1193     TList* headerList = genCocktailHeader->GetHeaders();
1194     for (Int_t i=0; i<headerList->GetEntries(); i++) {
1195       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1196       if (pythiaGenHeader)
1197         break;
1198     }
1199     if(!pythiaGenHeader){
1200       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1201       return 0;
1202     }
1203   }
1204   return pythiaGenHeader;
1205
1206 }
1207
1208 //_______________________________________________________________________
1209 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1210 {
1211   //MV: copied from AliESDtrack since method is not available in AliAODTrack
1212
1213   //
1214   // TPC cluster information
1215   // type 0: get fraction of found/findable clusters with neighbourhood definition
1216   //      1: findable clusters with neighbourhood definition
1217   //      2: found clusters
1218   //
1219   // definition of findable clusters:
1220   //            a cluster is defined as findable if there is another cluster
1221   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1222   //           effects with a very simple algorithm.
1223   //
1224
1225   TBits fTPCClusterMap = tr->GetTPCClusterMap(); 
1226   if (type==2) return fTPCClusterMap.CountBits();
1227
1228   Int_t found=0;
1229   Int_t findable=0;
1230   Int_t last=-nNeighbours;
1231
1232   for (Int_t i=row0; i<row1; ++i){
1233     //look to current row
1234     if (fTPCClusterMap[i]) {
1235       last=i;
1236       ++found;
1237       ++findable;
1238       continue;
1239     }
1240     //look to nNeighbours before
1241     if ((i-last)<=nNeighbours) {
1242       ++findable;
1243       continue;
1244     }
1245     //look to nNeighbours after
1246     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1247       if (fTPCClusterMap[j]){
1248         ++findable;
1249         break;
1250       }
1251     }
1252   }
1253   if (type==1) return findable;
1254
1255   if (type==0){
1256     Float_t fraction=0;
1257     if (findable>0)
1258       fraction=(Float_t)found/(Float_t)findable;
1259     else
1260       fraction=0;
1261     return fraction;
1262   }
1263   return 0;  // undefined type - default value
1264 }
1265
1266 void AliPWG4HighPtTrackQA::FillSystematicCutHist(AliESDtrack *track) {
1267
1268   fSystTrackCuts->Fill("noCut",1);
1269   if(TMath::Abs(track->Eta())>0.9) 
1270     fSystTrackCuts->Fill("eta",1);
1271   if(track->Pt()<0.15 || track->Pt()>1e10) 
1272     fSystTrackCuts->Fill("0.15<pT<1e10",1);
1273   if(track->GetKinkIndex(0)>0)
1274     fSystTrackCuts->Fill("kink",1);
1275   if(track->GetTPCclusters(0)<70)
1276     fSystTrackCuts->Fill("NClusterTPC",1);
1277   if(track->GetTPCclusters(0)>0.) {
1278     if(track->GetTPCchi2()/Float_t(track->GetTPCclusters(0))>4.)
1279       fSystTrackCuts->Fill("Chi2PerNClusTPC",1);
1280   }
1281   
1282   Float_t dcaToVertexXY = 0.;
1283   Float_t dcaToVertexZ  = 0.;
1284   track->GetImpactParameters(dcaToVertexXY,dcaToVertexZ);
1285   Float_t fCutMaxDCAToVertexXY = 2.4;
1286   Float_t fCutMaxDCAToVertexZ  = 3.2;
1287   Float_t dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1288   if(dcaToVertex>1.) 
1289     fSystTrackCuts->Fill("DCA2D",1);
1290   
1291 }
1292
1293 //________________________________________________________________________
1294 void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1295 {
1296   // The Terminate() function is the last function to be called during
1297   // a query. It always runs on the client, it can be used to present
1298   // the results graphically or save the results to file.
1299
1300 }
1301
1302 #endif