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(1e6),
86 fPtAllminPtMCvsPtMC(0),
87 fPtAllminPtMCvsPtAll(0),
88 fPtAllvsPtMCvsMult(0),
89 fPtAllminPtMCvsPtAllvsMult(0),
90 fPtAllminPtMCvsPtAllNPointTPC(0),
91 fPtAllminPtMCvsPtAllNPointTPCIter1(0),
92 fPtAllminPtMCvsPtAllChi2TPC(0),
93 fPtAllminPtMCvsPtAllChi2TPCIter1(0),
94 fPtAllminPtMCvsPtAllDCAR(0),
95 fPtAllminPtMCvsPtAllDCAZ(0),
96 fPtAllminPtMCvsPtAllPhi(0),
97 fPtAllminPtMCvsPtAllNPointITS(0),
98 fPtAllminPtMCvsPtAllNSigmaToVertex(0),
99 fPtAllminPtMCvsPtAllChi2C(0),
100 fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
106 fPtBinEdges[0][0] = 10.;
107 fPtBinEdges[0][1] = 1.;
108 fPtBinEdges[1][0] = 20.;
109 fPtBinEdges[1][1] = 2.;
110 fPtBinEdges[2][0] = 100.;
111 fPtBinEdges[2][1] = 10.;
114 //________________________________________________________________________
115 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
116 AliAnalysisTask(name,""),
124 fSigmaConstrainedMax(1e6),
140 fNPointTPCMultRec(0),
143 fPtAllminPtMCvsPtMC(0),
144 fPtAllminPtMCvsPtAll(0),
145 fPtAllvsPtMCvsMult(0),
146 fPtAllminPtMCvsPtAllvsMult(0),
147 fPtAllminPtMCvsPtAllNPointTPC(0),
148 fPtAllminPtMCvsPtAllNPointTPCIter1(0),
149 fPtAllminPtMCvsPtAllChi2TPC(0),
150 fPtAllminPtMCvsPtAllChi2TPCIter1(0),
151 fPtAllminPtMCvsPtAllDCAR(0),
152 fPtAllminPtMCvsPtAllDCAZ(0),
153 fPtAllminPtMCvsPtAllPhi(0),
154 fPtAllminPtMCvsPtAllNPointITS(0),
155 fPtAllminPtMCvsPtAllNSigmaToVertex(0),
156 fPtAllminPtMCvsPtAllChi2C(0),
157 fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
163 // Constructor. Initialization of Inputs and Outputs
165 AliDebug(2,Form("AliPWG4HighPtQAMC Calling Constructor"));
167 fPtBinEdges[0][0] = 10.;
168 fPtBinEdges[0][1] = 1.;
169 fPtBinEdges[1][0] = 20.;
170 fPtBinEdges[1][1] = 2.;
171 fPtBinEdges[2][0] = 100.;
172 fPtBinEdges[2][1] = 10.;
174 // Input slot #0 works with a TChain ESD
175 DefineInput(0, TChain::Class());
176 // Output slot #0, #1 write into a TList
177 DefineOutput(0, TList::Class());
180 //________________________________________________________________________
181 void AliPWG4HighPtQAMC::SetPtBinEdges(Int_t region, Double_t ptmax, Double_t ptBinWidth) {
184 fPtBinEdges[region][0] = ptmax;
185 fPtBinEdges[region][1] = ptBinWidth;
188 AliError("Only 3 regions alowed. Use region 0/1/2\n");
194 //________________________________________________________________________
195 void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
197 // Connect ESD and MC event handler here
199 AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
200 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
202 AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
206 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
209 AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
212 fESD = esdH->GetEvent();
214 AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
216 AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
219 fMC = eventHandler->MCEvent();
224 //________________________________________________________________________
225 void AliPWG4HighPtQAMC::CreateOutputObjects() {
226 //Create output objects
227 AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
229 Bool_t oldStatus = TH1::AddDirectoryStatus();
230 TH1::AddDirectory(kFALSE);
233 fHistList = new TList();
234 fHistList->SetOwner(kTRUE);
236 Float_t fgkPtMin = 0.;
238 //fPtBinEdges[region][0] = ptmax of region ; fPtBinEdges[region][1] = binWidth of region
239 const Float_t ptmin1 = fgkPtMin;
240 const Float_t ptmax1 = fPtBinEdges[0][0];
241 const Float_t ptmin2 = ptmax1 ;
242 const Float_t ptmax2 = fPtBinEdges[1][0];
243 const Float_t ptmin3 = ptmax2 ;
244 const Float_t ptmax3 = fPtBinEdges[2][0];//fgkPtMax;
245 const Int_t nbin11 = (int)((ptmax1-ptmin1)/fPtBinEdges[0][1]);
246 const Int_t nbin12 = (int)((ptmax2-ptmin2)/fPtBinEdges[1][1])+nbin11;
247 const Int_t nbin13 = (int)((ptmax3-ptmin3)/fPtBinEdges[2][1])+nbin12;
248 Int_t fgkNPtBins=nbin13;
249 //Create array with low edges of each bin
250 Double_t *binsPt=new Double_t[fgkNPtBins+1];
251 for(Int_t i=0; i<=fgkNPtBins; i++) {
252 if(i<=nbin11) binsPt[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
253 if(i<=nbin12 && i>nbin11) binsPt[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
254 if(i<=nbin13 && i>nbin12) binsPt[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
257 Int_t fgkNPhiBins = 18*6;
258 Float_t kMinPhi = 0.;
259 Float_t kMaxPhi = 2.*TMath::Pi();
260 Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
261 for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
263 Int_t fgkResPtBins=80;
264 Float_t kMinResPt = -1.;
265 Float_t kMaxResPt = 1.;
266 Double_t *binsResPt = new Double_t[fgkResPtBins+1];
267 for(Int_t i=0; i<=fgkResPtBins; i++) binsResPt[i]=(Double_t)kMinResPt + (kMaxResPt-kMinResPt)/fgkResPtBins*(Double_t)i ;
269 Int_t fgkMultBins=50;
270 Float_t kMinMult = 0.;
271 Float_t kMaxMult = 3000.;
272 Double_t *binsMult = new Double_t[fgkMultBins+1];
273 for(Int_t i=0; i<=fgkMultBins; i++) binsMult[i]=(Double_t)kMinMult + (kMaxMult-kMinMult)/fgkMultBins*(Double_t)i ;
275 Int_t fgkNEtaBins=20;
276 Float_t fgkEtaMin = -1.;
277 Float_t fgkEtaMax = 1.;
278 Double_t *binsEta=new Double_t[fgkNEtaBins+1];
279 for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
281 Int_t fgkNNClustersTPCBins=80;
282 Float_t fgkNClustersTPCMin = 0.5;
283 Float_t fgkNClustersTPCMax = 160.5;
284 Double_t *binsNClustersTPC=new Double_t[fgkNNClustersTPCBins+1];
285 for(Int_t i=0; i<=fgkNNClustersTPCBins; i++) binsNClustersTPC[i]=(Double_t)fgkNClustersTPCMin + (fgkNClustersTPCMax-fgkNClustersTPCMin)/fgkNNClustersTPCBins*(Double_t)i ;
287 Int_t fgkNDCA2DBins=80;
288 Float_t fgkDCA2DMin = -0.2;
289 Float_t fgkDCA2DMax = 0.2;
290 if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
294 Double_t *binsDCA2D=new Double_t[fgkNDCA2DBins+1];
295 for(Int_t i=0; i<=fgkNDCA2DBins; i++) binsDCA2D[i]=(Double_t)fgkDCA2DMin + (fgkDCA2DMax-fgkDCA2DMin)/fgkNDCA2DBins*(Double_t)i ;
297 Int_t fgkNDCAZBins=80;
298 Float_t fgkDCAZMin = -2.;
299 Float_t fgkDCAZMax = 2.;
300 if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
304 Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
305 for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
307 Int_t fgkNNPointITSBins=9;
308 Float_t fgkNPointITSMin = -0.5;
309 Float_t fgkNPointITSMax = 8.5;
310 Double_t *binsNPointITS=new Double_t[fgkNNPointITSBins+1];
311 for(Int_t i=0; i<=fgkNNPointITSBins; i++) binsNPointITS[i]=(Double_t)fgkNPointITSMin + (fgkNPointITSMax-fgkNPointITSMin)/fgkNNPointITSBins*(Double_t)i ;
313 Int_t fgkNNSigmaToVertexBins=20;
314 Float_t fgkNSigmaToVertexMin = 0.;
315 Float_t fgkNSigmaToVertexMax = 8.;
316 Double_t *binsNSigmaToVertex=new Double_t[fgkNNSigmaToVertexBins+1];
317 for(Int_t i=0; i<=fgkNNSigmaToVertexBins; i++) binsNSigmaToVertex[i]=(Double_t)fgkNSigmaToVertexMin + (fgkNSigmaToVertexMax-fgkNSigmaToVertexMin)/fgkNNSigmaToVertexBins*(Double_t)i ;
319 Int_t fgkNChi2CBins=20;
320 Float_t fgkChi2CMin = 0.;
321 Float_t fgkChi2CMax = 100.;
322 Double_t *binsChi2C=new Double_t[fgkNChi2CBins+1];
323 for(Int_t i=0; i<=fgkNChi2CBins; i++) binsChi2C[i]=(Double_t)fgkChi2CMin + (fgkChi2CMax-fgkChi2CMin)/fgkNChi2CBins*(Double_t)i ;
325 Int_t fgkNRel1PtUncertaintyBins=50;
326 Float_t fgkRel1PtUncertaintyMin = 0.;
327 Float_t fgkRel1PtUncertaintyMax = 1.;
328 Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
329 for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
331 Float_t fgkChi2PerClusMin = 0.;
332 Float_t fgkChi2PerClusMax = 4.;
333 Int_t fgkNChi2PerClusBins = (int)(fgkChi2PerClusMax*10.);
334 Double_t *binsChi2PerClus=new Double_t[fgkNChi2PerClusBins+1];
335 for(Int_t i=0; i<=fgkNChi2PerClusBins; i++) binsChi2PerClus[i]=(Double_t)fgkChi2PerClusMin + (fgkChi2PerClusMax-fgkChi2PerClusMin)/fgkNChi2PerClusBins*(Double_t)i ;
338 fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
339 fHistList->Add(fNEventAll);
340 fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
341 fHistList->Add(fNEventSel);
342 fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
344 fNEventReject->Fill("noESD",0);
345 fNEventReject->Fill("Trigger",0);
346 fNEventReject->Fill("noMCEvent",0);
347 fNEventReject->Fill("noStack",0);
348 fNEventReject->Fill("NTracks<2",0);
349 fNEventReject->Fill("noVTX",0);
350 fNEventReject->Fill("VtxStatus",0);
351 fNEventReject->Fill("NCont<2",0);
352 fNEventReject->Fill("ZVTX>10",0);
353 fHistList->Add(fNEventReject);
355 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
356 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
357 fHistList->Add(fh1Xsec);
359 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
360 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
361 fHistList->Add(fh1Trials);
363 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
364 fHistList->Add(fh1PtHard);
365 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
366 fHistList->Add(fh1PtHardTrials);
368 fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
369 fHistList->Add(fPtAll);
370 fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
371 fHistList->Add(fPtSel);
373 fPtSelFakes = new TH1F("fPtSelFakes","PtSelFakes",fgkNPtBins, binsPt);
374 fHistList->Add(fPtSelFakes);
376 fNPointTPCFakes = new TH1F("fNPointTPCFakes","fNPointTPCFakes",fgkNNClustersTPCBins, binsNClustersTPC);
377 fHistList->Add(fNPointTPCFakes);
379 fPtSelLargeLabel = new TH1F("fPtSelLargeLabel","PtSelLargeLabel",fgkNPtBins, binsPt);
380 fHistList->Add(fPtSelLargeLabel);
382 fMultRec = new TH1F("fMultRec","Multiple reconstruction of tracks",fgkMultBins, binsMult);
383 fHistList->Add(fMultRec);
385 fNPointTPCMultRec = new TH1F("fNPointTPCMultRec","Multiple reconstruction of tracks",fgkNNClustersTPCBins, binsNClustersTPC);
386 fHistList->Add(fNPointTPCMultRec);
388 fDeltaPtMultRec = new TH2F("fDeltaPtMultRec","Delta pT vs pT for multiple reconstructed tracks",50,0.,50.,40,-20.,20.);
389 fHistList->Add(fDeltaPtMultRec);
391 fPtAllvsPtMC = new TH2F("fPtAllvsPtMC","fPtAllvsPtMC;p_{T,MC};p_{T,rec}",fgkNPtBins, binsPt,fgkNPtBins, binsPt);
392 fHistList->Add(fPtAllvsPtMC);
394 fPtAllminPtMCvsPtMC = new TH2F("fPtAllminPtMCvsPtMC","PtAllminPtMCvsPtMC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt);
395 fPtAllminPtMCvsPtMC->SetXTitle("p_{t}^{MC}");
396 fPtAllminPtMCvsPtMC->SetYTitle("(1/p_{t}^{rec}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
397 fHistList->Add(fPtAllminPtMCvsPtMC);
399 fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, binsPt,fgkResPtBins,binsResPt);
400 fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{rec}");
401 fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{rec}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
402 fHistList->Add(fPtAllminPtMCvsPtAll);
404 fPtAllvsPtMCvsMult = new TH3F("fPtAllvsPtMCvsMult","fPtAllvsPtMCvsMult;p_{T,MC};p_{T,rec};N_{tracks}",fgkNPtBins, binsPt,fgkNPtBins, binsPt,fgkMultBins,binsMult);
405 fHistList->Add(fPtAllvsPtMCvsMult);
407 fPtAllminPtMCvsPtAllvsMult = new TH3F("fPtAllminPtMCvsPtAllvsMult","fPtAllminPtMCvsPtAllvsMult",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkMultBins,binsMult);
408 fPtAllminPtMCvsPtAllvsMult->SetXTitle("p_{t}^{MC}");
409 fPtAllminPtMCvsPtAllvsMult->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
410 fPtAllminPtMCvsPtAllvsMult->SetZTitle("N_{tracks}");
411 fHistList->Add(fPtAllminPtMCvsPtAllvsMult);
413 fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNClustersTPCBins, binsNClustersTPC);
414 fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
415 fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
416 fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{cls,TPC}");
417 fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
419 fPtAllminPtMCvsPtAllNPointTPCIter1 = new TH3F("fPtAllminPtMCvsPtAllNPointTPCIter1","PtAllminPtMCvsPtAllNPointTPCIter1",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNClustersTPCBins, binsNClustersTPC);
420 fPtAllminPtMCvsPtAllNPointTPCIter1->SetXTitle("p_{t}^{MC}");
421 fPtAllminPtMCvsPtAllNPointTPCIter1->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
422 fPtAllminPtMCvsPtAllNPointTPCIter1->SetZTitle("N_{clsIter1,TPC}");
423 fHistList->Add(fPtAllminPtMCvsPtAllNPointTPCIter1);
425 fPtAllminPtMCvsPtAllChi2TPC = new TH3F("fPtAllminPtMCvsPtAllChi2TPC","PtAllminPtMCvsPtAllChi2TPC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2PerClusBins, binsChi2PerClus);
426 fPtAllminPtMCvsPtAllChi2TPC->SetXTitle("p_{t}^{MC}");
427 fPtAllminPtMCvsPtAllChi2TPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
428 fPtAllminPtMCvsPtAllChi2TPC->SetZTitle("#chi^{2} TPC");
429 fHistList->Add(fPtAllminPtMCvsPtAllChi2TPC);
431 fPtAllminPtMCvsPtAllChi2TPCIter1 = new TH3F("fPtAllminPtMCvsPtAllChi2TPCIter1","PtAllminPtMCvsPtAllChi2TPCIter1",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2PerClusBins, binsChi2PerClus);
432 fPtAllminPtMCvsPtAllChi2TPCIter1->SetXTitle("p_{t}^{MC}");
433 fPtAllminPtMCvsPtAllChi2TPCIter1->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
434 fPtAllminPtMCvsPtAllChi2TPCIter1->SetZTitle("#chi^{2} TPC Iter1");
435 fHistList->Add(fPtAllminPtMCvsPtAllChi2TPCIter1);
437 fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNDCA2DBins,binsDCA2D);
438 fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
439 fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
440 fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
441 fHistList->Add(fPtAllminPtMCvsPtAllDCAR);
443 fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNDCAZBins,binsDCAZ);
444 fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
445 fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
446 fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
447 fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);
449 fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNPhiBins,binsPhi);
450 fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
451 fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
452 fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
453 fHistList->Add(fPtAllminPtMCvsPtAllPhi);
455 fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNPointITSBins,binsNPointITS);
456 fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
457 fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
458 fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
459 fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
461 fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
462 fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
463 fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
464 fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
465 fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);
467 fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2CBins,binsChi2C);
468 fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
469 fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
470 fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
471 fHistList->Add(fPtAllminPtMCvsPtAllChi2C);
473 fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
474 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
475 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
476 fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
477 fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);
480 fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, binsPt);
481 fHistList->Add(fPtAllMC);
482 fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, binsPt);
483 fHistList->Add(fPtSelMC);
485 TH1::AddDirectory(oldStatus);
487 PostData(0, fHistList);
489 if(binsPt) delete [] binsPt;
490 if(binsResPt) delete [] binsResPt;
491 if(binsPhi) delete [] binsPhi;
492 if(binsNClustersTPC) delete [] binsNClustersTPC;
493 if(binsDCA2D) delete [] binsDCA2D;
494 if(binsDCAZ) delete [] binsDCAZ;
495 if(binsNPointITS) delete [] binsNPointITS;
496 if(binsNSigmaToVertex) delete [] binsNSigmaToVertex;
497 if(binsChi2C) delete [] binsChi2C;
498 if(binsEta) delete [] binsEta;
499 if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
500 if(binsChi2PerClus) delete [] binsChi2PerClus;
501 if(binsMult) delete [] binsMult;
505 //________________________________________________________________________
506 Bool_t AliPWG4HighPtQAMC::SelectEvent() {
508 // Decide if event should be selected for analysis
511 // Checks following requirements:
513 // - trigger info from AliPhysicsSelection
514 // - MCevent available
515 // - number of reconstructed tracks > 1
516 // - primary vertex reconstructed
517 // - z-vertex < 10 cm
519 Bool_t selectEvent = kTRUE;
521 //fESD object available?
523 AliDebug(2,Form("ERROR: fInputEvent not available\n"));
524 fNEventReject->Fill("noESD",1);
525 selectEvent = kFALSE;
530 UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
531 if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
532 AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
533 fNEventReject->Fill("Trigger",1);
534 selectEvent = kFALSE;
541 AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
542 fStack = fMC->Stack(); //Particles Stack
544 AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
547 AliDebug(2,Form("ERROR: Could not retrieve MC eventHandler"));
548 fNEventReject->Fill("noStack",1);
549 selectEvent = kFALSE;
553 AliDebug(2,Form("ERROR: Could not retrieve stack"));
554 fNEventReject->Fill("noMCEvent",1);
555 selectEvent = kFALSE;
559 //Check if number of reconstructed tracks is larger than 1
560 if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2) {
561 fNEventReject->Fill("NTracks<2",1);
562 selectEvent = kFALSE;
566 //Check if vertex is reconstructed
567 if(fTrackType==1) fVtx = fESD->GetPrimaryVertexTPC();
568 else fVtx = fESD->GetPrimaryVertexSPD();
571 fNEventReject->Fill("noVTX",1);
572 selectEvent = kFALSE;
576 if(!fVtx->GetStatus()) {
577 fNEventReject->Fill("VtxStatus",1);
578 selectEvent = kFALSE;
583 // TString vtxName(fVtx->GetName());
584 if(fVtx->GetNContributors()<2) {
585 fNEventReject->Fill("NCont<2",1);
586 selectEvent = kFALSE;
590 //Check if z-vertex < 10 cm
592 fVtx->GetXYZ(primVtx);
593 if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
594 fNEventReject->Fill("ZVTX>10",1);
595 selectEvent = kFALSE;
599 AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",fVtx->GetTitle(), fVtx->GetStatus(), fVtx->GetNContributors()));
605 //________________________________________________________________________
606 void AliPWG4HighPtQAMC::Exec(Option_t *) {
608 // Called for each event
609 AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
610 // All events without selection
611 fNEventAll->Fill(0.);
615 PostData(0, fHistList);
619 // ---- Get MC Header information (for MC productions in pThard bins) ----
620 Double_t ptHard = 0.;
621 Double_t nTrials = 1; // trials for MC trigger weight for real data
624 AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(fMC);
626 nTrials = pythiaGenHeader->Trials();
627 ptHard = pythiaGenHeader->GetPtHard();
629 fh1PtHard->Fill(ptHard);
630 fh1PtHardTrials->Fill(ptHard,nTrials);
632 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
636 //Need to keep track of selected events
637 fNEventSel->Fill(0.);
639 Int_t nTracks = fESD->GetNumberOfTracks();
640 AliDebug(2,Form("nTracks ESD%d", nTracks));
642 int nMCtracks = fStack->GetNtrack();
651 Float_t nSigmaToVertex = 0.;
652 Float_t relUncertainty1Pt = 0.;
654 int mult = fTrackCuts->CountAcceptedTracks(fESD);
656 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
658 AliESDtrack *track = 0;
659 AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
660 if(!esdtrack) continue;
663 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
666 else if(fTrackType==2) {
667 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
670 AliExternalTrackParam exParam;
671 Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
676 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
678 else if(fTrackType==7) {
679 //use global constrained track
680 track = new AliESDtrack(*esdtrack);
682 // track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
690 //Cut on chi2 of constrained fit
691 if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
697 Int_t label = TMath::Abs(track->GetLabel());
698 if(label>=nMCtracks) {
699 if (fTrackCuts->AcceptTrack(track)) { fPtSelLargeLabel->Fill(pt); }
700 if(fTrackType==1 || fTrackType==2 || fTrackType==7) delete track;
704 TParticle *particle = fStack->Particle(label) ;
706 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
707 if(track) delete track;
712 ptMC = particle->Pt();
714 if(fTrackType==1 || fTrackType==2)
715 track->GetImpactParametersTPC(dca2D,dcaZ); //TPConly
717 track->GetImpactParameters(dca2D,dcaZ); //Global
719 UChar_t itsMap = track->GetITSClusterMap();
720 for (Int_t i=0; i < 6; i++) {
721 if (itsMap & (1 << i))
726 fPtAllMC->Fill(ptMC);
728 if (fTrackCuts->AcceptTrack(track)) {
731 if(fTrackCutsReject ) {
732 if(fTrackCutsReject->AcceptTrack(track) ) {
733 if(track) delete track;
738 if(esdtrack->GetConstrainedParam())
739 track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
745 nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
746 chi2C = track->GetConstrainedChi2();
747 relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
750 if(track->GetLabel()<0) {
751 fPtSelFakes->Fill(pt);
752 fNPointTPCFakes->Fill(track->GetTPCNcls());
754 fPtSelMC->Fill(ptMC);
755 fPtAllvsPtMC->Fill(ptMC,pt);
756 fPtAllminPtMCvsPtMC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
757 fPtAllminPtMCvsPtAll->Fill(pt,(1./pt-1./ptMC)/(1./ptMC) );
758 fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
759 fPtAllminPtMCvsPtAllNPointTPCIter1->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNclsIter1());
760 if(track->GetTPCNcls()>0.)
761 fPtAllminPtMCvsPtAllChi2TPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCchi2()/track->GetTPCNcls());
762 if(track->GetTPCNclsIter1()>0.)
763 fPtAllminPtMCvsPtAllChi2TPCIter1->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCchi2Iter1()/track->GetTPCNclsIter1());
764 fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
765 fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
766 fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
767 fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
768 fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
769 fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
770 fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
772 fPtAllvsPtMCvsMult->Fill(ptMC,pt,mult);
773 fPtAllminPtMCvsPtAllvsMult->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC), mult);
775 //Check if track is reconstructed multiple times
778 for (Int_t iTrack2 = iTrack+1; iTrack2 < nTracks; iTrack2++) {
779 // AliESDtrack *track2 = GetTrackForAnalysis(iTrack2);
781 AliESDtrack *esdtrack2 = fESD->GetTrack(iTrack2);
782 if(!esdtrack2) continue;
785 track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
786 else if(fTrackType==2) {
787 track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
791 AliExternalTrackParam exParam2;
792 Bool_t relate = track2->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam2);
797 track2->Set(exParam2.GetX(),exParam2.GetAlpha(),exParam2.GetParameter(),exParam2.GetCovariance());
805 if (fTrackCuts->AcceptTrack(track2)) {
806 Int_t label2 = TMath::Abs(track2->GetLabel());
808 fNPointTPCMultRec->Fill(track->GetTPCNcls());
809 fNPointTPCMultRec->Fill(track2->GetTPCNcls());
810 fDeltaPtMultRec->Fill(track->Pt(),track->Pt()-track2->Pt());
814 if(fTrackType==1 || fTrackType==2) delete track2;
816 if(multCounter>1) fMultRec->Fill(multCounter);
819 }//fTrackCuts selection
822 if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
823 if(track) delete track;
829 PostData(0, fHistList);
832 //________________________________________________________________________
833 Bool_t AliPWG4HighPtQAMC::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
835 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
836 // This is to called in Notify and should provide the path to the AOD/ESD file
837 // Copied from AliAnalysisTaskJetSpectrum2
840 TString file(currFile);
844 if(file.Contains("root_archive.zip#")){
845 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
846 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
847 file.Replace(pos+1,20,"");
850 // not an archive take the basename....
851 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
853 // Printf("%s",file.Data());
856 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
858 // next trial fetch the histgram file
859 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
861 // not a severe condition but inciate that we have no information
865 // find the tlist we want to be independtent of the name so use the Tkey
866 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
871 TList *list = dynamic_cast<TList*>(key->ReadObj());
876 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
877 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
880 } // no tree pyxsec.root
882 TTree *xtree = (TTree*)fxsec->Get("Xsection");
888 Double_t xsection = 0;
889 xtree->SetBranchAddress("xsection",&xsection);
890 xtree->SetBranchAddress("ntrials",&ntrials);
898 //________________________________________________________________________
899 Bool_t AliPWG4HighPtQAMC::Notify()
902 // Implemented Notify() to read the cross sections
903 // and number of trials from pyxsec.root
904 // Copied from AliAnalysisTaskJetSpectrum2
907 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
908 Float_t xsection = 0;
913 TFile *curfile = tree->GetCurrentFile();
915 Error("Notify","No current file");
918 if(!fh1Xsec||!fh1Trials){
919 // Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
922 PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
923 fh1Xsec->Fill("<#sigma>",xsection);
924 // construct a poor man average trials
925 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
926 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
931 //________________________________________________________________________
932 AliGenPythiaEventHeader* AliPWG4HighPtQAMC::GetPythiaEventHeader(AliMCEvent *mcEvent){
934 if(!mcEvent)return 0;
935 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
936 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
937 if(!pythiaGenHeader){
939 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
941 if (!genCocktailHeader) {
942 // AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
943 // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
946 TList* headerList = genCocktailHeader->GetHeaders();
947 for (Int_t i=0; i<headerList->GetEntries(); i++) {
948 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
952 if(!pythiaGenHeader){
953 AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
957 return pythiaGenHeader;
961 //________________________________________________________________________
962 void AliPWG4HighPtQAMC::Terminate(Option_t *)
964 // The Terminate() function is the last function to be called during
965 // a query. It always runs on the client, it can be used to present
966 // the results graphically or save the results to file.