Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskSEHFCJqa.cxx
CommitLineData
e8b1c5e0 1/**************************************************************************
2 * Copyright(c) 1998-2009, 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//
18// Class AliAnalysisTaskSEHFCJqa
19// AliAnalysisTaskSE for the extraction of the fraction of prompt charm
20// using the charm hadron impact parameter to the primary vertex
21//
22// Author: Andrea Rossi, andrea.rossi@pd.infn.it
23/////////////////////////////////////////////////////////////
24
25
26#include <TH1F.h>
27#include <TH2F.h>
28#include <TH3F.h>
29#include <TAxis.h>
30#include <THnSparse.h>
31#include <TDatabasePDG.h>
32#include <TMath.h>
33#include <TROOT.h>
34#include "AliAODEvent.h"
35#include "AliAODRecoDecayHF2Prong.h"
36#include "AliAODRecoDecayHF.h"
37#include "AliAODRecoDecay.h"
38#include "AliAnalysisDataSlot.h"
39#include "AliAnalysisDataContainer.h"
40#include "AliAODTrack.h"
41#include "AliAODHandler.h"
42#include "AliESDtrack.h"
43#include "AliAODVertex.h"
44#include "AliESDVertex.h"
45#include "AliVertexerTracks.h"
46#include "AliAODMCParticle.h"
47#include "AliAODPid.h"
48#include "AliTPCPIDResponse.h"
49#include "AliAODMCHeader.h"
50#include "AliAnalysisVertexingHF.h"
51#include "AliAnalysisTaskSEHFCJqa.h"
52#include "AliRDHFCutsD0toKpi.h"
53#include "AliAODInputHandler.h"
54#include "AliAnalysisManager.h"
55#include "AliNormalizationCounter.h"
56
57class TCanvas;
58class TTree;
59class TChain;
60class AliAnalysisTaskSE;
61
62
63ClassImp(AliAnalysisTaskSEHFCJqa)
64
65 AliAnalysisTaskSEHFCJqa::AliAnalysisTaskSEHFCJqa()
66: AliAnalysisTaskSE(),
67 fReadMC(),
68 ffilterbit(),
69 fKeepTrackNegID(),
70 fpidResp(),
71 fCuts(),
72 fhEventCounter(),
73 fhImpParResolITSsel(),
74 fhImpParResolITSselGoodTracks(),
75 fhSparseFilterMask(),
76 fhSparseFilterMaskTrackAcc(),
77 fhSparseFilterMaskImpPar(),
78 fhSparseEoverPeleTPC(),
79 fhSparseShowShapeEleTPC(),
80 fhnSigmaTPCTOFEle(),
81 fhnSigmaTPCTOFPion(),
82 fhnSigmaTPCTOFKaon(),
83 fhnSigmaTPCTOFProton(),
84 fhTrackEMCal(),
85 fSparseRecoJets(),
86 fLoadJet(),
87 fJetArrayString(),
88 fListTrackAndPID(),
89 fListJets()
90{// standard constructor
91
92}
93
94
95
96//________________________________________________________________________
97AliAnalysisTaskSEHFCJqa::AliAnalysisTaskSEHFCJqa(const char *name)
98 : AliAnalysisTaskSE(name),
99 fReadMC(),
100 ffilterbit(AliAODTrack::kTrkGlobalNoDCA),
101 fKeepTrackNegID(),
102 fpidResp(),
103 fCuts(),
104 fhEventCounter(),
105 fhImpParResolITSsel(),
106 fhImpParResolITSselGoodTracks(),
107 fhSparseFilterMask(),
108 fhSparseFilterMaskTrackAcc(),
109 fhSparseFilterMaskImpPar(),
110 fhSparseEoverPeleTPC(),
111 fhSparseShowShapeEleTPC(),
112 fhnSigmaTPCTOFEle(),
113 fhnSigmaTPCTOFPion(),
114 fhnSigmaTPCTOFKaon(),
115 fhnSigmaTPCTOFProton(),
116 fhTrackEMCal(),
117 fSparseRecoJets(),
118 fLoadJet(),
119 fJetArrayString(),
120 fListTrackAndPID(),
121 fListJets()
122{ // default constructor
123
124 DefineOutput(1, TH1F::Class()); // counter
125 DefineOutput(2, TList::Class()); // single track properties list and PID
126 DefineOutput(3, TList::Class()); // jet properties list
127
128}
129
130//________________________________________________________________________
131AliAnalysisTaskSEHFCJqa::~AliAnalysisTaskSEHFCJqa(){
132 // destructor
133 delete fCuts;
134 delete fpidResp;
135 delete fhEventCounter;
136 delete fhSparseFilterMask;
137 delete fhSparseFilterMaskTrackAcc;
138 delete fhSparseFilterMaskImpPar;
139 delete fhSparseEoverPeleTPC;
140 delete fhSparseShowShapeEleTPC;
141 delete fhnSigmaTPCTOFEle;
142 delete fhnSigmaTPCTOFPion;
143 delete fhnSigmaTPCTOFKaon;
144 delete fhnSigmaTPCTOFProton;
145 delete fhTrackEMCal;
146 delete fSparseRecoJets;
147 delete fhImpParResolITSsel;
148 delete fhImpParResolITSselGoodTracks;
149 delete fListTrackAndPID;
150 delete fListJets;
151}
152
153//________________________________________________________________________
154void AliAnalysisTaskSEHFCJqa::Init()
155{
156 // Initialization
157
158}
159
160
161//________________________________________________________________________
162void AliAnalysisTaskSEHFCJqa::UserCreateOutputObjects(){
163
164
165 //########## DEFINE THE TLISTS ##################
166 fListTrackAndPID=new TList();
167 fListTrackAndPID->SetOwner();
168 fListTrackAndPID->SetName("fListTrackAndPID");
169
170 fListJets=new TList();
171 fListJets->SetOwner();
172 fListJets->SetName("fListJets");
173
174
175 // ########### DEFINE THE EVENT COUNTER ############
176 fhEventCounter=new TH1F("fhEventCounter","Counter of event selected",20,-0.5,19.5);
177 fhEventCounter->GetXaxis()->SetBinLabel(1,"Events analyzed");
178 fhEventCounter->GetXaxis()->SetBinLabel(2,"Event selected");
179 fhEventCounter->GetXaxis()->SetBinLabel(3,"Jet array present");
180 fhEventCounter->GetXaxis()->SetBinLabel(4,"Vtx Track Ncont");
181
182
183
184 Int_t nbinsRecoJets[8]={50,20,20,20,5,5,10,60};
185 Double_t binlowRecoJets[8]={5.,-1.,-TMath::Pi(),0.99,0.,-0.5,0,0.};
186 Double_t binupRecoJets[8]={55.,1.,TMath::Pi(),20.99,4.99,4.5,2.,60.};
187 fSparseRecoJets=new THnSparseF("fSparseRecoJets","fSparseRecoJets;jetpt;eta;phi;ntrks;nEle;parton;partContr;ptPart;",8,nbinsRecoJets,binlowRecoJets,binupRecoJets);
188
189 fListJets->Add(fSparseRecoJets);
190
191 // Num axes: 10 filter bits + ID + TPCrefit,ITSrefit,bothTPCandITSrefit + kAny,kFirst,kBoth+ 20Nclust TPC+ 20 NTPC crossed padRows + 20 DCA
192 Int_t nbinsFilterMask[16]={2,2,2,2,2,2,2,2,2,2,2,3,3,20,20,70};
193 Double_t binlowFilterMask[16]={-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,0,0,-3.5};
194 Double_t binupFilterMask[16]={1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,2.5,2.5,160,160,3.5};
195 fhSparseFilterMask=new THnSparseF("fhSparseFilterMask","fhSparseFilterMask;fb0;fb1;fb2;fb3;fb4;fb5;fb6;fb7;fb8;fb9;trackID;refitting;SPD;NTPCclust;NTPCcrossRows;DCA",16,nbinsFilterMask,binlowFilterMask,binupFilterMask);
196 fListTrackAndPID->Add(fhSparseFilterMask);
197
198
199// Num axes: 10 filter bits + ID*isSelected 5 + kAny,kFirst,kBoth + phi+ eta +pt
200 Int_t nbinsFilterMaskTrackAcc[15]={2,2,2,2,2,2,2,2,2,2,5,3,36,30,30};
201 Double_t binlowFilterMaskTrackAcc[15]={-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-2.5,-0.5,0.,-1.5,0.};
202 Double_t binupFilterMaskTrackAcc[15]={1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,2.5,2.5,TMath::Pi()*2.,1.5,15.};
203 fhSparseFilterMaskTrackAcc=new THnSparseF("fhSparseFilterMaskTrackAcc","fhSparseFilterMaskTrackAcc;fb0;fb1;fb2;fb3;fb4;fb5;fb6;fb7;fb8;fb9;trackIDandsel;SPD;#phi (rad);#eta;p_{T} (GeV/c)",15,nbinsFilterMaskTrackAcc,binlowFilterMaskTrackAcc,binupFilterMaskTrackAcc);
204 fListTrackAndPID->Add(fhSparseFilterMaskTrackAcc);
205
206
207// Num axes: ID*isSelected 5 + kAny,kFirst,kBoth + phi+ eta +pt + imp par
208 Int_t nbinsFilterMaskImpPar[6]={5,3,36,30,30,200};
209 Double_t binlowFilterMaskImpPar[6]={-2.5,-0.5,0.,-1.5,0.,-300.};
210 Double_t binupFilterMaskImpPar[6]={2.5,2.5,TMath::Pi()*2.,1.5,15.,300.};
211 fhSparseFilterMaskImpPar=new THnSparseF("fhSparseFilterMaskImpPar","fhSparseFilterMaskImpPar;fb0;fb1;fb2;fb3;fb4;fb5;fb6;fb7;fb8;fb9;trackIDandsel;SPD;#phi (rad);#eta;p_{T} (GeV/c; imp par (#mum)",6,nbinsFilterMaskImpPar,binlowFilterMaskImpPar,binupFilterMaskImpPar);
212 fListTrackAndPID->Add(fhSparseFilterMaskImpPar);
213
214
215fhImpParResolITSsel=new TH3F("fhImpParResolITSsel","fhImpParResolITSsel;p_{T} (GeV/c);imp. par (#mum);ITS clust",50.,0.,10.,200,-800.,800.,38,-0.5,37.5);
216 // the convention for ITS clust axis: number between 0 and 48, first digit (units) = number of clust in ITS, second digits (x10): SPD status: 0 = none, 1=kFirst, 2=kSecond; 3= kBoth
217 fListTrackAndPID->Add(fhImpParResolITSsel);
218
219 fhImpParResolITSselGoodTracks=new TH3F("fhImpParResolITSselGoodTracks","fhImpParResolITSselGoodTracks;p_{T} (GeV/c);imp. par (#mum);ITS clust",50.,0.,10.,200,-800.,800.,38,-0.5,37.5);
220 // the convention for ITS clust axis: number between 0 and 48, first digit (units) = number of clust in ITS, second digits (x10): SPD status: 0 = none, 1=kFirst, 2=kSecond; 3= kBoth
221 fListTrackAndPID->Add(fhImpParResolITSselGoodTracks);
222
223
224 // PID PLOTS
225 fhnSigmaTPCTOFEle=new TH3F("fhnSigmaTPCTOFEle","fhnSigmaTPCTOFEle;p (GeV/c);n#sigma_{TPC}^{ele};n#sigma_{TOF}^{ele}",120,0.,30,64,-8.,8.,64,-8.,8.);
226 fhnSigmaTPCTOFPion=new TH3F("fhnSigmaTPCTOFPion","fhnSigmaTPCTOFPion;p (GeV/c);n#sigma_{TPC}^{ele};n#sigma_{TOF}^{ele}",120,0.,30,64,-8.,8.,64,-8.,8.);
227 fhnSigmaTPCTOFKaon=new TH3F("fhnSigmaTPCTOFKaon","fhnSigmaTPCTOFKaon;p (GeV/c);n#sigma_{TPC}^{ele};n#sigma_{TOF}^{ele}",120,0.,30,64,-8.,8.,64,-8.,8.);
228 fhnSigmaTPCTOFProton=new TH3F("fhnSigmaTPCTOFProton","fhnSigmaTPCTOFProton;p (GeV/c);n#sigma_{TPC}^{ele};n#sigma_{TOF}^{ele}",120,0.,30,64,-8.,8.,64,-8.,8.);
229
230 fListTrackAndPID->Add(fhnSigmaTPCTOFEle);
231 fListTrackAndPID->Add(fhnSigmaTPCTOFPion);
232 fListTrackAndPID->Add(fhnSigmaTPCTOFKaon);
233 fListTrackAndPID->Add(fhnSigmaTPCTOFProton);
234
235
236 /// EMCAL PLOTS
237 //study of NsigmaTPC, e/p, p
238 Int_t nbinsEoP[4]={202,400,100,400};
239 Double_t binlowEoP[4]= {-1., -20.,-20., -1};
240 Double_t binupEoP[4]= {100., 20, 20.,9};
241 fhSparseEoverPeleTPC = new THnSparseF("fhSparseEoP", "fhSparseEoP; p;nsigmatpc;nsigmaElePIDresp;E/p;",4, nbinsEoP, binlowEoP, binupEoP);
242 fListTrackAndPID->Add(fhSparseEoverPeleTPC);
243 //fSparseEoverPallHadr = new THnSparseF("fSparseEoP", "fSparseEoP; p;nsigmatpc;E/p;",3, nbinsEoP, binlowEoP, binupEoP);
244
245 Int_t nbinsEmShSh[7]={35,120,100,50,50,50,50};
246 Double_t binlowEmShSh[7]= {5., -20., -1,0,0,0,0};
247 Double_t binupEmShSh[7]= {40., 10, 9,1,1,2,50};
248 fhSparseShowShapeEleTPC = new THnSparseF("fhSparseShowShapeEleTPC", "fhSparseShowShapeEleTPC; pt;nsigmatpc;E/p;M02;M20;disp;Nclust",7, nbinsEmShSh, binlowEmShSh, binupEmShSh);
249 fListTrackAndPID->Add(fhSparseShowShapeEleTPC);
250
251
252
253 Int_t nbinsTrEM[6]={124,20,124,25,20,20};
254 Double_t binlowTrEM[6]={-1,-1.,-1.0,0.,0.};
255 Double_t binupTrEM[6]={31,1.,31.,5.,1.,1.};
256
257 fhTrackEMCal=new THnSparseF("fhTrackEMCal","fhTrackEMCal;p(GeV/c);eta;clusterE(GeV/c);E/p;trkDistX;trkDistZ",6,nbinsTrEM,binlowTrEM,binupTrEM);
258
259
260
261 //fhPhotonicEle=new TH3F("fhPhotonicEle","fhPhotonicEle; pele;pgamma;mass",20,0.,20.,20,0.,20.,50,0.,5);
262 //fhPhotonicEleLS=new TH3F("fhPhotonicEleLS","fhPhotonicEleLS; pele;pgamma;mass",20,0.,20.,20,0.,20.,50,0.,5);
263
264
265 PostData(1,fhEventCounter);
266 PostData(2,fListTrackAndPID);
267 PostData(3,fListJets);
268
269}
270
271
272
273//________________________________________________________________________
274void AliAnalysisTaskSEHFCJqa::UserExec(Option_t */*option*/){
275
276
277
278// Execute analysis for current event:
279 // heavy flavor candidates association to MC truth
280
281 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
7ae7961c 282 if(aod){
283 fhEventCounter->Fill(0);
e8b1c5e0 284
7ae7961c 285 if(!fCuts->IsEventSelected(aod)){
286 PostData(1,fhEventCounter);
287 return;
288 }
289 fhEventCounter->Fill(1);
e8b1c5e0 290 }
e8b1c5e0 291
7ae7961c 292 TClonesArray *arrayJets=0x0;
e8b1c5e0 293 if(!aod && AODEvent() && IsStandardAOD()) {
294 // In case there is an AOD handler writing a standard AOD, use the AOD
295 // event in memory rather than the input (ESD) event.
296 aod = dynamic_cast<AliAODEvent*> (AODEvent());
297 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
298 // have to taken from the AOD event hold by the AliAODExtension
7ae7961c 299 if (!aod) {
300 Printf("ERROR: aod not available");
301 return;
302 }
303 else {
304 fhEventCounter->Fill(0);
305
306 if(!fCuts->IsEventSelected(aod)){
307 PostData(1,fhEventCounter);
308 return;
309 }
310 fhEventCounter->Fill(1);
311 }
312
e8b1c5e0 313 AliAODHandler* aodHandler = (AliAODHandler*)
314 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
315
316 if(fLoadJet>=1&&aodHandler->GetExtensions()) {
317 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.Jets.root");
318 AliAODEvent* aodFromExt = ext->GetAOD();
319
320 // load jet array
321 arrayJets = (TClonesArray*)aodFromExt->GetList()->FindObject(fJetArrayString.Data());
322 if(!arrayJets) {
323 Printf("AliAnalysisTaskSEHFCJqa::UserExec: desired branch %s not found!",fJetArrayString.Data());
324 PostData(1,fhEventCounter);
325 return;
326
327 }
328
329 }
330 } else {
331 // load jet array
332 if(fLoadJet>=1){
333 arrayJets = (TClonesArray*)aod->GetList()->FindObject(fJetArrayString.Data());
334 if(!arrayJets) {
335 Printf("AliAnalysisTaskSEHFCJqa::UserExec: desired Jet branch %s not found!",fJetArrayString.Data());
336 PostData(1,fhEventCounter);
337 return;
338 }
339
340 }
341 }
342
343 if(fLoadJet!=0 && !arrayJets) {
344 printf("AliAnalysisTaskSEHFCJqa::UserExec: desired jet input branch not found!\n");
345 PostData(1,fhEventCounter);
346 return;
347 }
348
349 fhEventCounter->Fill(2);
350
351 // AOD primary vertex
352 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
353 TString primTitle = vtx1->GetTitle();
354 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
355 fhEventCounter->Fill(3);
356
357 }
358 else {
359 PostData(1,fhEventCounter);
360 return;
361
362 }
363 // Convert primary vertex to esd vertex (needed for later usage of RelateToVertex)
364 Double_t pos[3],cov[6];
365 vtx1->GetXYZ(pos);
366 vtx1->GetCovarianceMatrix(cov);
367 const AliESDVertex vESD(pos,cov,100.,100);
368 Double_t magfield=aod->GetMagneticField();
369
370 TClonesArray *arrayMC=0x0;
371 AliAODMCHeader *aodmcHeader=0x0;
372 Double_t vtxTrue[3];
373
374
375 if(fReadMC){
376 // load MC particles
377 arrayMC =
378 (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
379 if(!arrayMC) {
380 Printf("AliAnalysisTaskSEHFCJqa::UserExec: MC particles branch not found!\n");
381 return;
382 }
383 // load MC header
384 aodmcHeader =
385 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
386 if(!aodmcHeader) {
387 Printf("AliAnalysisTaskSEHFCJqa::UserExec: MC header branch not found!\n");
388 return;
389 }
390 // MC primary vertex
391 aodmcHeader->GetVertex(vtxTrue);
392
393 }
394
395 // Starting the fun part
396 SetupPIDresponse();// this sets the pid reponse to fPPIDRespons; could also get from the cut object (AliAODPidHF::GetPidResponse)
397
398 // Looping over aod tracks
399 for(Int_t j=0;j<aod->GetNTracks();j++){
400
401 AliAODTrack *aodtrack=aod->GetTrack(j);
402 // CHECK FILTER MAPS
7ae7961c 403 if(!FillTrackHistosAndSelectTrack(aodtrack,&vESD,magfield))continue;
e8b1c5e0 404 // if(j%100==0)
405
406 Double_t p=aodtrack->P();
407 Double_t eta=aodtrack->Eta();
408 // START PID: TPC
409
410 Double_t nsigmaEleTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kElectron);
411 Double_t nsigmaPionTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kPion);
412 Double_t nsigmaKaonTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kKaon);
413 Double_t nsigmaProtonTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kProton);
414
415
416 // TOF
417 Double_t nsigmaEleTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kElectron);
418 Double_t nsigmaPionTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kPion);
419 Double_t nsigmaKaonTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kKaon);
420 Double_t nsigmaProtonTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kProton);
421
422
423 fhnSigmaTPCTOFEle->Fill(p,nsigmaEleTPC,nsigmaEleTOF);
424 fhnSigmaTPCTOFPion->Fill(p,nsigmaPionTPC,nsigmaPionTOF);
425 fhnSigmaTPCTOFKaon->Fill(p,nsigmaKaonTPC,nsigmaKaonTOF);
426 fhnSigmaTPCTOFProton->Fill(p,nsigmaProtonTPC,nsigmaProtonTOF);
427 // if(j%100==0)
428
429
430 // NOW EMCAL
431 // CHECK WHETHER THERE IS A EMCAL CLUSTER
432 Int_t nClsId = aodtrack->GetEMCALcluster();
433 if(nClsId >0) {
434 AliVCluster *cluster = aod->GetCaloCluster(nClsId);
435 Double_t clsE = cluster->E();
436 if(!cluster->IsEMCAL()) continue;
437
438 // Check whether it is an electron candidate: do not reject here hadrons to allow QA checks of (E/p,nsigmaTPC,pt)
439 Double_t nEoverP = clsE/p;
440 Double_t eOverPpidResp;
441 Double_t showerShape[4];
442 Double_t nsigmaEleEMCal=fpidResp->NumberOfSigmasEMCAL(aodtrack,AliPID::kElectron,eOverPpidResp,showerShape);
443 Double_t poix[4]={p, nsigmaEleTPC, nsigmaEleEMCal, nEoverP};
444 fhSparseEoverPeleTPC->Fill(poix);
445
446
447 Double_t emcTrackDx=cluster->GetTrackDx();
448 Double_t emcTrackDz=cluster->GetTrackDz();
449 Double_t pointET[6]={p,eta,clsE,nEoverP,emcTrackDx,emcTrackDz};
450 fhTrackEMCal->Fill(pointET);
451
2942f542 452 Double_t pointEmShSh[7]={aodtrack->Pt(), static_cast<Double_t>(nsigmaEleTPC), static_cast<Double_t>(nEoverP),cluster->GetM02(),cluster->GetM20(),cluster->GetDispersion(), static_cast<Double_t>(cluster->GetNCells())};
e8b1c5e0 453
454 fhSparseShowShapeEleTPC->Fill(pointEmShSh);
455
456 }
457 }
458
459
460 // NOW LOOP OVER JETS
461
462 Int_t nJets=arrayJets->GetEntries();//calcolo numero di jet nell'event
463 AliAODJet *jet;
464 for(Int_t jetcand=0;jetcand<nJets;jetcand++){
465 // if(jetcand%100==0)
466
467 jet=(AliAODJet*)arrayJets->UncheckedAt(jetcand);
468
469 Double_t contribution=0,ptpart=-1;
470 Int_t partonnat=0;
471 if(fReadMC){
472
473 AliAODMCParticle *parton=IsMCJet(arrayMC,jet,contribution);
474 if(parton){
475 Int_t pdg=TMath::Abs(parton->PdgCode());
476 //printf("pdg parton: %d \n",pdg);
477 if(pdg==21)partonnat=1;
478 else if(pdg<4)partonnat=2;
479 else if(pdg==4)partonnat=3;
480 else if(pdg==5)partonnat=4;
481 ptpart=parton->Pt();
482 }
483
484 }
485
486
487 FillJetRecoHisto(jet,partonnat,contribution,ptpart);
488 }
489
490
491
492 PostData(1,fhEventCounter);
493 PostData(2,fListTrackAndPID);
494 PostData(3,fListJets);
495
496 return;
497}
498
499
500
501void AliAnalysisTaskSEHFCJqa::SetupPIDresponse(){
502
503 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
504 AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
505 fpidResp=inputHandler->GetPIDResponse();
506 if(!fpidResp)AliFatal("No PID response could be set");
507}
508
509//_______________________________________________________________
7ae7961c 510Bool_t AliAnalysisTaskSEHFCJqa::FillTrackHistosAndSelectTrack(AliAODTrack *aodtr, const AliESDVertex *primary, const Double_t magfield){
e8b1c5e0 511
512 Bool_t retval=kTRUE;
513 // THnSparse for filter bits
514 Double_t point[16]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,999.,-999.,-999.,-999.,-999.,-999.};
515 Double_t pointAcc[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,999.,-999.,-999.,-999.,-999.};
516 Double_t pointImpPar[6]={999.,-999.,-999.,-999.,-999.,-999.};
517
518 for(Int_t j=0;j<10;j++){
519 if(aodtr->TestFilterBit(TMath::Power(2,j))){
520 point[j]=1;
521 pointAcc[j]=1;
522 }
523 }
524 // check ID
525 Int_t trid=aodtr->GetID();
526 if(aodtr->GetID()>0)point[10]=1.;
527 else if(aodtr->GetID()>0)point[10]=0.;
528 if(aodtr->GetID()==0)point[10]=-1.;
529
530 Float_t iparxy,iparz;
531 Float_t pt=aodtr->Pt();
532 Int_t clustITS=aodtr->GetITSNcls();// for histo
533 AliESDtrack esdtrack(aodtr);
534 // needed to calculate impact parameters
535
536
537 // check refit status
538 Int_t refit=-1;
539 UInt_t status = esdtrack.GetStatus();
540 if(status&AliESDtrack::kTPCrefit)refit+=1;
541 if(status&AliESDtrack::kITSrefit)refit+=2;
542 point[11]=refit;
543 // CHECK SPD
544 point[12]=-1;
545 if(aodtr->HasPointOnITSLayer(0)){
546 clustITS+=10;
547 point[12]+=1;
548 }
549 if(aodtr->HasPointOnITSLayer(1)){
550 point[12]+=2;
551 clustITS+=20;
552 }
553 point[13]=aodtr->GetTPCNcls();
554 point[14]=aodtr->GetTPCNCrossedRows();
7ae7961c 555 esdtrack.RelateToVertex(primary,magfield,4.);// CHECK THIS : I put 4.. usually we set it to 3
e8b1c5e0 556 esdtrack.GetImpactParameters(iparxy,iparz);
557 point[15]=iparxy;
558
559
560
561
562 fhSparseFilterMask->Fill(point);
563 if(!aodtr->TestBit(ffilterbit))retval =kFALSE;
564
565 if(aodtr->GetID()<0&&!fKeepTrackNegID)retval = kFALSE;
566 if(retval) fhImpParResolITSsel->Fill(pt,iparxy*10000.,clustITS);
567
568 AliESDtrackCuts *cuts=fCuts->GetTrackCuts();
569 if(!cuts->IsSelected(&esdtrack)) retval = kFALSE;
570
571 if(fCuts->GetUseKinkRejection()){
572 AliAODVertex *maybeKink=aodtr->GetProdVertex();
573 if(maybeKink->GetType()==AliAODVertex::kKink) retval=kFALSE;
574 }
575
576 if(retval){
577 pointAcc[10]=1*trid;
578 }
579 else {
580 if(trid!=0)
581 pointAcc[10]=2*trid;
582 else pointAcc[10]=-3;
583 }
584
585 pointAcc[11]=point[12];
586 pointAcc[12]=aodtr->Phi();
587 if(pointAcc[12]<0.)pointAcc[12]+=2.*TMath::Pi();
588 pointAcc[13]=aodtr->Eta();
589 pointAcc[14]=pt;
590 fhSparseFilterMaskTrackAcc->Fill(pointAcc);
591
592 pointImpPar[0]=pointAcc[10];
593 pointImpPar[1]=pointAcc[11];
594 pointImpPar[2]=pointAcc[12];
595 pointImpPar[3]=pointAcc[13];
596 pointImpPar[4]=pointAcc[14];
597 pointImpPar[5]=iparxy;
598 fhSparseFilterMaskImpPar->Fill(pointImpPar);
599
600 if(retval) {
601 fhImpParResolITSselGoodTracks->Fill(pt,iparxy*10000.,clustITS);
602 }
603
604 return retval;
605
606}
607
608
609//---------------------------------------------------------------
610AliAODMCParticle* AliAnalysisTaskSEHFCJqa::IsMCJet(TClonesArray *arrayMC,const AliAODJet *jet, Double_t &contribution){// assignment of parton ID to jet
611 // method by L. Feldkamp
612 std::vector< int > idx;
613 std::vector< int > idx2;
614 std::vector< double > weight;
615
616 int counter =0;
617 int num = jet->GetRefTracks()->GetEntries();
618
619 for(int k=0;k<num;++k){
620
621 AliAODTrack * track = (AliAODTrack*)(jet->GetRefTracks()->At(k));
622 if (track->GetLabel() >=0){
623 AliAODMCParticle* part = (AliAODMCParticle*) arrayMC->At(track->GetLabel());
624 if(!part)continue;
625
626 int label =0 ;
627 AliAODMCParticle* motherParton=GetMCPartonOrigin(arrayMC,part, label);
628 if (!motherParton) Printf("no mother");
629 else {
630 counter++;
631 idx.push_back(label); //! Label of Mother
632 idx2.push_back(label);
633 weight.push_back(track->Pt()); //! Weight : P_t trak / P_t jet ... the latter used at the end
634 }
635 }///END LOOP OVER REFERENCED TRACKS
636 }
637 //! Remove duplicate indices for counting
638 sort( idx2.begin(), idx2.end() );
639 idx2.erase( unique( idx2.begin(), idx2.end() ), idx2.end() );
640 if (idx2.size() == 0) return 0x0;
7ae7961c 641 Double_t* arrayOfWeights = new Double_t[(UInt_t)idx2.size()];
642 if(!arrayOfWeights){
643 return 0x0;
644 }
645 for(UInt_t ii=0;ii<(UInt_t)idx2.size();ii++)arrayOfWeights[ii]=0;
e8b1c5e0 646
647 for (unsigned int idxloop =0 ;idxloop<idx2.size();idxloop++){
648 for (unsigned int z=0; z< idx.size() ; ++z){
649 int a = idx.at(z);
650 double w = weight.at(z);
651 if(a == idx2.at(idxloop)) arrayOfWeights[idxloop] += w;
652 }
653 }
654
655 int winner = -1;
656 double c=-1.;
657 for (unsigned int idxloop =0 ;idxloop<idx2.size();idxloop++){
658 if(c < arrayOfWeights[idxloop]){
659 winner =idxloop;
660 c=arrayOfWeights[idxloop];
661 }
662 }
663
664 AliAODMCParticle *parton = 0x0;
7ae7961c 665 if(winner>0){
666 parton=(AliAODMCParticle*)arrayMC->At(idx.at(winner));
667 contribution = arrayOfWeights[winner]/jet->Pt();
668 }
669 else {
670
671 if(arrayOfWeights) delete [] arrayOfWeights;
672 return 0x0;
673 }
674 if(arrayOfWeights) delete [] arrayOfWeights;
e8b1c5e0 675
676 return parton;
677
678
679}
680
681//---------------------------------------------------------------
682AliAODMCParticle *AliAnalysisTaskSEHFCJqa::GetMCPartonOrigin(TClonesArray *arrayMC,AliAODMCParticle *p, Int_t &idx)
683{ //Follows chain of track mothers until q/g or idx = -1
684 AliAODMCParticle *p2=0x0;
685 Int_t mlabel = TMath::Abs(p->GetMother()) ;
686 Double_t pz=0.;
687 while(mlabel > 1){
688 p2 = (AliAODMCParticle*)arrayMC->At(mlabel);
689 pz=TMath::Abs(p2->Pz());
690 //printf("Mother label %d, pdg %d, pz %f\n",mlabel,p2->PdgCode(),pz);
691 if( p2->PdgCode() == 21 || (p2->PdgCode() != 0 && abs(p2->PdgCode()) <6) )
692 {
693 idx = mlabel;
694 return p2;
695 }
696 mlabel = TMath::Abs(p2->GetMother());
697 }
698 idx=-1;
699 return 0x0;
700
701}
702
703//_____________________________________________________________
704void AliAnalysisTaskSEHFCJqa::FillJetRecoHisto(const AliAODJet *jet,Int_t partonnat,Double_t contribution,Double_t ptpart){
705//FIll sparse with reco jets properties
2942f542 706 Double_t point[8]={jet->Pt(),jet->Eta(),jet->Phi()-TMath::Pi(), static_cast<Double_t>(jet->GetRefTracks()->GetEntriesFast()),0, static_cast<Double_t>(partonnat),contribution,ptpart};
e8b1c5e0 707 fSparseRecoJets->Fill(point);
708}
709
710
711//_______________________________________________________________
712//void AliAnalysisTaskSEHFCJqa::FillTrackHistosPID(AliAODTrack *aodtr){
713//
714//}
715
716//_______________________________________________________________
717void AliAnalysisTaskSEHFCJqa::Terminate(const Option_t*){
718 //TERMINATE METHOD: NOTHING TO DO
719
720
721
722}
723