cfedea758a3e28d8c578b2409eea98ea952127cc
[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   fPtMax(100.),
71   fIsPbPb(0),
72   fCentClass(10),
73   fNVariables(13),
74   fVariables(0x0),
75   fAvgTrials(1),
76   fNEventAll(0),
77   fNEventSel(0),
78   fNEventReject(0),
79   fh1Centrality(0x0),
80   fh1Xsec(0),
81   fh1Trials(0),
82   fh1PtHard(0),
83   fh1PtHardTrials(0),
84   fh1NTracksAll(0x0),
85   fh1NTracksReject(0x0),
86   fh1NTracksSel(0x0),
87   fPtAll(0),  
88   fPtSel(0),    
89   fPtPhi(0x0),
90   fPtEta(0x0),
91   fPtDCA2D(0x0),
92   fPtDCAZ(0x0),
93   fPtNClustersTPC(0x0),
94   fPtNPointITS(0x0),
95   fPtChi2C(0x0),
96   fPtNSigmaToVertex(0x0),
97   fPtRelUncertainty1Pt(0x0),
98   fPtChi2PerClusterTPC(0x0),
99   fPtNCrossedRows(0x0),
100   fPtNCrossedRowsNClusF(0x0),
101   fPtNCrRNCrRNClusF(0x0),
102   fHistList(0)
103 {
104   SetNVariables(13);
105 }
106 //________________________________________________________________________
107 AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA(const char *name): 
108   AliAnalysisTaskSE(name), 
109   fDataType(kESD),
110   fEvent(0x0),
111   fESD(0x0),
112   fVtx(0x0),
113   fTrackCuts(),
114   fTrackType(0),
115   fFilterMask(0),
116   fPtMax(100.),
117   fIsPbPb(0),
118   fCentClass(10),
119   fNVariables(13),
120   fVariables(0x0),
121   fAvgTrials(1),
122   fNEventAll(0),
123   fNEventSel(0),
124   fNEventReject(0),
125   fh1Centrality(0x0),
126   fh1Xsec(0),
127   fh1Trials(0),
128   fh1PtHard(0),
129   fh1PtHardTrials(0),
130   fh1NTracksAll(0x0),
131   fh1NTracksReject(0x0),
132   fh1NTracksSel(0x0),
133   fPtAll(0),
134   fPtSel(0),
135   fPtPhi(0x0),
136   fPtEta(0x0),
137   fPtDCA2D(0x0),
138   fPtDCAZ(0x0),
139   fPtNClustersTPC(0x0),
140   fPtNPointITS(0x0),
141   fPtChi2C(0x0),
142   fPtNSigmaToVertex(0x0),
143   fPtRelUncertainty1Pt(0x0),
144   fPtChi2PerClusterTPC(0x0),
145   fPtNCrossedRows(0x0),
146   fPtNCrossedRowsNClusF(0x0),
147   fPtNCrRNCrRNClusF(0x0),
148   fHistList(0)
149 {
150   //
151   // Constructor. Initialization of Inputs and Outputs
152   //
153   AliDebug(2,Form("AliPWG4HighPtTrackQA Calling Constructor"));
154
155   SetNVariables(13);
156
157   // Input slot #0 works with a TChain ESD
158   DefineInput(0, TChain::Class());
159   // Output slot #1 write into a TList
160   DefineOutput(1, TList::Class());
161 }
162
163 //________________________________________________________________________
164 void AliPWG4HighPtTrackQA::UserCreateOutputObjects() {
165   //Create output objects
166   AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserCreateOutputObjects \n"));
167
168   Bool_t oldStatus = TH1::AddDirectoryStatus();
169   TH1::AddDirectory(kFALSE); 
170
171   OpenFile(1);
172   fHistList = new TList();
173   fHistList->SetOwner(kTRUE);
174   
175   Float_t fgkPtMin = 0.;
176   Float_t fgkPtMax = fPtMax;
177
178   Float_t ptBinEdges[2][2];
179   ptBinEdges[0][0] = 10.;
180   ptBinEdges[0][1] = 1.;
181   ptBinEdges[1][0] = 20.;
182   ptBinEdges[1][1] = 2.;
183   Float_t binWidth3 = 5.;
184   if(fPtMax>100.) {
185     ptBinEdges[0][0] = 100.;
186     ptBinEdges[0][1] = 5.;
187     ptBinEdges[1][0] = 300.;
188     ptBinEdges[1][1] = 10.;
189     binWidth3 = 20.;
190   }
191
192   const Float_t ptmin1 =  fgkPtMin;
193   const Float_t ptmax1 =  ptBinEdges[0][0];
194   const Float_t ptmin2 =  ptmax1 ;
195   const Float_t ptmax2 =  ptBinEdges[1][0];
196   const Float_t ptmin3 =  ptmax2 ;
197   const Float_t ptmax3 =  fgkPtMax;
198   const Int_t nbin11 = (int)((ptmax1-ptmin1)/ptBinEdges[0][1]);
199   const Int_t nbin12 = (int)((ptmax2-ptmin2)/ptBinEdges[1][1])+nbin11;
200   const Int_t nbin13 = (int)((ptmax3-ptmin3)/binWidth3)+nbin12;
201   Int_t fgkNPtBins=nbin13;
202   //Create array with low edges of each bin
203   Double_t *binsPt=new Double_t[fgkNPtBins+1];
204   for(Int_t i=0; i<=fgkNPtBins; i++) {
205     if(i<=nbin11) binsPt[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
206     if(i<=nbin12 && i>nbin11) binsPt[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
207     if(i<=nbin13 && i>nbin12) binsPt[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
208   }
209
210   Int_t fgkNPhiBins = 18*6;
211   Float_t kMinPhi   = 0.;
212   Float_t kMaxPhi   = 2.*TMath::Pi();
213   Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
214   for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
215
216   Int_t fgkNEtaBins=20;
217   Float_t fgkEtaMin = -1.;
218   Float_t fgkEtaMax = 1.;
219   Double_t *binsEta=new Double_t[fgkNEtaBins+1];
220   for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
221
222   Int_t fgkNNClustersTPCBins=80;
223   Float_t fgkNClustersTPCMin = 0.5;
224   Float_t fgkNClustersTPCMax = 160.5;
225   Double_t *binsNClustersTPC=new Double_t[fgkNNClustersTPCBins+1];
226   for(Int_t i=0; i<=fgkNNClustersTPCBins; i++) binsNClustersTPC[i]=(Double_t)fgkNClustersTPCMin + (fgkNClustersTPCMax-fgkNClustersTPCMin)/fgkNNClustersTPCBins*(Double_t)i ;
227
228   Int_t fgkNDCA2DBins=80;
229   Float_t fgkDCA2DMin = -0.2;
230   Float_t fgkDCA2DMax = 0.2;
231   Double_t *binsDCA2D=new Double_t[fgkNDCA2DBins+1];
232   for(Int_t i=0; i<=fgkNDCA2DBins; i++) binsDCA2D[i]=(Double_t)fgkDCA2DMin + (fgkDCA2DMax-fgkDCA2DMin)/fgkNDCA2DBins*(Double_t)i ;
233
234   Int_t fgkNDCAZBins=80;
235   Float_t fgkDCAZMin = -2.;
236   Float_t fgkDCAZMax = 2.;
237   Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
238   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
239
240   Int_t fgkNNPointITSBins=9;
241   Float_t fgkNPointITSMin = -0.5;
242   Float_t fgkNPointITSMax = 8.5;
243   Double_t *binsNPointITS=new Double_t[fgkNNPointITSBins+1];
244   for(Int_t i=0; i<=fgkNNPointITSBins; i++) binsNPointITS[i]=(Double_t)fgkNPointITSMin + (fgkNPointITSMax-fgkNPointITSMin)/fgkNNPointITSBins*(Double_t)i ;
245
246   Int_t fgkNNSigmaToVertexBins=40;
247   Float_t fgkNSigmaToVertexMin = 0.;
248   Float_t fgkNSigmaToVertexMax = 8.;
249   Double_t *binsNSigmaToVertex=new Double_t[fgkNNSigmaToVertexBins+1];
250   for(Int_t i=0; i<=fgkNNSigmaToVertexBins; i++) binsNSigmaToVertex[i]=(Double_t)fgkNSigmaToVertexMin + (fgkNSigmaToVertexMax-fgkNSigmaToVertexMin)/fgkNNSigmaToVertexBins*(Double_t)i ;
251
252   Int_t fgkNChi2CBins=20;
253   Float_t fgkChi2CMin = 0.;
254   Float_t fgkChi2CMax = 10.;
255   Double_t *binsChi2C=new Double_t[fgkNChi2CBins+1];
256   for(Int_t i=0; i<=fgkNChi2CBins; i++) binsChi2C[i]=(Double_t)fgkChi2CMin + (fgkChi2CMax-fgkChi2CMin)/fgkNChi2CBins*(Double_t)i ;
257
258   Int_t fgkNRel1PtUncertaintyBins=30;
259   Float_t fgkRel1PtUncertaintyMin = 0.;
260   Float_t fgkRel1PtUncertaintyMax = 0.3;
261   Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
262   for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
263
264   Float_t fgkChi2PerClusMin = 0.;
265   Float_t fgkChi2PerClusMax = 4.;
266   Int_t fgkNChi2PerClusBins = (int)(fgkChi2PerClusMax*10.);
267   Double_t *binsChi2PerClus=new Double_t[fgkNChi2PerClusBins+1];
268   for(Int_t i=0; i<=fgkNChi2PerClusBins; i++) binsChi2PerClus[i]=(Double_t)fgkChi2PerClusMin + (fgkChi2PerClusMax-fgkChi2PerClusMin)/fgkNChi2PerClusBins*(Double_t)i ;
269
270   Int_t fgkNCrossedRowsNClusFBins  = 50;
271   Float_t fgkNCrossedRowsNClusFMin = 0.;
272   Float_t fgkNCrossedRowsNClusFMax = 1.;
273   Double_t *binsNCrossedRowsNClusF=new Double_t[fgkNCrossedRowsNClusFBins+1];
274   for(Int_t i=0; i<=fgkNCrossedRowsNClusFBins; i++) binsNCrossedRowsNClusF[i]=(Double_t)fgkNCrossedRowsNClusFMin + (fgkNCrossedRowsNClusFMax-fgkNCrossedRowsNClusFMin)/fgkNCrossedRowsNClusFBins*(Double_t)i ;
275
276   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
277   fHistList->Add(fNEventAll);
278   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
279   fHistList->Add(fNEventSel);
280   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
281   //Set labels
282   fNEventReject->Fill("noESD",0);
283   fNEventReject->Fill("Trigger",0);
284   fNEventReject->Fill("NTracks<2",0);
285   fNEventReject->Fill("noVTX",0);
286   fNEventReject->Fill("VtxStatus",0);
287   fNEventReject->Fill("NCont<2",0);
288   fNEventReject->Fill("ZVTX>10",0);
289   fNEventReject->Fill("cent",0);
290   fNEventReject->Fill("cent>90",0);
291   fHistList->Add(fNEventReject);
292
293   fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
294   fHistList->Add(fh1Centrality);
295
296   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
297   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
298   fHistList->Add(fh1Xsec);
299
300   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
301   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
302   fHistList->Add(fh1Trials);
303
304   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
305   fHistList->Add(fh1PtHard);
306   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
307   fHistList->Add(fh1PtHardTrials);
308
309   fh1NTracksAll = new TH1F("fh1NTracksAll","fh1NTracksAll",1,-0.5,0.5);
310   fHistList->Add(fh1NTracksAll);
311
312   fh1NTracksReject = new TH1F("fh1NTracksReject","fh1NTracksReject",1,-0.5,0.5);
313   fh1NTracksReject->Fill("noESDtrack",0);
314   fh1NTracksReject->Fill("noTPCInner",0);
315   fh1NTracksReject->Fill("FillTPC",0);
316   fh1NTracksReject->Fill("noTPConly",0);
317   fh1NTracksReject->Fill("relate",0);
318   fh1NTracksReject->Fill("trackCuts",0);
319   fh1NTracksReject->Fill("laser",0);
320   fHistList->Add(fh1NTracksReject);
321
322   fh1NTracksSel = new TH1F("fh1NTracksSel","fh1NTracksSel",1,-0.5,0.5);
323   fHistList->Add(fh1NTracksSel);
324
325
326   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
327   fHistList->Add(fPtAll);
328   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
329   fHistList->Add(fPtSel);
330
331   fPtPhi = new TH2F("fPtPhi","fPtPhi",fgkNPtBins,binsPt,fgkNPhiBins,binsPhi);
332   fHistList->Add(fPtPhi);
333  
334   fPtEta = new TH2F("fPtEta","fPtEta",fgkNPtBins,binsPt,fgkNEtaBins,binsEta);
335   fHistList->Add(fPtEta);
336  
337   fPtDCA2D = new TH2F("fPtDCA2D","fPtDCA2D",fgkNPtBins,binsPt,fgkNDCA2DBins,binsDCA2D);
338   fHistList->Add(fPtDCA2D);
339  
340   fPtDCAZ = new TH2F("fPtDCAZ","fPtDCAZ",fgkNPtBins,binsPt,fgkNDCAZBins,binsDCAZ);
341   fHistList->Add(fPtDCAZ);
342  
343   fPtNClustersTPC = new TH2F("fPtNClustersTPC","fPtNClustersTPC",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
344   fHistList->Add(fPtNClustersTPC);
345  
346   fPtNPointITS = new TH2F("fPtNPointITS","fPtNPointITS",fgkNPtBins,binsPt,fgkNNPointITSBins,binsNPointITS);
347   fHistList->Add(fPtNPointITS);
348  
349   fPtChi2C = new TH2F("fPtChi2C","fPtChi2C",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
350   fHistList->Add(fPtChi2C);
351  
352   fPtNSigmaToVertex = new TH2F("fPtNSigmaToVertex","fPtNSigmaToVertex",fgkNPtBins,binsPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
353   fHistList->Add(fPtNSigmaToVertex);
354
355   fPtRelUncertainty1Pt = new TH2F("fPtRelUncertainty1Pt","fPtRelUncertainty1Pt",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
356   fHistList->Add(fPtRelUncertainty1Pt);
357  
358   fPtChi2PerClusterTPC = new TH2F("fPtChi2PerClusterTPC","fPtChi2PerClusterTPC",fgkNPtBins,binsPt,fgkNChi2PerClusBins,binsChi2PerClus);
359   fHistList->Add(fPtChi2PerClusterTPC);
360
361   fPtNCrossedRows = new TH2F("fPtNCrossedRows","fPtNCrossedRows",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
362   fHistList->Add(fPtNCrossedRows);
363
364   fPtNCrossedRowsNClusF = new TH2F("fPtNCrossedRowsNClusF","fPtNCrossedRowsNClusF",fgkNPtBins,binsPt,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
365   fHistList->Add(fPtNCrossedRowsNClusF);
366  
367   fPtNCrRNCrRNClusF = new TH3F("fPtNCrRNCrRNClusF","fPtNCrRNCrRNClusF",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
368   fHistList->Add(fPtNCrRNCrRNClusF);
369  
370   TH1::AddDirectory(oldStatus); 
371
372   PostData(1, fHistList);
373
374   if(binsPhi)               delete [] binsPhi;
375   if(binsPt)                delete [] binsPt;
376   if(binsNClustersTPC)      delete [] binsNClustersTPC;
377   if(binsDCA2D)             delete [] binsDCA2D;
378   if(binsDCAZ)              delete [] binsDCAZ;
379   if(binsNPointITS)         delete [] binsNPointITS;
380   if(binsNSigmaToVertex)    delete [] binsNSigmaToVertex;
381   if(binsChi2C)             delete [] binsChi2C;
382   if(binsEta)               delete [] binsEta;
383   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
384   if(binsChi2PerClus)       delete [] binsChi2PerClus;
385   if(binsChi2PerClus)       delete [] binsNCrossedRowsNClusF;
386 }
387
388 //________________________________________________________________________
389 Bool_t AliPWG4HighPtTrackQA::SelectEvent() {
390   //
391   // Decide if event should be selected for analysis
392   //
393
394   // Checks following requirements:
395   // - fEvent available
396   // - trigger info from AliPhysicsSelection
397   // - MCevent available
398   // - number of reconstructed tracks > 1
399   // - primary vertex reconstructed
400   // - z-vertex < 10 cm
401   // - centrality in case of PbPb
402
403   Bool_t selectEvent = kTRUE;
404
405   //fEvent object available?
406   if (!fEvent) {
407     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
408     fNEventReject->Fill("noAliVEvent",1);
409     selectEvent = kFALSE;
410     return selectEvent;
411   }
412
413   //Check if number of reconstructed tracks is larger than 1
414   if(!fEvent->GetNumberOfTracks() || fEvent->GetNumberOfTracks()<2)  {
415     fNEventReject->Fill("NTracks<2",1);
416     selectEvent = kFALSE;
417     return selectEvent;
418   }
419
420   //Check if vertex is reconstructed
421   if(fDataType==kESD&&dynamic_cast<AliESDEvent*>(fEvent)) {
422     fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexSPD();
423
424     if(!fVtx) {
425       fNEventReject->Fill("noVTX",1);
426       selectEvent = kFALSE;
427       return selectEvent;
428     }
429     
430     if(!fVtx->GetStatus()) {
431       fNEventReject->Fill("VtxStatus",1);
432       selectEvent = kFALSE;
433       return selectEvent;
434     }
435     
436     // Need vertex cut
437     if(fVtx->GetNContributors()<2) {
438       fNEventReject->Fill("NCont<2",1);
439       selectEvent = kFALSE;
440       return selectEvent;
441     }
442     
443     //Check if z-vertex < 10 cm
444     double primVtx[3];
445     fVtx->GetXYZ(primVtx);
446     if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
447       fNEventReject->Fill("ZVTX>10",1);
448       selectEvent = kFALSE;
449       return selectEvent;
450     }
451   }
452   else if(fDataType==kAOD&&dynamic_cast<AliAODEvent*>(fEvent)) {
453     const AliAODVertex *vtx = ((AliAODEvent*)fEvent)->GetPrimaryVertexSPD();
454     if(!vtx) {
455       fNEventReject->Fill("noVTX",1);
456       selectEvent = kFALSE;
457       return selectEvent;
458     }
459     
460     // Need vertex cut
461     if(vtx->GetNContributors()<2) {
462       fNEventReject->Fill("NCont<2",1);
463       selectEvent = kFALSE;
464       return selectEvent;
465     }
466     
467     //Check if z-vertex < 10 cm
468     double primVtx[3];
469     vtx->GetXYZ(primVtx);
470     if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
471       fNEventReject->Fill("ZVTX>10",1);
472       selectEvent = kFALSE;
473       return selectEvent;
474     }
475
476   }
477
478   //Centrality selection should only be done in case of PbPb
479   if(IsPbPb()) {
480     Float_t cent = 0.;
481     if(fCentClass!=CalculateCentrality(fEvent) && fCentClass!=10) {
482       fNEventReject->Fill("cent",1);
483       selectEvent = kFALSE;
484       return selectEvent;
485     }
486     else {
487       if(fDataType==kESD) {
488         if(dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()) {
489           cent = dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()->GetCentralityPercentile("V0M");
490         }
491       }
492       else if(fDataType==kAOD) {
493         if(dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality())
494           cent = dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality();
495        }
496       if(cent>90.) {
497         fNEventReject->Fill("cent>90",1);
498         selectEvent = kFALSE;
499         return selectEvent;     
500       }
501       fh1Centrality->Fill(cent);
502     }
503   }
504  
505   return selectEvent;
506
507 }
508
509 //________________________________________________________________________
510 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliVEvent *ev){
511   if(fDataType==kESD)
512     return CalculateCentrality(dynamic_cast<AliESDEvent*>(ev));
513   else if(fDataType==kAOD)
514     return CalculateCentrality(dynamic_cast<AliAODEvent*>(ev));
515   else
516     return 5;
517 }
518
519 //________________________________________________________________________
520 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliESDEvent *esd){
521
522
523   Float_t cent = 999;
524
525   if(esd){
526     if(esd->GetCentrality()){
527       cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
528       if(fDebug>3) cout << "centrality: " << cent << endl;
529     }
530   }
531
532   if(cent>80)return 4;
533   if(cent>50)return 3;
534   if(cent>30)return 2;
535   if(cent>10)return 1;
536   return 0;
537
538 }
539
540 //________________________________________________________________________
541 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliAODEvent *aod){
542
543   if(!aod)return 4;
544   Float_t cent = aod->GetHeader()->GetCentrality();
545   cout << "centrality: " << cent << endl;
546   if(cent>80)return 4;
547   if(cent>50)return 3;
548   if(cent>30)return 2;
549   if(cent>10)return 1;
550   return 0;
551
552 }
553
554 //________________________________________________________________________
555 void AliPWG4HighPtTrackQA::UserExec(Option_t *) {  
556   // Main loop
557   // Called for each event
558   AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserExec \n"));
559   
560   fEvent = InputEvent();
561   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
562
563   // All events without selection
564   fNEventAll->Fill(0.);
565
566   if(!SelectEvent()) {
567     // Post output data
568     PostData(1, fHistList);
569     return;
570   }
571
572
573   //Need to keep track of selected events
574   fNEventSel->Fill(0.);
575
576   fVariables = new TArrayF(fNVariables);
577   
578   if(fDataType==kESD) DoAnalysisESD();
579   if(fDataType==kAOD) DoAnalysisAOD();
580
581   //Delete old fVariables
582   if(fVariables) delete fVariables;
583
584   // Post output data
585   PostData(1, fHistList);
586
587 }
588
589 //________________________________________________________________________
590 void AliPWG4HighPtTrackQA::DoAnalysisESD() {
591
592   if(!fESD) {
593     PostData(1, fHistList);
594     return;
595   }
596
597   // ---- Get MC Header information (for MC productions in pThard bins) ----
598   Double_t ptHard = 0.;
599   Double_t nTrials = 1; // trials for MC trigger weight for real data
600   
601   AliMCEventHandler *eventHandlerMC = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
602   if (eventHandlerMC) {
603     
604     if(eventHandlerMC->MCEvent()){
605       AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(eventHandlerMC->MCEvent());
606       if(pythiaGenHeader){
607         nTrials = pythiaGenHeader->Trials();
608         ptHard  = pythiaGenHeader->GetPtHard();
609         
610         fh1PtHard->Fill(ptHard);
611         fh1PtHardTrials->Fill(ptHard,nTrials);
612         
613         fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
614       }
615     }
616   }
617
618   Int_t nTracks = fESD->GetNumberOfTracks();
619   AliDebug(2,Form("nTracks ESD%d", nTracks));
620
621   /*
622     Variables to be put in fVariables
623     0: pt
624     1: phi
625     2: eta
626     3: dca2D
627     4: dcaZ 
628     5: nClustersTPC
629     6: nPointITS   
630     7: chi2C       
631     8: nSigmaToVertex
632     9: relUncertainty1Pt
633     10: chi2PerClusterTPC
634     11: #crossed rows
635     12: (#crossed rows)/(#findable clusters)
636   */
637
638   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
639     fh1NTracksAll->Fill(0.);
640
641     //Get track for analysis
642     AliESDtrack *track;
643     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
644     if(!esdtrack) {
645       fh1NTracksReject->Fill("noESDtrack",1);
646       continue;
647     }
648
649     if(fTrackType==1)
650       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
651     else if(fTrackType==2) {
652       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
653       if(!track) {
654         fh1NTracksReject->Fill("noTPConly",1);
655         continue;
656       }
657       AliExternalTrackParam exParam;
658       Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
659       if( !relate ) {
660         fh1NTracksReject->Fill("relate",1);
661         delete track;
662         continue;
663       }
664       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
665     }
666     else
667       track = esdtrack;
668     
669     if(!track) {
670       if(fTrackType==1 || fTrackType==2) delete track;
671       continue;
672     }
673
674     fPtAll->Fill(track->Pt());
675
676     if (!(fTrackCuts->AcceptTrack(track))) {
677       fh1NTracksReject->Fill("trackCuts",1);
678       if(fTrackType==1 || fTrackType==2) delete track;
679       continue;
680     }
681     if(track->GetTPCsignal()<10) { //Cut on laser tracks
682       fh1NTracksReject->Fill("laser",1);
683       if(fTrackType==1 || fTrackType==2) delete track;
684       continue; 
685     }
686
687     fh1NTracksSel->Fill(0.);
688
689     fVariables->Reset(0.);
690       
691     fVariables->SetAt(track->Pt(),0);
692     fVariables->SetAt(track->Phi(),1);
693     fVariables->SetAt(track->Eta(),2);
694
695     Float_t dca2D = 0.;
696     Float_t dcaz  = 0.;
697     if(fTrackType==0) {       //Global
698       track->GetImpactParameters(dca2D,dcaz);
699     }
700     else if(fTrackType==1 || fTrackType==2) {  //TPConly
701       track->GetImpactParametersTPC(dca2D,dcaz);
702     }
703     fVariables->SetAt(dca2D,3);
704     fVariables->SetAt(dcaz,5);
705
706     fVariables->SetAt((float)track->GetTPCNcls(),5);
707
708     Int_t nPointITS = 0;
709     UChar_t itsMap = track->GetITSClusterMap();
710     for (Int_t i=0; i < 6; i++) {
711       if (itsMap & (1 << i))
712         nPointITS ++;
713     }
714     fVariables->SetAt((float)nPointITS,6);
715     fVariables->SetAt(track->GetConstrainedChi2(),7);
716     fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
717   
718     fVariables->SetAt(TMath::Sqrt(track->GetSigma1Pt2())*fVariables->At(0),9);
719   
720     if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
721     
722     //cout << "#crossed rows (1): " << track->GetTPCClusterInfo(1) << endl;
723     //cout << "#crossed rows (2): " << track->GetTPCClusterInfo(2) << endl;
724     fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
725     Float_t crossedRowsTPCNClsF = track->GetTPCClusterInfo(2,0);
726     //if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
727     fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
728     
729     FillHistograms();
730   
731     //      int mult = fTrackCuts->CountAcceptedTracks(fESD);
732
733     if(fTrackType==1  || fTrackType==2) delete track;
734     
735   }//track loop
736
737 }
738
739 //________________________________________________________________________
740 void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
741
742   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
743   if(!aod)return;
744   for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
745
746     AliAODTrack *aodtrack = aod->GetTrack(iTrack);
747     if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
748
749     fVariables->Reset(0.);
750
751     fVariables->SetAt(aodtrack->Pt(),0);
752     fVariables->SetAt(aodtrack->Phi(),1);
753     fVariables->SetAt(aodtrack->Eta(),2);
754
755     Double_t dca[2] = {1e6,1e6};
756     Double_t covar[3] = {1e6,1e6,1e6};
757     if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
758       fVariables->SetAt(dca[0],3);
759       fVariables->SetAt(dca[1],4);
760     }
761
762     fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
763     fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
764     fVariables->SetAt(0.,7);
765     fVariables->SetAt(0.,8);
766     fVariables->SetAt(0.,9);
767     fVariables->SetAt(aodtrack->Chi2perNDF(),10);
768     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2),11);
769     Float_t crossedRowsTPCNClsF = 0.;
770     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
771     fVariables->SetAt(crossedRowsTPCNClsF,12);
772     
773     fPtAll->Fill(fVariables->At(0));
774
775     FillHistograms();
776
777   }
778
779 }
780
781 //________________________________________________________________________
782 void AliPWG4HighPtTrackQA::FillHistograms() {
783
784   fPtSel->Fill(fVariables->At(0));
785   fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
786   fPtEta->Fill(fVariables->At(0),fVariables->At(2));
787   fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
788   fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
789   fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
790   fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
791   if(fDataType==kESD) {
792     fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
793     fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
794     fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(9));
795   }
796   fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
797   fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
798   fPtNCrossedRowsNClusF->Fill(fVariables->At(0),fVariables->At(12));
799   fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
800 }
801
802 //________________________________________________________________________
803 Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
804   //
805   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
806   // This is to called in Notify and should provide the path to the AOD/ESD file
807   // Copied from AliAnalysisTaskJetSpectrum2
808   //
809
810   TString file(currFile);  
811   fXsec = 0;
812   fTrials = 1;
813
814   if(file.Contains("root_archive.zip#")){
815     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
816     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
817     file.Replace(pos+1,20,"");
818   }
819   else {
820     // not an archive take the basename....
821     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
822   }
823   //  Printf("%s",file.Data());
824    
825
826   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
827   if(!fxsec){
828     // next trial fetch the histgram file
829     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
830     if(!fxsec){
831         // not a severe condition but inciate that we have no information
832       return kFALSE;
833     }
834     else{
835       // find the tlist we want to be independtent of the name so use the Tkey
836       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
837       if(!key){
838         fxsec->Close();
839         return kFALSE;
840       }
841       TList *list = dynamic_cast<TList*>(key->ReadObj());
842       if(!list){
843         fxsec->Close();
844         return kFALSE;
845       }
846       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
847       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
848       fxsec->Close();
849     }
850   } // no tree pyxsec.root
851   else {
852     TTree *xtree = (TTree*)fxsec->Get("Xsection");
853     if(!xtree){
854       fxsec->Close();
855       return kFALSE;
856     }
857     UInt_t   ntrials  = 0;
858     Double_t  xsection  = 0;
859     xtree->SetBranchAddress("xsection",&xsection);
860     xtree->SetBranchAddress("ntrials",&ntrials);
861     xtree->GetEntry(0);
862     fTrials = ntrials;
863     fXsec = xsection;
864     fxsec->Close();
865   }
866   return kTRUE;
867 }
868 //________________________________________________________________________
869 Bool_t AliPWG4HighPtTrackQA::Notify()
870 {
871   //
872   // Implemented Notify() to read the cross sections
873   // and number of trials from pyxsec.root
874   // Copied from AliAnalysisTaskJetSpectrum2
875   // 
876
877   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
878   Float_t xsection = 0;
879   Float_t ftrials  = 1;
880
881   fAvgTrials = 1;
882   if(tree){
883     TFile *curfile = tree->GetCurrentFile();
884     if (!curfile) {
885       Error("Notify","No current file");
886       return kFALSE;
887     }
888     if(!fh1Xsec||!fh1Trials){
889       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
890       return kFALSE;
891     }
892     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
893     fh1Xsec->Fill("<#sigma>",xsection);
894     // construct a poor man average trials 
895     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
896     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
897   }  
898   return kTRUE;
899 }
900
901 //________________________________________________________________________
902 AliGenPythiaEventHeader*  AliPWG4HighPtTrackQA::GetPythiaEventHeader(AliMCEvent *mcEvent){
903   
904   if(!mcEvent)return 0;
905   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
906   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
907   if(!pythiaGenHeader){
908     // cocktail ??
909     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
910     
911     if (!genCocktailHeader) {
912       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
913       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
914       return 0;
915     }
916     TList* headerList = genCocktailHeader->GetHeaders();
917     for (Int_t i=0; i<headerList->GetEntries(); i++) {
918       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
919       if (pythiaGenHeader)
920         break;
921     }
922     if(!pythiaGenHeader){
923       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
924       return 0;
925     }
926   }
927   return pythiaGenHeader;
928
929 }
930
931 //_______________________________________________________________________
932 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
933 {
934   //MV: copied from AliESDtrack since method is not available in AliAODTrack
935
936   //
937   // TPC cluster information
938   // type 0: get fraction of found/findable clusters with neighbourhood definition
939   //      1: findable clusters with neighbourhood definition
940   //      2: found clusters
941   //
942   // definition of findable clusters:
943   //            a cluster is defined as findable if there is another cluster
944   //           within +- nNeighbours pad rows. The idea is to overcome threshold
945   //           effects with a very simple algorithm.
946   //
947
948   TBits fTPCClusterMap = tr->GetTPCClusterMap(); 
949   if (type==2) return fTPCClusterMap.CountBits();
950
951   Int_t found=0;
952   Int_t findable=0;
953   Int_t last=-nNeighbours;
954
955   for (Int_t i=row0; i<row1; ++i){
956     //look to current row
957     if (fTPCClusterMap[i]) {
958       last=i;
959       ++found;
960       ++findable;
961       continue;
962     }
963     //look to nNeighbours before
964     if ((i-last)<=nNeighbours) {
965       ++findable;
966       continue;
967     }
968     //look to nNeighbours after
969     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
970       if (fTPCClusterMap[j]){
971         ++findable;
972         break;
973       }
974     }
975   }
976   if (type==1) return findable;
977
978   if (type==0){
979     Float_t fraction=0;
980     if (findable>0)
981       fraction=(Float_t)found/(Float_t)findable;
982     else
983       fraction=0;
984     return fraction;
985   }
986   return 0;  // undefined type - default value
987 }
988
989
990 //________________________________________________________________________
991 void AliPWG4HighPtTrackQA::Terminate(Option_t *)
992 {
993   // The Terminate() function is the last function to be called during
994   // a query. It always runs on the client, it can be used to present
995   // the results graphically or save the results to file.
996
997 }
998
999 #endif