1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 //-----------------------------------------------------------------------
24 #ifndef ALIPWG4HighPtQAMC_CXX
25 #define ALIPWG4HighPtQAMC_CXX
27 #include "AliPWG4HighPtQAMC.h"
40 #include "AliAnalysisTask.h"
41 #include "AliAnalysisManager.h"
42 #include "AliESDInputHandler.h"
43 #include "AliMCEvent.h"
44 #include "AliMCEventHandler.h"
46 #include "AliESDtrack.h"
47 #include "AliESDtrackCuts.h"
48 #include "AliExternalTrackParam.h"
50 #include "AliGenPythiaEventHeader.h"
51 #include "AliGenCocktailEventHeader.h"
52 //#include "AliAnalysisHelperJetTasks.h"
54 using namespace std; //required for resolving the 'cout' symbol
56 ClassImp(AliPWG4HighPtQAMC)
58 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC()
59 : AliAnalysisTask("AliPWG4HighPtQAMC", ""),
67 fSigmaConstrainedMax(5.),
86 fPtAllminPtMCvsPtAll(0),
87 fPtAllvsPtMCvsMult(0),
88 fPtAllminPtMCvsPtAllvsMult(0),
89 fPtAllminPtMCvsPtAllNPointTPC(0),
90 fPtAllminPtMCvsPtAllDCAR(0),
91 fPtAllminPtMCvsPtAllDCAZ(0),
92 fPtAllminPtMCvsPtAllPhi(0),
93 fPtAllminPtMCvsPtAllNPointITS(0),
94 fPtAllminPtMCvsPtAllNSigmaToVertex(0),
95 fPtAllminPtMCvsPtAllChi2C(0),
96 fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
102 fPtITSminPtMCvsPtITS(0),
103 fPtITSminPtMCvsPtITSNPointTPC(0),
104 fPtITSminPtMCvsPtITSDCAR(0),
105 fPtITSminPtMCvsPtITSDCAZ(0),
106 fPtITSminPtMCvsPtITSPhi(0),
107 fPtITSminPtMCvsPtITSNPointITS(0),
108 fPtITSminPtMCvsPtITSNSigmaToVertex(0),
109 fPtITSminPtMCvsPtITSChi2C(0),
110 fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
115 //________________________________________________________________________
116 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
117 AliAnalysisTask(name,""),
125 fSigmaConstrainedMax(5.),
141 fNPointTPCMultRec(0),
144 fPtAllminPtMCvsPtAll(0),
145 fPtAllvsPtMCvsMult(0),
146 fPtAllminPtMCvsPtAllvsMult(0),
147 fPtAllminPtMCvsPtAllNPointTPC(0),
148 fPtAllminPtMCvsPtAllDCAR(0),
149 fPtAllminPtMCvsPtAllDCAZ(0),
150 fPtAllminPtMCvsPtAllPhi(0),
151 fPtAllminPtMCvsPtAllNPointITS(0),
152 fPtAllminPtMCvsPtAllNSigmaToVertex(0),
153 fPtAllminPtMCvsPtAllChi2C(0),
154 fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
160 fPtITSminPtMCvsPtITS(0),
161 fPtITSminPtMCvsPtITSNPointTPC(0),
162 fPtITSminPtMCvsPtITSDCAR(0),
163 fPtITSminPtMCvsPtITSDCAZ(0),
164 fPtITSminPtMCvsPtITSPhi(0),
165 fPtITSminPtMCvsPtITSNPointITS(0),
166 fPtITSminPtMCvsPtITSNSigmaToVertex(0),
167 fPtITSminPtMCvsPtITSChi2C(0),
168 fPtITSminPtMCvsPtITSRel1PtUncertainty(0),
172 // Constructor. Initialization of Inputs and Outputs
174 AliDebug(2,Form("AliPWG4HighPtQAMC Calling Constructor"));
175 // Input slot #0 works with a TChain ESD
176 DefineInput(0, TChain::Class());
177 // Output slot #0, #1 write into a TList
178 DefineOutput(0, TList::Class());
179 DefineOutput(1, TList::Class());
182 //________________________________________________________________________
183 void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
185 // Connect ESD and MC event handler here
187 AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
188 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
190 AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
194 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
197 AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
200 fESD = esdH->GetEvent();
202 AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
204 AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
207 fMC = eventHandler->MCEvent();
212 //________________________________________________________________________
213 void AliPWG4HighPtQAMC::CreateOutputObjects() {
214 //Create output objects
215 AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
217 Bool_t oldStatus = TH1::AddDirectoryStatus();
218 TH1::AddDirectory(kFALSE);
221 fHistList = new TList();
222 fHistList->SetOwner(kTRUE);
224 fHistListITS = new TList();
225 fHistListITS->SetOwner(kTRUE);
227 Int_t fgkNPhiBins=18;
228 Float_t kMinPhi = 0.;
229 Float_t kMaxPhi = 2.*TMath::Pi();
231 Int_t fgkNPtBins=100;
232 Float_t fgkPtMin=0.;//2.;
233 Float_t fgkPtMax=fPtMax;
234 Int_t fgkResPtBins=80;
236 Int_t fgkNMultBins = 50;
237 Float_t fgkMultMin = 0;
238 Float_t fgkMultMax = 500;
240 fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
241 fHistList->Add(fNEventAll);
242 fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
243 fHistList->Add(fNEventSel);
244 fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
246 fNEventReject->Fill("noESD",0);
247 fNEventReject->Fill("Trigger",0);
248 fNEventReject->Fill("noMCEvent",0);
249 fNEventReject->Fill("noStack",0);
250 fNEventReject->Fill("NTracks<2",0);
251 fNEventReject->Fill("noVTX",0);
252 fNEventReject->Fill("VtxStatus",0);
253 fNEventReject->Fill("NCont<2",0);
254 fNEventReject->Fill("ZVTX>10",0);
255 fHistList->Add(fNEventReject);
257 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
258 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
259 fHistList->Add(fh1Xsec);
261 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
262 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
263 fHistList->Add(fh1Trials);
265 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
266 fHistList->Add(fh1PtHard);
267 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
268 fHistList->Add(fh1PtHardTrials);
270 fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
271 fHistList->Add(fPtAll);
272 fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
273 fHistList->Add(fPtSel);
275 fPtSelFakes = new TH1F("fPtSelFakes","PtSelFakes",fgkNPtBins, fgkPtMin, fgkPtMax);
276 fHistList->Add(fPtSelFakes);
278 fNPointTPCFakes = new TH1F("fNPointTPCFakes","fNPointTPCFakes",160,0.5,160.5);
279 fHistList->Add(fNPointTPCFakes);
281 fPtSelLargeLabel = new TH1F("fPtSelLargeLabel","PtSelLargeLabel",fgkNPtBins, fgkPtMin, fgkPtMax);
282 fHistList->Add(fPtSelLargeLabel);
284 fMultRec = new TH1F("fMultRec","Multiple reconstruction of tracks",20,0,20);
285 fHistList->Add(fMultRec);
287 fNPointTPCMultRec = new TH1F("fNPointTPCMultRec","Multiple reconstruction of tracks",160,0.5,160.5);
288 fHistList->Add(fNPointTPCMultRec);
290 fDeltaPtMultRec = new TH2F("fDeltaPtMultRec","Delta pT vs pT for multiple reconstructed tracks",100,0.,50.,100,-20.,20.);
291 fHistList->Add(fDeltaPtMultRec);
293 fPtAllvsPtMC = new TH2F("fPtAllvsPtMC","fPtAllvsPtMC;p_{T,MC};p_{T,rec}",fgkNPtBins, fgkPtMin,fgkPtMax,fgkNPtBins, fgkPtMin,fgkPtMax);
294 fHistList->Add(fPtAllvsPtMC);
296 fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
297 fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{MC}");
298 fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
299 fHistList->Add(fPtAllminPtMCvsPtAll);
301 fPtAllvsPtMCvsMult = new TH3F("fPtAllvsPtMCvsMult","fPtAllvsPtMCvsMult;p_{T,MC};p_{T,rec};N_{tracks}",fgkNPtBins, fgkPtMin,fgkPtMax,fgkNPtBins, fgkPtMin,fgkPtMax,fgkNMultBins,fgkMultMin,fgkMultMax);
302 fHistList->Add(fPtAllvsPtMCvsMult);
304 fPtAllminPtMCvsPtAllvsMult = new TH3F("fPtAllminPtMCvsPtAllvsMult","fPtAllminPtMCvsPtAllvsMult",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNMultBins,fgkMultMin,fgkMultMax);
305 fPtAllminPtMCvsPtAllvsMult->SetXTitle("p_{t}^{MC}");
306 fPtAllminPtMCvsPtAllvsMult->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
307 fPtAllminPtMCvsPtAllvsMult->SetZTitle("N_{tracks}");
308 fHistList->Add(fPtAllminPtMCvsPtAllvsMult);
310 fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
311 fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
312 fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
313 fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{point,TPC}");
314 fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
316 fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
317 fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
318 fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
319 fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
320 fHistList->Add(fPtAllminPtMCvsPtAllDCAR);
322 fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
323 fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
324 fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
325 fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
326 fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);
328 fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
329 fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
330 fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
331 fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
332 fHistList->Add(fPtAllminPtMCvsPtAllPhi);
334 fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
335 fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
336 fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
337 fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
338 fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
340 fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
341 fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
342 fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
343 fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
344 fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);
346 fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
347 fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
348 fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
349 fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
350 fHistList->Add(fPtAllminPtMCvsPtAllChi2C);
352 fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
353 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
354 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
355 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
356 fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);
359 fPtSelITS = new TH1F("fPtSelITSrefit","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
360 fHistListITS->Add(fPtSelITS);
362 fPtITSminPtMCvsPtITS = new TH2F("fPtITSminPtMCvsPtITS","PtITSminPtMCvsPtITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
363 fPtITSminPtMCvsPtITS->SetXTitle("p_{t}^{MC}");
364 fPtITSminPtMCvsPtITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
365 fHistListITS->Add(fPtITSminPtMCvsPtITS);
367 fPtITSminPtMCvsPtITSNPointTPC = new TH3F("fPtITSminPtMCvsPtITSNPointTPC","PtITSminPtMCvsPtITSNPointTPC",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,160,0.5,160.5);
368 fPtITSminPtMCvsPtITSNPointTPC->SetXTitle("p_{t}^{MC}");
369 fPtITSminPtMCvsPtITSNPointTPC->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
370 fPtITSminPtMCvsPtITSNPointTPC->SetZTitle("N_{point,TPC}");
371 fHistListITS->Add(fPtITSminPtMCvsPtITSNPointTPC);
373 fPtITSminPtMCvsPtITSDCAR = new TH3F("fPtITSminPtMCvsPtITSDCAR","PtITSminPtMCvsPtITSDCAR",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-0.2,0.2);
374 fPtITSminPtMCvsPtITSDCAR->SetXTitle("p_{t}^{MC}");
375 fPtITSminPtMCvsPtITSDCAR->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
376 fPtITSminPtMCvsPtITSDCAR->SetZTitle("DCA_{R}");
377 fHistListITS->Add(fPtITSminPtMCvsPtITSDCAR);
379 fPtITSminPtMCvsPtITSDCAZ = new TH3F("fPtITSminPtMCvsPtITSDCAZ","PtITSminPtMCvsPtITSDCAZ",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,80,-2.,2.);
380 fPtITSminPtMCvsPtITSDCAZ->SetXTitle("p_{t}^{MC}");
381 fPtITSminPtMCvsPtITSDCAZ->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
382 fPtITSminPtMCvsPtITSDCAZ->SetZTitle("DCA_{Z}");
383 fHistListITS->Add(fPtITSminPtMCvsPtITSDCAZ);
385 fPtITSminPtMCvsPtITSPhi = new TH3F("fPtITSminPtMCvsPtITSPhi","PtITSminPtMCvsPtITSPhi",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,fgkNPhiBins,kMinPhi,kMaxPhi);
386 fPtITSminPtMCvsPtITSPhi->SetXTitle("p_{t}^{MC}");
387 fPtITSminPtMCvsPtITSPhi->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
388 fPtITSminPtMCvsPtITSPhi->SetZTitle("#phi");
389 fHistListITS->Add(fPtITSminPtMCvsPtITSPhi);
391 fPtITSminPtMCvsPtITSNPointITS = new TH3F("fPtITSminPtMCvsPtITSNPointITS","PtITSminPtMCvsPtITSNPointITS",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,9,-0.5,8.5);
392 fPtITSminPtMCvsPtITSNPointITS->SetXTitle("p_{t}^{MC}");
393 fPtITSminPtMCvsPtITSNPointITS->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})}");
394 fPtITSminPtMCvsPtITSNPointITS->SetZTitle("N_{point,ITS}}");
395 fHistListITS->Add(fPtITSminPtMCvsPtITSNPointITS);
397 fPtITSminPtMCvsPtITSNSigmaToVertex = new TH3F("fPtITSminPtMCvsPtITSNSigmaToVertex","PtITSminPtMCvsPtITSNSigmaToVertex",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,40,0.,8.);
398 fPtITSminPtMCvsPtITSNSigmaToVertex->SetXTitle("p_{t}^{MC}");
399 fPtITSminPtMCvsPtITSNSigmaToVertex->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
400 fPtITSminPtMCvsPtITSNSigmaToVertex->SetZTitle("N#sigma to vertex");
401 fHistListITS->Add(fPtITSminPtMCvsPtITSNSigmaToVertex);
403 fPtITSminPtMCvsPtITSChi2C = new TH3F("fPtITSminPtMCvsPtITSChi2C","PtITSminPtMCvsPtITSChi2C",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,20,0.,10.);
404 fPtITSminPtMCvsPtITSChi2C->SetXTitle("p_{t}^{MC}");
405 fPtITSminPtMCvsPtITSChi2C->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
406 fPtITSminPtMCvsPtITSChi2C->SetZTitle("Constrained #chi^{2}");
407 fHistListITS->Add(fPtITSminPtMCvsPtITSChi2C);
409 fPtITSminPtMCvsPtITSRel1PtUncertainty = new TH3F("fPtITSminPtMCvsPtITSRel1PtUncertainty","PtITSminPtMCvsPtITSRel1PtUncertainty",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.,30,0.,0.3);
410 fPtITSminPtMCvsPtITSRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
411 fPtITSminPtMCvsPtITSRel1PtUncertainty->SetYTitle("(1/p_{t}^{ITS}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
412 fPtITSminPtMCvsPtITSRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
413 fHistListITS->Add(fPtITSminPtMCvsPtITSRel1PtUncertainty);
415 fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
416 fHistList->Add(fPtAllMC);
417 fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
418 fHistList->Add(fPtSelMC);
419 fPtSelMCITS = new TH1F("fPtSelMCITS","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
420 fHistList->Add(fPtSelMCITS);
422 TH1::AddDirectory(oldStatus);
424 PostData(0, fHistList);
425 PostData(1, fHistListITS);
429 //________________________________________________________________________
430 Bool_t AliPWG4HighPtQAMC::SelectEvent() {
432 // Decide if event should be selected for analysis
435 // Checks following requirements:
437 // - trigger info from AliPhysicsSelection
438 // - MCevent available
439 // - number of reconstructed tracks > 1
440 // - primary vertex reconstructed
441 // - z-vertex < 10 cm
443 Bool_t selectEvent = kTRUE;
445 //fESD object available?
447 AliDebug(2,Form("ERROR: fInputEvent not available\n"));
448 fNEventReject->Fill("noESD",1);
449 selectEvent = kFALSE;
454 UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
455 if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
456 AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
457 fNEventReject->Fill("Trigger",1);
458 selectEvent = kFALSE;
465 AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
466 fStack = fMC->Stack(); //Particles Stack
468 AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
471 AliDebug(2,Form("ERROR: Could not retrieve MC eventHandler"));
472 fNEventReject->Fill("noStack",1);
473 selectEvent = kFALSE;
477 AliDebug(2,Form("ERROR: Could not retrieve stack"));
478 fNEventReject->Fill("noMCEvent",1);
479 selectEvent = kFALSE;
483 //Check if number of reconstructed tracks is larger than 1
484 if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2) {
485 fNEventReject->Fill("NTracks<2",1);
486 selectEvent = kFALSE;
490 //Check if vertex is reconstructed
491 if(fTrackType==1) fVtx = fESD->GetPrimaryVertexTPC();
492 else fVtx = fESD->GetPrimaryVertexSPD();
495 fNEventReject->Fill("noVTX",1);
496 selectEvent = kFALSE;
500 if(!fVtx->GetStatus()) {
501 fNEventReject->Fill("VtxStatus",1);
502 selectEvent = kFALSE;
507 // TString vtxName(fVtx->GetName());
508 if(fVtx->GetNContributors()<2) {
509 fNEventReject->Fill("NCont<2",1);
510 selectEvent = kFALSE;
514 //Check if z-vertex < 10 cm
516 fVtx->GetXYZ(primVtx);
517 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
518 fNEventReject->Fill("ZVTX>10",1);
519 selectEvent = kFALSE;
523 AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",fVtx->GetTitle(), fVtx->GetStatus(), fVtx->GetNContributors()));
529 //________________________________________________________________________
530 void AliPWG4HighPtQAMC::Exec(Option_t *) {
532 // Called for each event
533 AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
534 // All events without selection
535 fNEventAll->Fill(0.);
539 PostData(0, fHistList);
540 PostData(1, fHistListITS);
544 // ---- Get MC Header information (for MC productions in pThard bins) ----
545 Double_t ptHard = 0.;
546 Double_t nTrials = 1; // trials for MC trigger weight for real data
549 AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(fMC);
551 nTrials = pythiaGenHeader->Trials();
552 ptHard = pythiaGenHeader->GetPtHard();
554 fh1PtHard->Fill(ptHard);
555 fh1PtHardTrials->Fill(ptHard,nTrials);
557 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
561 //Need to keep track of selected events
562 fNEventSel->Fill(0.);
564 Int_t nTracks = fESD->GetNumberOfTracks();
565 AliDebug(2,Form("nTracks ESD%d", nTracks));
567 int nMCtracks = fStack->GetNtrack();
576 Float_t nSigmaToVertex = 0.;
577 Float_t relUncertainty1Pt = 0.;
579 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
582 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
583 if(!esdtrack) continue;
585 if(fTrackType==1 || fTrackType==3)
586 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
587 else if(fTrackType==2) {
588 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
593 AliExternalTrackParam exParam;
594 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
599 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
605 if(fTrackType==1 || fTrackType==2) delete track;
610 //Cut on chi2 of constrained fit
611 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
617 Int_t label = TMath::Abs(track->GetLabel());
618 if(label>=nMCtracks) {
619 if (fTrackCuts->AcceptTrack(track)) { fPtSelLargeLabel->Fill(pt); }
620 if(fTrackType==1 || fTrackType==2) delete track;
623 TParticle *particle = fStack->Particle(label) ;
625 if(fTrackType==1 || fTrackType==2) delete track;
629 ptMC = particle->Pt();
631 if(fTrackType==0) { //Global
634 track->GetImpactParameters(dca2D,dcaZ);
636 else if(fTrackType==1 || fTrackType==2) { //TPConly
639 track->GetImpactParametersTPC(dca2D,dcaZ);
644 UChar_t itsMap = track->GetITSClusterMap();
645 for (Int_t i=0; i < 6; i++) {
646 if (itsMap & (1 << i))
649 nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
650 chi2C = track->GetConstrainedChi2();
651 relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
654 fPtAllMC->Fill(ptMC);
656 if (fTrackCuts->AcceptTrack(track)) {
659 if(track->GetLabel()<0) {
660 fPtSelFakes->Fill(pt);
661 fNPointTPCFakes->Fill(track->GetTPCNcls());
663 fPtSelMC->Fill(ptMC);
664 fPtAllvsPtMC->Fill(ptMC,pt);
665 fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
666 fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
667 fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
668 fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
669 fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
670 fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
671 fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
672 fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
673 fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
675 int mult = fTrackCuts->CountAcceptedTracks(fESD);
676 fPtAllvsPtMCvsMult->Fill(ptMC,pt,mult);
677 fPtAllminPtMCvsPtAllvsMult->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC), mult);
679 //Check if track is reconstructed multiple times
681 for (Int_t iTrack2 = iTrack+1; iTrack2 < nTracks; iTrack2++) {
682 // AliESDtrack *track2 = GetTrackForAnalysis(iTrack2);
684 AliESDtrack *esdtrack2 = fESD->GetTrack(iTrack2);
685 if(!esdtrack2) continue;
688 track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
689 else if(fTrackType==2) {
690 track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
695 AliExternalTrackParam exParam2;
696 Bool_t relate = track2->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam2);
701 track2->Set(exParam2.GetX(),exParam2.GetAlpha(),exParam2.GetParameter(),exParam2.GetCovariance());
706 if(fTrackType==1 || fTrackType==2) delete track2;
710 if (fTrackCuts->AcceptTrack(track2)) {
711 Int_t label2 = TMath::Abs(track2->GetLabel());
713 fNPointTPCMultRec->Fill(track->GetTPCNcls());
714 fNPointTPCMultRec->Fill(track2->GetTPCNcls());
715 fDeltaPtMultRec->Fill(track->Pt(),track->Pt()-track2->Pt());
719 if(fTrackType==1 || fTrackType==2) delete track2;
721 if(multCounter>1) fMultRec->Fill(multCounter);
723 }//fTrackCuts selection
726 if (fTrackCutsITS->AcceptTrack(track) && fTrackType==0) {
729 fPtSelMCITS->Fill(ptMC);
730 fPtITSminPtMCvsPtITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
731 fPtITSminPtMCvsPtITSNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
732 fPtITSminPtMCvsPtITSDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
733 fPtITSminPtMCvsPtITSDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
734 fPtITSminPtMCvsPtITSPhi->Fill(ptMC,(pt-ptMC)/(pt),phi);
735 fPtITSminPtMCvsPtITSNPointITS->Fill(ptMC,(pt-ptMC)/(pt),nPointITS);
736 fPtITSminPtMCvsPtITSNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
737 fPtITSminPtMCvsPtITSChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
738 fPtITSminPtMCvsPtITSRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
739 }//fTrackCutsITS loop
741 if(fTrackType==1 || fTrackType==2) delete track;
746 PostData(0, fHistList);
747 PostData(1, fHistListITS);
750 //________________________________________________________________________
751 Bool_t AliPWG4HighPtQAMC::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
753 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
754 // This is to called in Notify and should provide the path to the AOD/ESD file
755 // Copied from AliAnalysisTaskJetSpectrum2
758 TString file(currFile);
762 if(file.Contains("root_archive.zip#")){
763 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
764 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
765 file.Replace(pos+1,20,"");
768 // not an archive take the basename....
769 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
771 // Printf("%s",file.Data());
774 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
776 // next trial fetch the histgram file
777 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
779 // not a severe condition but inciate that we have no information
783 // find the tlist we want to be independtent of the name so use the Tkey
784 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
789 TList *list = dynamic_cast<TList*>(key->ReadObj());
794 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
795 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
798 } // no tree pyxsec.root
800 TTree *xtree = (TTree*)fxsec->Get("Xsection");
806 Double_t xsection = 0;
807 xtree->SetBranchAddress("xsection",&xsection);
808 xtree->SetBranchAddress("ntrials",&ntrials);
816 //________________________________________________________________________
817 Bool_t AliPWG4HighPtQAMC::Notify()
820 // Implemented Notify() to read the cross sections
821 // and number of trials from pyxsec.root
822 // Copied from AliAnalysisTaskJetSpectrum2
825 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
826 Float_t xsection = 0;
831 TFile *curfile = tree->GetCurrentFile();
833 Error("Notify","No current file");
836 if(!fh1Xsec||!fh1Trials){
837 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
840 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
841 fh1Xsec->Fill("<#sigma>",xsection);
842 // construct a poor man average trials
843 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
844 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
849 //________________________________________________________________________
850 AliGenPythiaEventHeader* AliPWG4HighPtQAMC::GetPythiaEventHeader(AliMCEvent *mcEvent){
852 if(!mcEvent)return 0;
853 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
854 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
855 if(!pythiaGenHeader){
857 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
859 if (!genCocktailHeader) {
860 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
861 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
864 TList* headerList = genCocktailHeader->GetHeaders();
865 for (Int_t i=0; i<headerList->GetEntries(); i++) {
866 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
870 if(!pythiaGenHeader){
871 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
875 return pythiaGenHeader;
879 //________________________________________________________________________
880 void AliPWG4HighPtQAMC::Terminate(Option_t *)
882 // The Terminate() function is the last function to be called during
883 // a query. It always runs on the client, it can be used to present
884 // the results graphically or save the results to file.