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