]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliPWG4HighPtQAMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
- Inherits now from TNamed,
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliPWG4HighPtQAMC.cxx
CommitLineData
fdceab34 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
df943115 24#ifndef ALIPWG4HighPtQAMC_CXX
25#define ALIPWG4HighPtQAMC_CXX
26
fdceab34 27#include "AliPWG4HighPtQAMC.h"
28
29#include "TH1.h"
30#include "TH2.h"
31#include "TH3.h"
b4691ee7 32#include "TProfile.h"
fdceab34 33#include "TList.h"
b4691ee7 34#include "TFile.h"
fdceab34 35#include "TChain.h"
36#include "TH3F.h"
b4691ee7 37#include "TKey.h"
38#include "TSystem.h"
39
40#include "AliAnalysisTask.h"
fdceab34 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"
df943115 49#include "AliLog.h"
b4691ee7 50#include "AliGenPythiaEventHeader.h"
51#include "AliGenCocktailEventHeader.h"
52//#include "AliAnalysisHelperJetTasks.h"
fdceab34 53
54using namespace std; //required for resolving the 'cout' symbol
55
56ClassImp(AliPWG4HighPtQAMC)
57
b4691ee7 58AliPWG4HighPtQAMC::AliPWG4HighPtQAMC()
59: AliAnalysisTask("AliPWG4HighPtQAMC", ""),
fdceab34 60 fESD(0),
b1cd0099 61 fMC(0),
fdceab34 62 fTrackCuts(0),
63 fTrackCutsITS(0),
71e77a79 64 fTrackType(0),
b1cd0099 65 fPtMax(100.),
b4691ee7 66 fAvgTrials(1),
8f0faa80 67 fNEventAll(0),
68 fNEventSel(0),
b4691ee7 69 fh1Xsec(0),
70 fh1Trials(0),
71 fh1PtHard(0),
72 fh1PtHardTrials(0),
df943115 73 fPtAll(0),
74 fPtSel(0),
fdceab34 75 fPtAllminPtMCvsPtAll(0),
76 fPtAllminPtMCvsPtAllNPointTPC(0),
77 fPtAllminPtMCvsPtAllDCAR(0),
78 fPtAllminPtMCvsPtAllDCAZ(0),
79 fPtAllminPtMCvsPtAllPhi(0),
80 fPtAllminPtMCvsPtAllNPointITS(0),
81 fPtAllminPtMCvsPtAllNSigmaToVertex(0),
82 fPtAllminPtMCvsPtAllChi2C(0),
83 fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
84 fPtAllMC(0),
85 fPtSelMC(0),
86 fPtSelMCITS(0),
87 fHistList(0),
88 fPtSelITS(0),
89 fPtITSminPtMCvsPtITS(0),
90 fPtITSminPtMCvsPtITSNPointTPC(0),
91 fPtITSminPtMCvsPtITSDCAR(0),
92 fPtITSminPtMCvsPtITSDCAZ(0),
93 fPtITSminPtMCvsPtITSPhi(0),
94 fPtITSminPtMCvsPtITSNPointITS(0),
95 fPtITSminPtMCvsPtITSNSigmaToVertex(0),
96 fPtITSminPtMCvsPtITSChi2C(0),
97 fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
98 fHistListITS(0)
99{
df943115 100
fdceab34 101}
102//________________________________________________________________________
103AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
b4691ee7 104 AliAnalysisTask(name,""),
fdceab34 105 fESD(0),
b1cd0099 106 fMC(0),
fdceab34 107 fTrackCuts(),
108 fTrackCutsITS(),
71e77a79 109 fTrackType(0),
b1cd0099 110 fPtMax(100.),
b4691ee7 111 fAvgTrials(1),
8f0faa80 112 fNEventAll(0),
113 fNEventSel(0),
b4691ee7 114 fh1Xsec(0),
115 fh1Trials(0),
116 fh1PtHard(0),
117 fh1PtHardTrials(0),
df943115 118 fPtAll(0),
119 fPtSel(0),
fdceab34 120 fPtAllminPtMCvsPtAll(0),
121 fPtAllminPtMCvsPtAllNPointTPC(0),
122 fPtAllminPtMCvsPtAllDCAR(0),
123 fPtAllminPtMCvsPtAllDCAZ(0),
124 fPtAllminPtMCvsPtAllPhi(0),
125 fPtAllminPtMCvsPtAllNPointITS(0),
126 fPtAllminPtMCvsPtAllNSigmaToVertex(0),
127 fPtAllminPtMCvsPtAllChi2C(0),
128 fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
129 fPtAllMC(0),
130 fPtSelMC(0),
131 fPtSelMCITS(0),
132 fHistList(0),
133 fPtSelITS(0),
134 fPtITSminPtMCvsPtITS(0),
135 fPtITSminPtMCvsPtITSNPointTPC(0),
136 fPtITSminPtMCvsPtITSDCAR(0),
137 fPtITSminPtMCvsPtITSDCAZ(0),
138 fPtITSminPtMCvsPtITSPhi(0),
139 fPtITSminPtMCvsPtITSNPointITS(0),
140 fPtITSminPtMCvsPtITSNSigmaToVertex(0),
141 fPtITSminPtMCvsPtITSChi2C(0),
142 fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
143 fHistListITS(0)
144{
145 //
146 // Constructor. Initialization of Inputs and Outputs
147 //
f51451be 148 AliDebug(2,Form("AliPWG4HighPtQAMC Calling Constructor"));
fdceab34 149 // Input slot #0 works with a TChain ESD
150 DefineInput(0, TChain::Class());
b4691ee7 151 // Output slot #0, #1 write into a TList
fdceab34 152 DefineOutput(0, TList::Class());
fdceab34 153 DefineOutput(1, TList::Class());
df943115 154}
fdceab34 155
156//________________________________________________________________________
157void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
158{
b1cd0099 159 // Connect ESD and MC event handler here
fdceab34 160 // Called once
df943115 161 AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
fdceab34 162 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
163 if (!tree) {
b1cd0099 164 AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
165 return;
fdceab34 166 }
b1cd0099 167
168 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
169
170 if (!esdH) {
171 AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
172 return;
173 } else
174 fESD = esdH->GetEvent();
175
176 AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
177 if (!eventHandler) {
178 AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
179 }
180 else
181 fMC = eventHandler->MCEvent();
182
fdceab34 183}
184
b4691ee7 185
fdceab34 186//________________________________________________________________________
187void AliPWG4HighPtQAMC::CreateOutputObjects() {
188 //Create output objects
df943115 189 AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
190
191 Bool_t oldStatus = TH1::AddDirectoryStatus();
192 TH1::AddDirectory(kFALSE);
193
fdceab34 194 OpenFile(0);
195 fHistList = new TList();
b4691ee7 196 fHistList->SetOwner(kTRUE);
fdceab34 197 OpenFile(1);
198 fHistListITS = new TList();
b4691ee7 199 fHistListITS->SetOwner(kTRUE);
fdceab34 200
201 Int_t fgkNPhiBins=18;
202 Float_t kMinPhi = 0.;
203 Float_t kMaxPhi = 2.*TMath::Pi();
204
205 Int_t fgkNPtBins=98;
206 Float_t fgkPtMin=2.;
b1cd0099 207 Float_t fgkPtMax=fPtMax;
8f0faa80 208 Int_t fgkResPtBins=80;
fdceab34 209
8f0faa80 210 fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
211 fHistList->Add(fNEventAll);
212 fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
213 fHistList->Add(fNEventSel);
b4691ee7 214
215 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
216 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
217 fHistList->Add(fh1Xsec);
218
219 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
220 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
221 fHistList->Add(fh1Trials);
222
223 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
224 fHistList->Add(fh1PtHard);
225 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
226 fHistList->Add(fh1PtHardTrials);
227
fdceab34 228 fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
229 fHistList->Add(fPtAll);
230 fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
231 fHistList->Add(fPtSel);
232
233 fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
234 fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{MC}");
235 fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
236 fHistList->Add(fPtAllminPtMCvsPtAll);
237
238 fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
239 fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
240 fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
241 fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{point,TPC}");
242 fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
243
8f0faa80 244 fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
fdceab34 245 fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
246 fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
247 fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
248 fHistList->Add(fPtAllminPtMCvsPtAllDCAR);
249
8f0faa80 250 fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
fdceab34 251 fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
252 fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
253 fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
254 fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);
255
256 fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
257 fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
258 fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
259 fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
260 fHistList->Add(fPtAllminPtMCvsPtAllPhi);
261
262 fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
263 fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
264 fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
265 fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
266 fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
267
268 fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
269 fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
270 fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
271 fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
272 fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);
273
274 fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
275 fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
276 fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
277 fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
278 fHistList->Add(fPtAllminPtMCvsPtAllChi2C);
279
280 fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
281 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
282 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
283 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
284 fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);
285
286 //ITSrefit
287 fPtSelITS = new TH1F("fPtSelITSrefit","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
288 fHistListITS->Add(fPtSelITS);
289
290 fPtITSminPtMCvsPtITS = new TH2F("fPtITSminPtMCvsPtITS","PtITSminPtMCvsPtITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
291 fPtITSminPtMCvsPtITS->SetXTitle("p_{t}^{MC}");
292 fPtITSminPtMCvsPtITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
293 fHistListITS->Add(fPtITSminPtMCvsPtITS);
294
295 fPtITSminPtMCvsPtITSNPointTPC = new TH3F("fPtITSminPtMCvsPtITSNPointTPC","PtITSminPtMCvsPtITSNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
296 fPtITSminPtMCvsPtITSNPointTPC->SetXTitle("p_{t}^{MC}");
297 fPtITSminPtMCvsPtITSNPointTPC->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
298 fPtITSminPtMCvsPtITSNPointTPC->SetZTitle("N_{point,TPC}");
299 fHistListITS->Add(fPtITSminPtMCvsPtITSNPointTPC);
300
8f0faa80 301 fPtITSminPtMCvsPtITSDCAR = new TH3F("fPtITSminPtMCvsPtITSDCAR","PtITSminPtMCvsPtITSDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
fdceab34 302 fPtITSminPtMCvsPtITSDCAR->SetXTitle("p_{t}^{MC}");
303 fPtITSminPtMCvsPtITSDCAR->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
304 fPtITSminPtMCvsPtITSDCAR->SetZTitle("DCA_{R}");
305 fHistListITS->Add(fPtITSminPtMCvsPtITSDCAR);
306
8f0faa80 307 fPtITSminPtMCvsPtITSDCAZ = new TH3F("fPtITSminPtMCvsPtITSDCAZ","PtITSminPtMCvsPtITSDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
fdceab34 308 fPtITSminPtMCvsPtITSDCAZ->SetXTitle("p_{t}^{MC}");
309 fPtITSminPtMCvsPtITSDCAZ->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
310 fPtITSminPtMCvsPtITSDCAZ->SetZTitle("DCA_{Z}");
311 fHistListITS->Add(fPtITSminPtMCvsPtITSDCAZ);
312
313 fPtITSminPtMCvsPtITSPhi = new TH3F("fPtITSminPtMCvsPtITSPhi","PtITSminPtMCvsPtITSPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
314 fPtITSminPtMCvsPtITSPhi->SetXTitle("p_{t}^{MC}");
315 fPtITSminPtMCvsPtITSPhi->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
316 fPtITSminPtMCvsPtITSPhi->SetZTitle("#phi");
317 fHistListITS->Add(fPtITSminPtMCvsPtITSPhi);
318
319 fPtITSminPtMCvsPtITSNPointITS = new TH3F("fPtITSminPtMCvsPtITSNPointITS","PtITSminPtMCvsPtITSNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
320 fPtITSminPtMCvsPtITSNPointITS->SetXTitle("p_{t}^{MC}");
321 fPtITSminPtMCvsPtITSNPointITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})}");
322 fPtITSminPtMCvsPtITSNPointITS->SetZTitle("N_{point,ITS}}");
323 fHistListITS->Add(fPtITSminPtMCvsPtITSNPointITS);
324
325 fPtITSminPtMCvsPtITSNSigmaToVertex = new TH3F("fPtITSminPtMCvsPtITSNSigmaToVertex","PtITSminPtMCvsPtITSNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
326 fPtITSminPtMCvsPtITSNSigmaToVertex->SetXTitle("p_{t}^{MC}");
327 fPtITSminPtMCvsPtITSNSigmaToVertex->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
328 fPtITSminPtMCvsPtITSNSigmaToVertex->SetZTitle("N#sigma to vertex");
329 fHistListITS->Add(fPtITSminPtMCvsPtITSNSigmaToVertex);
330
331 fPtITSminPtMCvsPtITSChi2C = new TH3F("fPtITSminPtMCvsPtITSChi2C","PtITSminPtMCvsPtITSChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
332 fPtITSminPtMCvsPtITSChi2C->SetXTitle("p_{t}^{MC}");
333 fPtITSminPtMCvsPtITSChi2C->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
334 fPtITSminPtMCvsPtITSChi2C->SetZTitle("Constrained #chi^{2}");
335 fHistListITS->Add(fPtITSminPtMCvsPtITSChi2C);
336
337 fPtITSminPtMCvsPtITSRel1PtUncertainty = new TH3F("fPtITSminPtMCvsPtITSRel1PtUncertainty","PtITSminPtMCvsPtITSRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
338 fPtITSminPtMCvsPtITSRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
339 fPtITSminPtMCvsPtITSRel1PtUncertainty->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
340 fPtITSminPtMCvsPtITSRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
341 fHistListITS->Add(fPtITSminPtMCvsPtITSRel1PtUncertainty);
342
343 fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
344 fHistList->Add(fPtAllMC);
345 fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
346 fHistList->Add(fPtSelMC);
347 fPtSelMCITS = new TH1F("fPtSelMCITS","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
348 fHistList->Add(fPtSelMCITS);
349
df943115 350 TH1::AddDirectory(oldStatus);
351
fdceab34 352}
353//________________________________________________________________________
354void AliPWG4HighPtQAMC::Exec(Option_t *) {
355 // Main loop
356 // Called for each event
df943115 357 AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
8f0faa80 358 // All events without selection
359 fNEventAll->Fill(0.);
360
fdceab34 361 if (!fESD) {
b4691ee7 362 AliDebug(2,Form("ERROR: fInputEvent not available\n"));
8f0faa80 363 PostData(0, fHistList);
364 PostData(1, fHistListITS);
fdceab34 365 return;
366 }
367
bfd58011 368 UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
369 if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
8f0faa80 370 AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
371 // Post output data
372 PostData(0, fHistList);
373 PostData(1, fHistListITS);
374 return;
375 }
376
b1cd0099 377 AliStack* stack = 0x0;
378
379 if(fMC) {
380 AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
381 stack = fMC->Stack(); //Particles Stack
382 AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
383 } else {
384 AliDebug(2,Form("ERROR: Could not retrieve MC eventHandler"));
8f0faa80 385 PostData(0, fHistList);
386 PostData(1, fHistListITS);
fdceab34 387 return;
388 }
389
b4691ee7 390 //___ get MC information __________________________________________________________________
391
392 Double_t ptHard = 0.;
393 Double_t nTrials = 1; // trials for MC trigger weight for real data
394
395 if(fMC){
396 AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(fMC);
397 if(!pythiaGenHeader){
398 AliDebug(2,Form("ERROR: Could not retrieve AliGenPythiaEventHeader"));
399 PostData(0, fHistList);
400 PostData(1, fHistListITS);
401 return;
402 } else {
403 nTrials = pythiaGenHeader->Trials();
404 ptHard = pythiaGenHeader->GetPtHard();
405
406 fh1PtHard->Fill(ptHard);
407 fh1PtHardTrials->Fill(ptHard,nTrials);
408
409 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
410 }
411 }
412
fdceab34 413 const AliESDVertex *vtx = fESD->GetPrimaryVertex();
fdceab34 414 // Need vertex cut
75946d5d 415 TString vtxName(vtx->GetName());
416 if(vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
417 // SPD vertex
418 vtx = fESD->GetPrimaryVertexSPD();
419 if(vtx->GetNContributors()<2) {
420 vtx = 0x0;
421 // Post output data
422 PostData(0, fHistList);
423 PostData(1, fHistListITS);
424 return;
425 }
8f0faa80 426 }
75946d5d 427
428 double primVtx[3];
8f0faa80 429 vtx->GetXYZ(primVtx);
8f0faa80 430 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
431 // Post output data
432 PostData(0, fHistList);
433 PostData(1, fHistListITS);
434 return;
435 }
df943115 436
437 AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
fdceab34 438
8f0faa80 439 if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2) {
440 PostData(0, fHistList);
441 PostData(1, fHistListITS);
442 return;
443 }
fdceab34 444
b4691ee7 445 //Need to keep track of selected events
446 fNEventSel->Fill(0.);
447
fdceab34 448 Int_t nTracks = fESD->GetNumberOfTracks();
b1cd0099 449 AliDebug(2,Form("nTracks ESD%d", nTracks));
df943115 450
cba109a0 451 int nMCtracks = stack->GetNtrack();
452
71e77a79 453 Float_t pt = 0.;
454 Float_t ptMC = 0.;
455 Float_t phi = 0.;
456 Float_t dca2D = 0.;
457 Float_t dcaZ = 0.;
458 Int_t nPointITS = 0;
459 Float_t chi2C = 0.;
460 Float_t nSigmaToVertex = 0.;
461 Float_t relUncertainty1Pt = 0.;
462
fdceab34 463 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
464
71e77a79 465 AliESDtrack *track;
466 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
467 if(!esdtrack) continue;
468 AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)esdtrack->GetTPCInnerParam();
469
470 if(fTrackType==1)
471 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
472 else
473 track = esdtrack;
474
475
fdceab34 476 if(!track) continue;
71e77a79 477
fdceab34 478 Int_t label = TMath::Abs(track->GetLabel());
cba109a0 479 if(label>=nMCtracks)continue;
fdceab34 480 TParticle *particle = stack->Particle(label) ;
481 if(!particle) continue;
482
71e77a79 483 ptMC = particle->Pt();
484
485 if(fTrackType==0) { //Global
486 pt = track->Pt();
487 phi = track->Phi();
488 track->GetImpactParameters(dca2D,dcaZ);
489 }
490 else if(fTrackType==1) { //TPConly
491 pt = trackTPC->Pt();
492 phi = trackTPC->Phi();
493 track->GetImpactParametersTPC(dca2D,dcaZ);
494 }
495 else {continue;}
496
497
fdceab34 498 UChar_t itsMap = track->GetITSClusterMap();
fdceab34 499 for (Int_t i=0; i < 6; i++) {
500 if (itsMap & (1 << i))
501 nPointITS ++;
502 }
71e77a79 503 nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
504 chi2C = track->GetConstrainedChi2();
505 relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
506
fdceab34 507 fPtAll->Fill(pt);
508 fPtAllMC->Fill(ptMC);
71e77a79 509
fdceab34 510 if (fTrackCuts->AcceptTrack(track)) {
511
512 fPtSel->Fill(pt);
513
514 fPtSelMC->Fill(ptMC);
515 fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
516 fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
517 fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
518 fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
519 fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
520 fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
521 fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
522 fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
523 fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
524 }//fTrackCuts selection
525
526
71e77a79 527
fdceab34 528 //ITSrefit selection
71e77a79 529 if (fTrackCutsITS->AcceptTrack(track) && fTrackType==0) {
fdceab34 530
531 fPtSelITS->Fill(pt);
532 fPtSelMCITS->Fill(ptMC);
533 fPtITSminPtMCvsPtITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
534 fPtITSminPtMCvsPtITSNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
535 fPtITSminPtMCvsPtITSDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
536 fPtITSminPtMCvsPtITSDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
537 fPtITSminPtMCvsPtITSPhi->Fill(ptMC,(pt-ptMC)/(pt),phi);
538 fPtITSminPtMCvsPtITSNPointITS->Fill(ptMC,(pt-ptMC)/(pt),nPointITS);
539 fPtITSminPtMCvsPtITSNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
540 fPtITSminPtMCvsPtITSChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
541 fPtITSminPtMCvsPtITSRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
542 }//fTrackCutsITS loop
71e77a79 543
fdceab34 544 }//ESD track loop
545
546 // Post output data
547 PostData(0, fHistList);
548 PostData(1, fHistListITS);
549
550}
b4691ee7 551//________________________________________________________________________
552Bool_t AliPWG4HighPtQAMC::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
553 //
554 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
555 // This is to called in Notify and should provide the path to the AOD/ESD file
556 // Copied from AliAnalysisTaskJetSpectrum2
557 //
558
559 TString file(currFile);
560 fXsec = 0;
561 fTrials = 1;
562
563 if(file.Contains("root_archive.zip#")){
564 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
565 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
566 file.Replace(pos+1,20,"");
567 }
568 else {
569 // not an archive take the basename....
570 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
571 }
572 Printf("%s",file.Data());
573
574
575 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
576 if(!fxsec){
577 // next trial fetch the histgram file
578 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
579 if(!fxsec){
580 // not a severe condition but inciate that we have no information
581 return kFALSE;
582 }
583 else{
584 // find the tlist we want to be independtent of the name so use the Tkey
585 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
586 if(!key){
587 fxsec->Close();
588 return kFALSE;
589 }
590 TList *list = dynamic_cast<TList*>(key->ReadObj());
591 if(!list){
592 fxsec->Close();
593 return kFALSE;
594 }
595 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
596 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
597 fxsec->Close();
598 }
599 } // no tree pyxsec.root
600 else {
601 TTree *xtree = (TTree*)fxsec->Get("Xsection");
602 if(!xtree){
603 fxsec->Close();
604 return kFALSE;
605 }
606 UInt_t ntrials = 0;
607 Double_t xsection = 0;
608 xtree->SetBranchAddress("xsection",&xsection);
609 xtree->SetBranchAddress("ntrials",&ntrials);
610 xtree->GetEntry(0);
611 fTrials = ntrials;
612 fXsec = xsection;
613 fxsec->Close();
614 }
615 return kTRUE;
616}
617//________________________________________________________________________
618Bool_t AliPWG4HighPtQAMC::Notify()
619{
620 //
621 // Implemented Notify() to read the cross sections
622 // and number of trials from pyxsec.root
623 // Copied from AliAnalysisTaskJetSpectrum2
624 //
625
626 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
627 Float_t xsection = 0;
628 Float_t ftrials = 1;
629
630 fAvgTrials = 1;
631 if(tree){
632 TFile *curfile = tree->GetCurrentFile();
633 if (!curfile) {
634 Error("Notify","No current file");
635 return kFALSE;
636 }
637 if(!fh1Xsec||!fh1Trials){
638 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
639 return kFALSE;
640 }
641 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
642 fh1Xsec->Fill("<#sigma>",xsection);
643 // construct a poor man average trials
644 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
645 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
646 }
647 return kTRUE;
648}
649
650//________________________________________________________________________
651AliGenPythiaEventHeader* AliPWG4HighPtQAMC::GetPythiaEventHeader(AliMCEvent *mcEvent){
652
653 if(!mcEvent)return 0;
654 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
655 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
656 if(!pythiaGenHeader){
657 // cocktail ??
658 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
659
660 if (!genCocktailHeader) {
661 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
662 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
663 return 0;
664 }
665 TList* headerList = genCocktailHeader->GetHeaders();
666 for (Int_t i=0; i<headerList->GetEntries(); i++) {
667 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
668 if (pythiaGenHeader)
669 break;
670 }
671 if(!pythiaGenHeader){
672 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
673 return 0;
674 }
675 }
676 return pythiaGenHeader;
677
678}
679
fdceab34 680//________________________________________________________________________
681void AliPWG4HighPtQAMC::Terminate(Option_t *)
682{
df943115 683 // The Terminate() function is the last function to be called during
684 // a query. It always runs on the client, it can be used to present
685 // the results graphically or save the results to file.
fdceab34 686
fdceab34 687}
df943115 688
689#endif