ca1cfbd75e7572b891fb2c8d75757a9e5daf4504
[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         if(track) 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             if(!track) {
900               fh1NTracksReject->Fill("noTPConly",1);
901               continue;
902             }
903             AliExternalTrackParam exParam;
904             Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
905             if( !relate ) {
906               fh1NTracksReject->Fill("relate",1);
907               if(track) delete track;
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
914             track = esdtrack;
915             track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
916
917           }
918         }
919       }
920     }
921     else
922       track = esdtrack;
923     
924     if(!track) {
925       continue;
926     }
927
928     if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
929       //Cut on chi2 of constrained fit
930       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
931         fh1NTracksReject->Fill("chi2",1);
932         if(track) delete track;
933         continue;
934       }
935     }
936
937     if(fTrackType!=4)
938       FillSystematicCutHist(track);
939
940     fPtAll->Fill(track->Pt());
941
942     if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
943       fh1NTracksReject->Fill("trackCuts",1);
944       if(fTrackType==1 || fTrackType==2) {
945         if(track) delete track;
946       }
947       continue;
948     }
949
950     fh1NTracksSel->Fill(0.);
951
952     fVariables->Reset(0.);
953       
954     fVariables->SetAt(track->Pt(),0); 
955     fVariables->SetAt(track->Phi(),1); 
956     fVariables->SetAt(track->Eta(),2); 
957
958     Float_t dca2D = 0.;
959     Float_t dcaz  = 0.;
960     if(fTrackType==0) {       //Global
961       track->GetImpactParameters(dca2D,dcaz);
962     }
963     else if(fTrackType==1 || fTrackType==2 || fTrackType==4) {  //TPConly
964       track->GetImpactParametersTPC(dca2D,dcaz);
965     }
966     fVariables->SetAt(dca2D,3);
967     fVariables->SetAt(dcaz,5);
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);
978     Float_t chi2C = (float)track->GetConstrainedChi2();
979     if(fTrackType==1 || fTrackType==2 || fTrackType==4)
980       chi2C = (float)track->GetConstrainedChi2TPC();
981     fVariables->SetAt(chi2C,7);
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     
988     fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
989     Float_t crossedRowsTPCNClsF = 1.;//track->GetTPCClusterInfo(2,0);
990     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
991     fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
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);
997     
998     FillHistograms();
999   
1000     //      int mult = fTrackCuts->CountAcceptedTracks(fESD);
1001
1002     if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5) {
1003       if(track) delete track;
1004     }
1005     
1006   }//track loop
1007
1008 }
1009
1010 //________________________________________________________________________
1011 void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
1012
1013   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1014   if(!aod)return;
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);
1035     fVariables->SetAt(aodtrack->Chi2perNDF(),7);
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 //________________________________________________________________________
1053 void 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));
1066     fPtUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)*fVariables->At(9));
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)));
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)));
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 //________________________________________________________________________
1088 Bool_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 //________________________________________________________________________
1154 Bool_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 //________________________________________________________________________
1187 AliGenPythiaEventHeader*  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 //_______________________________________________________________________
1217 Float_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
1274 void 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 }
1300
1301 //________________________________________________________________________
1302 void 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