adding new QA tasks from Marta
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliPWG4HighPtQAMC.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 compares the global reconstruction with the MC information
18 // Momentum resolution is stored as function of track cuts and pt.
19 // Output: Histograms for different set of cuts
20 //-----------------------------------------------------------------------
21 // Author : Marta Verweij - UU
22 //-----------------------------------------------------------------------
23
24 #include "AliPWG4HighPtQAMC.h"
25
26 #include "TH1.h"
27 #include "TH2.h"
28 #include "TH3.h"
29 #include "TList.h"
30 #include "TChain.h"
31 #include "TH3F.h"
32 #include "AliAnalysisManager.h"
33 #include "AliESDInputHandler.h"
34 #include "AliMCEvent.h"
35 #include "AliMCEventHandler.h"
36 #include "AliStack.h"
37 #include "AliESDtrack.h"
38 #include "AliESDtrackCuts.h"
39 #include "AliExternalTrackParam.h"
40
41
42 using namespace std; //required for resolving the 'cout' symbol
43
44 ClassImp(AliPWG4HighPtQAMC)
45
46 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(): AliAnalysisTask("AliPWG4HighPtQAMC", ""), 
47   fESD(0), 
48   fTrackCuts(0), 
49   fTrackCutsITS(0),
50   fNEvent(0), // just to avoid warnings, inititialized in InitPointers too
51   fPtAll(0),  //
52   fPtSel(0),  //
53   fPtAllminPtMCvsPtAll(0),
54   fPtAllminPtMCvsPtAllNPointTPC(0),
55   fPtAllminPtMCvsPtAllDCAR(0),
56   fPtAllminPtMCvsPtAllDCAZ(0),
57   fPtAllminPtMCvsPtAllPhi(0),
58   fPtAllminPtMCvsPtAllNPointITS(0),
59   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
60   fPtAllminPtMCvsPtAllChi2C(0),
61   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
62   fPtAllMC(0),
63   fPtSelMC(0),
64   fPtSelMCITS(0),
65   fHistList(0),
66   fPtSelITS(0),
67   fPtITSminPtMCvsPtITS(0),
68   fPtITSminPtMCvsPtITSNPointTPC(0),
69   fPtITSminPtMCvsPtITSDCAR(0),
70   fPtITSminPtMCvsPtITSDCAZ(0),
71   fPtITSminPtMCvsPtITSPhi(0),
72   fPtITSminPtMCvsPtITSNPointITS(0),
73   fPtITSminPtMCvsPtITSNSigmaToVertex(0),
74   fPtITSminPtMCvsPtITSChi2C(0),
75   fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
76   fHistListITS(0)
77 {
78   InitHistPointers();
79 }
80 //________________________________________________________________________
81 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name): 
82   AliAnalysisTask(name, ""), 
83   fESD(0),
84   fTrackCuts(),
85   fTrackCutsITS(),
86   fNEvent(0), // just to avoid warnings, inititialized in InitPointers too
87   fPtAll(0),  //
88   fPtSel(0),  //
89   fPtAllminPtMCvsPtAll(0),
90   fPtAllminPtMCvsPtAllNPointTPC(0),
91   fPtAllminPtMCvsPtAllDCAR(0),
92   fPtAllminPtMCvsPtAllDCAZ(0),
93   fPtAllminPtMCvsPtAllPhi(0),
94   fPtAllminPtMCvsPtAllNPointITS(0),
95   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
96   fPtAllminPtMCvsPtAllChi2C(0),
97   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
98   fPtAllMC(0),
99   fPtSelMC(0),
100   fPtSelMCITS(0),
101   fHistList(0),
102   fPtSelITS(0),
103   fPtITSminPtMCvsPtITS(0),
104   fPtITSminPtMCvsPtITSNPointTPC(0),
105   fPtITSminPtMCvsPtITSDCAR(0),
106   fPtITSminPtMCvsPtITSDCAZ(0),
107   fPtITSminPtMCvsPtITSPhi(0),
108   fPtITSminPtMCvsPtITSNPointITS(0),
109   fPtITSminPtMCvsPtITSNSigmaToVertex(0),
110   fPtITSminPtMCvsPtITSChi2C(0),
111   fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
112   fHistListITS(0)
113 {
114   //
115   // Constructor. Initialization of Inputs and Outputs
116   //
117   Info("AliPWG4HighPtQAMC","Calling Constructor");
118   // Input slot #0 works with a TChain ESD
119   DefineInput(0, TChain::Class());
120   // Output slot #0 writes into a TList
121   DefineOutput(0, TList::Class());
122   // Output slot #1 writes into a TList
123   DefineOutput(1, TList::Class());
124   // Output slot #2 writes into a TList
125   DefineOutput(2, TList::Class());
126   InitHistPointers();
127   //  TH1::AddDirectory(kFALSE);
128   //  TH2::AddDirectory(kFALSE);
129   // TH3::AddDirectory(kFALSE);
130   }
131
132 //________________________________________________________________________
133 void AliPWG4HighPtQAMC::ConnectInputData(Option_t *) 
134 {
135   // Connect ESD here
136   // Called once
137   printf(">> AliPWG4HighPtQATPConly::ConnectInputData \n");
138   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
139   if (!tree) {
140     Printf("ERROR: Could not read chain from input slot 0");
141   } else {
142     
143     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
144     
145     if (!esdH) {
146       Printf("ERROR: Could not get ESDInputHandler");
147     } else
148       fESD = esdH->GetEvent();
149   }
150 }
151
152 //________________________________________________________________________
153 void AliPWG4HighPtQAMC::InitHistPointers() {
154   //Initialize histograms
155   fNEvent = 0;
156   fPtAll = 0;
157   fPtSel = 0;
158   //Global tracking compared to MC pt
159   fPtAllMC = 0;
160   fPtSelMC = 0;
161   fPtAllminPtMCvsPtAll = 0;
162   fPtAllminPtMCvsPtAllNPointTPC = 0;
163   fPtAllminPtMCvsPtAllDCAR = 0;
164   fPtAllminPtMCvsPtAllDCAZ = 0;
165   fPtAllminPtMCvsPtAllPhi = 0;
166   fPtAllminPtMCvsPtAllNPointITS = 0;
167   fPtAllminPtMCvsPtAllNSigmaToVertex = 0;
168   fPtAllminPtMCvsPtAllChi2C = 0;
169   fPtAllminPtMCvsPtAllRel1PtUncertainty = 0;
170   //ITSrefit histos compared to MC pt
171   fPtSelITS = 0;
172   fPtITSminPtMCvsPtITS = 0;
173   fPtITSminPtMCvsPtITSNPointTPC = 0;
174   fPtITSminPtMCvsPtITSDCAR = 0;
175   fPtITSminPtMCvsPtITSDCAZ = 0;
176   fPtITSminPtMCvsPtITSPhi = 0;
177   fPtITSminPtMCvsPtITSNPointITS = 0;
178   fPtITSminPtMCvsPtITSNSigmaToVertex = 0;
179   fPtITSminPtMCvsPtITSChi2C = 0;
180   fPtITSminPtMCvsPtITSRel1PtUncertainty = 0;
181 }
182 //________________________________________________________________________
183 void AliPWG4HighPtQAMC::CreateOutputObjects() {
184   //Create output objects
185   printf(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n");
186   OpenFile(0);
187   fHistList = new TList();
188   OpenFile(1);
189   fHistListITS = new TList();
190
191   Int_t fgkNPhiBins=18;
192   Float_t kMinPhi = 0.;
193   Float_t kMaxPhi = 2.*TMath::Pi();
194   
195   Int_t fgkNPtBins=98;
196   Float_t fgkPtMin=2.;
197   Float_t fgkPtMax=100.;
198   Int_t fgkResPtBins=40;
199
200   fNEvent = new TH1F("fNEvent","NEvent",1,-0.5,0.5);
201   fHistList->Add(fNEvent);
202   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
203   fHistList->Add(fPtAll);
204   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
205   fHistList->Add(fPtSel);
206   
207   fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
208   fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{MC}");
209   fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
210   fHistList->Add(fPtAllminPtMCvsPtAll);
211   
212   fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
213   fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
214   fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
215   fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{point,TPC}");
216   fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
217
218   fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,-1.,1.);
219   fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
220   fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
221   fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
222   fHistList->Add(fPtAllminPtMCvsPtAllDCAR);
223
224   fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,-2.,2.);
225   fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
226   fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
227   fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
228   fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);
229
230   fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
231   fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
232   fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
233   fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
234   fHistList->Add(fPtAllminPtMCvsPtAllPhi);
235
236   fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
237   fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
238   fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
239   fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
240   fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
241   
242   fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
243   fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
244   fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
245   fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
246   fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);
247
248   fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
249   fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
250   fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
251   fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
252   fHistList->Add(fPtAllminPtMCvsPtAllChi2C);
253
254   fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
255   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
256   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
257   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
258   fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);
259
260   //ITSrefit
261   fPtSelITS = new TH1F("fPtSelITSrefit","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
262   fHistListITS->Add(fPtSelITS);
263   
264   fPtITSminPtMCvsPtITS = new TH2F("fPtITSminPtMCvsPtITS","PtITSminPtMCvsPtITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
265   fPtITSminPtMCvsPtITS->SetXTitle("p_{t}^{MC}");
266   fPtITSminPtMCvsPtITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
267   fHistListITS->Add(fPtITSminPtMCvsPtITS);
268   
269   fPtITSminPtMCvsPtITSNPointTPC = new TH3F("fPtITSminPtMCvsPtITSNPointTPC","PtITSminPtMCvsPtITSNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
270   fPtITSminPtMCvsPtITSNPointTPC->SetXTitle("p_{t}^{MC}");
271   fPtITSminPtMCvsPtITSNPointTPC->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
272   fPtITSminPtMCvsPtITSNPointTPC->SetZTitle("N_{point,TPC}");
273   fHistListITS->Add(fPtITSminPtMCvsPtITSNPointTPC);
274     
275   fPtITSminPtMCvsPtITSDCAR = new TH3F("fPtITSminPtMCvsPtITSDCAR","PtITSminPtMCvsPtITSDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,-1.,1.);
276   fPtITSminPtMCvsPtITSDCAR->SetXTitle("p_{t}^{MC}");
277   fPtITSminPtMCvsPtITSDCAR->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
278   fPtITSminPtMCvsPtITSDCAR->SetZTitle("DCA_{R}");
279   fHistListITS->Add(fPtITSminPtMCvsPtITSDCAR);
280   
281   fPtITSminPtMCvsPtITSDCAZ = new TH3F("fPtITSminPtMCvsPtITSDCAZ","PtITSminPtMCvsPtITSDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,-2.,2.);
282   fPtITSminPtMCvsPtITSDCAZ->SetXTitle("p_{t}^{MC}");
283   fPtITSminPtMCvsPtITSDCAZ->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
284   fPtITSminPtMCvsPtITSDCAZ->SetZTitle("DCA_{Z}");
285   fHistListITS->Add(fPtITSminPtMCvsPtITSDCAZ);
286   
287   fPtITSminPtMCvsPtITSPhi = new TH3F("fPtITSminPtMCvsPtITSPhi","PtITSminPtMCvsPtITSPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
288   fPtITSminPtMCvsPtITSPhi->SetXTitle("p_{t}^{MC}");
289   fPtITSminPtMCvsPtITSPhi->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
290   fPtITSminPtMCvsPtITSPhi->SetZTitle("#phi");
291   fHistListITS->Add(fPtITSminPtMCvsPtITSPhi);
292   
293   fPtITSminPtMCvsPtITSNPointITS = new TH3F("fPtITSminPtMCvsPtITSNPointITS","PtITSminPtMCvsPtITSNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
294   fPtITSminPtMCvsPtITSNPointITS->SetXTitle("p_{t}^{MC}");
295   fPtITSminPtMCvsPtITSNPointITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})}");
296   fPtITSminPtMCvsPtITSNPointITS->SetZTitle("N_{point,ITS}}");
297   fHistListITS->Add(fPtITSminPtMCvsPtITSNPointITS); 
298   
299   fPtITSminPtMCvsPtITSNSigmaToVertex = new TH3F("fPtITSminPtMCvsPtITSNSigmaToVertex","PtITSminPtMCvsPtITSNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
300   fPtITSminPtMCvsPtITSNSigmaToVertex->SetXTitle("p_{t}^{MC}");
301   fPtITSminPtMCvsPtITSNSigmaToVertex->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
302   fPtITSminPtMCvsPtITSNSigmaToVertex->SetZTitle("N#sigma to vertex");
303   fHistListITS->Add(fPtITSminPtMCvsPtITSNSigmaToVertex);
304
305   fPtITSminPtMCvsPtITSChi2C = new TH3F("fPtITSminPtMCvsPtITSChi2C","PtITSminPtMCvsPtITSChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
306   fPtITSminPtMCvsPtITSChi2C->SetXTitle("p_{t}^{MC}");
307   fPtITSminPtMCvsPtITSChi2C->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
308   fPtITSminPtMCvsPtITSChi2C->SetZTitle("Constrained #chi^{2}");
309   fHistListITS->Add(fPtITSminPtMCvsPtITSChi2C);
310
311   fPtITSminPtMCvsPtITSRel1PtUncertainty = new TH3F("fPtITSminPtMCvsPtITSRel1PtUncertainty","PtITSminPtMCvsPtITSRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
312   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
313   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
314   fPtITSminPtMCvsPtITSRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
315   fHistListITS->Add(fPtITSminPtMCvsPtITSRel1PtUncertainty);
316
317   fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
318   fHistList->Add(fPtAllMC);
319   fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
320   fHistList->Add(fPtSelMC);
321   fPtSelMCITS = new TH1F("fPtSelMCITS","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
322   fHistList->Add(fPtSelMCITS);
323   
324 }
325 //________________________________________________________________________
326 void AliPWG4HighPtQAMC::Exec(Option_t *) {  
327   // Main loop
328   // Called for each event
329   printf(">> AliPWG4HighPtQATPConly::Exec \n");  
330   
331   if (!fESD) {
332     Printf("ERROR: fESD not available");
333     return;
334   }
335
336  AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
337   if (!eventHandler) {
338     Printf("ERROR: Could not retrieve MC event handler");
339     return;
340   }
341
342   AliMCEvent* mcEvent = eventHandler->MCEvent();
343   if (!mcEvent) {
344     Printf("ERROR: Could not retrieve MC event");
345     return;
346   }
347
348   Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
349
350   if (!fESD) {
351     Printf("ERROR: fESD not available");
352     return;
353   }
354
355   AliStack* stack = mcEvent->Stack();                //Particles Stack
356
357   Printf("MC particles stack: %d", stack->GetNtrack());
358
359   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
360
361   // Need vertex cut
362   if (vtx->GetNContributors() < 2)
363     return;
364
365   printf("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors());
366   // Need to keep track of evts without vertex
367   fNEvent->Fill(0.);
368
369   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2) return;
370   Int_t nTracks = fESD->GetNumberOfTracks();
371   printf("nTracks %d\n", nTracks);
372   
373   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
374     
375     AliESDtrack *track = fESD->GetTrack(iTrack);
376     if(!track) continue;
377     Int_t label = TMath::Abs(track->GetLabel());
378     TParticle *particle = stack->Particle(label) ;
379     if(!particle) continue;
380
381     Float_t pt = track->Pt();
382     Float_t ptMC = particle->Pt();
383     Float_t phi = track->Phi();
384     Float_t dca2D, dcaZ;
385     track->GetImpactParameters(dca2D,dcaZ);
386     UChar_t itsMap = track->GetITSClusterMap();
387     Int_t nPointITS = 0;
388     for (Int_t i=0; i < 6; i++) {
389       if (itsMap & (1 << i))
390         nPointITS ++;
391     }
392     Float_t nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
393     Float_t chi2C = track->GetConstrainedChi2();
394     Float_t relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
395
396     fPtAll->Fill(pt);
397     fPtAllMC->Fill(ptMC);
398     
399     if (fTrackCuts->AcceptTrack(track)) {
400
401       fPtSel->Fill(pt);
402       
403       fPtSelMC->Fill(ptMC);
404       fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
405       fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
406       fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
407       fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
408       fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
409       fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
410       fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
411       fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
412       fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
413     }//fTrackCuts selection
414     
415     
416     //ITSrefit selection
417     if (fTrackCutsITS->AcceptTrack(track)) {
418       
419       fPtSelITS->Fill(pt);
420       fPtSelMCITS->Fill(ptMC);
421       fPtITSminPtMCvsPtITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
422       fPtITSminPtMCvsPtITSNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
423       fPtITSminPtMCvsPtITSDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
424       fPtITSminPtMCvsPtITSDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
425       fPtITSminPtMCvsPtITSPhi->Fill(ptMC,(pt-ptMC)/(pt),phi);
426       fPtITSminPtMCvsPtITSNPointITS->Fill(ptMC,(pt-ptMC)/(pt),nPointITS);
427       fPtITSminPtMCvsPtITSNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
428       fPtITSminPtMCvsPtITSChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
429       fPtITSminPtMCvsPtITSRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
430     }//fTrackCutsITS loop
431       
432   }//ESD track loop
433    
434   // Post output data
435   PostData(0, fHistList);
436   PostData(1, fHistListITS);
437
438 }
439 //________________________________________________________________________
440 void AliPWG4HighPtQAMC::Terminate(Option_t *)
441 {
442   printf("->AliPWG4HighPtQATPConly::Terminate \n");
443
444   printf("-<AliPWG4HighPtQATPConly::Terminate \n");
445 }