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