]>
Commit | Line | Data |
---|---|---|
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 | ||
57 | class TCanvas; | |
58 | class TTree; | |
59 | class TChain; | |
60 | class AliAnalysisTaskSE; | |
61 | ||
62 | ||
63 | ClassImp(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 | //________________________________________________________________________ | |
97 | AliAnalysisTaskSEHFCJqa::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 | //________________________________________________________________________ | |
131 | AliAnalysisTaskSEHFCJqa::~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 | //________________________________________________________________________ | |
154 | void AliAnalysisTaskSEHFCJqa::Init() | |
155 | { | |
156 | // Initialization | |
157 | ||
158 | } | |
159 | ||
160 | ||
161 | //________________________________________________________________________ | |
162 | void 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 | ||
215 | fhImpParResolITSsel=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 | //________________________________________________________________________ | |
274 | void 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 | |
4640c275 | 399 | for(Int_t j=0;j<aod->GetNumberOfTracks();j++){ |
e8b1c5e0 | 400 | |
f15c1f69 | 401 | AliAODTrack *aodtrack=dynamic_cast<AliAODTrack*>(aod->GetTrack(j)); |
402 | if(!aodtrack) AliFatal("Not a standard AOD"); | |
e8b1c5e0 | 403 | // CHECK FILTER MAPS |
7ae7961c | 404 | if(!FillTrackHistosAndSelectTrack(aodtrack,&vESD,magfield))continue; |
e8b1c5e0 | 405 | // if(j%100==0) |
406 | ||
407 | Double_t p=aodtrack->P(); | |
408 | Double_t eta=aodtrack->Eta(); | |
409 | // START PID: TPC | |
410 | ||
411 | Double_t nsigmaEleTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kElectron); | |
412 | Double_t nsigmaPionTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kPion); | |
413 | Double_t nsigmaKaonTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kKaon); | |
414 | Double_t nsigmaProtonTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kProton); | |
415 | ||
416 | ||
417 | // TOF | |
418 | Double_t nsigmaEleTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kElectron); | |
419 | Double_t nsigmaPionTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kPion); | |
420 | Double_t nsigmaKaonTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kKaon); | |
421 | Double_t nsigmaProtonTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kProton); | |
422 | ||
423 | ||
424 | fhnSigmaTPCTOFEle->Fill(p,nsigmaEleTPC,nsigmaEleTOF); | |
425 | fhnSigmaTPCTOFPion->Fill(p,nsigmaPionTPC,nsigmaPionTOF); | |
426 | fhnSigmaTPCTOFKaon->Fill(p,nsigmaKaonTPC,nsigmaKaonTOF); | |
427 | fhnSigmaTPCTOFProton->Fill(p,nsigmaProtonTPC,nsigmaProtonTOF); | |
428 | // if(j%100==0) | |
429 | ||
430 | ||
431 | // NOW EMCAL | |
432 | // CHECK WHETHER THERE IS A EMCAL CLUSTER | |
433 | Int_t nClsId = aodtrack->GetEMCALcluster(); | |
434 | if(nClsId >0) { | |
435 | AliVCluster *cluster = aod->GetCaloCluster(nClsId); | |
436 | Double_t clsE = cluster->E(); | |
437 | if(!cluster->IsEMCAL()) continue; | |
438 | ||
439 | // Check whether it is an electron candidate: do not reject here hadrons to allow QA checks of (E/p,nsigmaTPC,pt) | |
440 | Double_t nEoverP = clsE/p; | |
441 | Double_t eOverPpidResp; | |
442 | Double_t showerShape[4]; | |
443 | Double_t nsigmaEleEMCal=fpidResp->NumberOfSigmasEMCAL(aodtrack,AliPID::kElectron,eOverPpidResp,showerShape); | |
444 | Double_t poix[4]={p, nsigmaEleTPC, nsigmaEleEMCal, nEoverP}; | |
445 | fhSparseEoverPeleTPC->Fill(poix); | |
446 | ||
447 | ||
448 | Double_t emcTrackDx=cluster->GetTrackDx(); | |
449 | Double_t emcTrackDz=cluster->GetTrackDz(); | |
450 | Double_t pointET[6]={p,eta,clsE,nEoverP,emcTrackDx,emcTrackDz}; | |
451 | fhTrackEMCal->Fill(pointET); | |
452 | ||
2942f542 | 453 | 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 | 454 | |
455 | fhSparseShowShapeEleTPC->Fill(pointEmShSh); | |
456 | ||
457 | } | |
458 | } | |
459 | ||
460 | ||
461 | // NOW LOOP OVER JETS | |
462 | ||
463 | Int_t nJets=arrayJets->GetEntries();//calcolo numero di jet nell'event | |
464 | AliAODJet *jet; | |
465 | for(Int_t jetcand=0;jetcand<nJets;jetcand++){ | |
466 | // if(jetcand%100==0) | |
467 | ||
468 | jet=(AliAODJet*)arrayJets->UncheckedAt(jetcand); | |
469 | ||
470 | Double_t contribution=0,ptpart=-1; | |
471 | Int_t partonnat=0; | |
472 | if(fReadMC){ | |
473 | ||
474 | AliAODMCParticle *parton=IsMCJet(arrayMC,jet,contribution); | |
475 | if(parton){ | |
476 | Int_t pdg=TMath::Abs(parton->PdgCode()); | |
477 | //printf("pdg parton: %d \n",pdg); | |
478 | if(pdg==21)partonnat=1; | |
479 | else if(pdg<4)partonnat=2; | |
480 | else if(pdg==4)partonnat=3; | |
481 | else if(pdg==5)partonnat=4; | |
482 | ptpart=parton->Pt(); | |
483 | } | |
484 | ||
485 | } | |
486 | ||
487 | ||
488 | FillJetRecoHisto(jet,partonnat,contribution,ptpart); | |
489 | } | |
490 | ||
491 | ||
492 | ||
493 | PostData(1,fhEventCounter); | |
494 | PostData(2,fListTrackAndPID); | |
495 | PostData(3,fListJets); | |
496 | ||
497 | return; | |
498 | } | |
499 | ||
500 | ||
501 | ||
502 | void AliAnalysisTaskSEHFCJqa::SetupPIDresponse(){ | |
503 | ||
504 | AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); | |
505 | AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler(); | |
506 | fpidResp=inputHandler->GetPIDResponse(); | |
507 | if(!fpidResp)AliFatal("No PID response could be set"); | |
508 | } | |
509 | ||
510 | //_______________________________________________________________ | |
7ae7961c | 511 | Bool_t AliAnalysisTaskSEHFCJqa::FillTrackHistosAndSelectTrack(AliAODTrack *aodtr, const AliESDVertex *primary, const Double_t magfield){ |
e8b1c5e0 | 512 | |
513 | Bool_t retval=kTRUE; | |
514 | // THnSparse for filter bits | |
515 | Double_t point[16]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,999.,-999.,-999.,-999.,-999.,-999.}; | |
516 | Double_t pointAcc[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,999.,-999.,-999.,-999.,-999.}; | |
517 | Double_t pointImpPar[6]={999.,-999.,-999.,-999.,-999.,-999.}; | |
518 | ||
519 | for(Int_t j=0;j<10;j++){ | |
520 | if(aodtr->TestFilterBit(TMath::Power(2,j))){ | |
521 | point[j]=1; | |
522 | pointAcc[j]=1; | |
523 | } | |
524 | } | |
525 | // check ID | |
526 | Int_t trid=aodtr->GetID(); | |
527 | if(aodtr->GetID()>0)point[10]=1.; | |
528 | else if(aodtr->GetID()>0)point[10]=0.; | |
529 | if(aodtr->GetID()==0)point[10]=-1.; | |
530 | ||
531 | Float_t iparxy,iparz; | |
532 | Float_t pt=aodtr->Pt(); | |
533 | Int_t clustITS=aodtr->GetITSNcls();// for histo | |
534 | AliESDtrack esdtrack(aodtr); | |
535 | // needed to calculate impact parameters | |
536 | ||
537 | ||
538 | // check refit status | |
539 | Int_t refit=-1; | |
540 | UInt_t status = esdtrack.GetStatus(); | |
541 | if(status&AliESDtrack::kTPCrefit)refit+=1; | |
542 | if(status&AliESDtrack::kITSrefit)refit+=2; | |
543 | point[11]=refit; | |
544 | // CHECK SPD | |
545 | point[12]=-1; | |
546 | if(aodtr->HasPointOnITSLayer(0)){ | |
547 | clustITS+=10; | |
548 | point[12]+=1; | |
549 | } | |
550 | if(aodtr->HasPointOnITSLayer(1)){ | |
551 | point[12]+=2; | |
552 | clustITS+=20; | |
553 | } | |
554 | point[13]=aodtr->GetTPCNcls(); | |
555 | point[14]=aodtr->GetTPCNCrossedRows(); | |
7ae7961c | 556 | esdtrack.RelateToVertex(primary,magfield,4.);// CHECK THIS : I put 4.. usually we set it to 3 |
e8b1c5e0 | 557 | esdtrack.GetImpactParameters(iparxy,iparz); |
558 | point[15]=iparxy; | |
559 | ||
560 | ||
561 | ||
562 | ||
563 | fhSparseFilterMask->Fill(point); | |
564 | if(!aodtr->TestBit(ffilterbit))retval =kFALSE; | |
565 | ||
566 | if(aodtr->GetID()<0&&!fKeepTrackNegID)retval = kFALSE; | |
567 | if(retval) fhImpParResolITSsel->Fill(pt,iparxy*10000.,clustITS); | |
568 | ||
569 | AliESDtrackCuts *cuts=fCuts->GetTrackCuts(); | |
570 | if(!cuts->IsSelected(&esdtrack)) retval = kFALSE; | |
571 | ||
572 | if(fCuts->GetUseKinkRejection()){ | |
573 | AliAODVertex *maybeKink=aodtr->GetProdVertex(); | |
574 | if(maybeKink->GetType()==AliAODVertex::kKink) retval=kFALSE; | |
575 | } | |
576 | ||
577 | if(retval){ | |
578 | pointAcc[10]=1*trid; | |
579 | } | |
580 | else { | |
581 | if(trid!=0) | |
582 | pointAcc[10]=2*trid; | |
583 | else pointAcc[10]=-3; | |
584 | } | |
585 | ||
586 | pointAcc[11]=point[12]; | |
587 | pointAcc[12]=aodtr->Phi(); | |
588 | if(pointAcc[12]<0.)pointAcc[12]+=2.*TMath::Pi(); | |
589 | pointAcc[13]=aodtr->Eta(); | |
590 | pointAcc[14]=pt; | |
591 | fhSparseFilterMaskTrackAcc->Fill(pointAcc); | |
592 | ||
593 | pointImpPar[0]=pointAcc[10]; | |
594 | pointImpPar[1]=pointAcc[11]; | |
595 | pointImpPar[2]=pointAcc[12]; | |
596 | pointImpPar[3]=pointAcc[13]; | |
597 | pointImpPar[4]=pointAcc[14]; | |
598 | pointImpPar[5]=iparxy; | |
599 | fhSparseFilterMaskImpPar->Fill(pointImpPar); | |
600 | ||
601 | if(retval) { | |
602 | fhImpParResolITSselGoodTracks->Fill(pt,iparxy*10000.,clustITS); | |
603 | } | |
604 | ||
605 | return retval; | |
606 | ||
607 | } | |
608 | ||
609 | ||
610 | //--------------------------------------------------------------- | |
611 | AliAODMCParticle* AliAnalysisTaskSEHFCJqa::IsMCJet(TClonesArray *arrayMC,const AliAODJet *jet, Double_t &contribution){// assignment of parton ID to jet | |
612 | // method by L. Feldkamp | |
613 | std::vector< int > idx; | |
614 | std::vector< int > idx2; | |
615 | std::vector< double > weight; | |
616 | ||
617 | int counter =0; | |
618 | int num = jet->GetRefTracks()->GetEntries(); | |
619 | ||
620 | for(int k=0;k<num;++k){ | |
621 | ||
622 | AliAODTrack * track = (AliAODTrack*)(jet->GetRefTracks()->At(k)); | |
623 | if (track->GetLabel() >=0){ | |
624 | AliAODMCParticle* part = (AliAODMCParticle*) arrayMC->At(track->GetLabel()); | |
625 | if(!part)continue; | |
626 | ||
627 | int label =0 ; | |
628 | AliAODMCParticle* motherParton=GetMCPartonOrigin(arrayMC,part, label); | |
629 | if (!motherParton) Printf("no mother"); | |
630 | else { | |
631 | counter++; | |
632 | idx.push_back(label); //! Label of Mother | |
633 | idx2.push_back(label); | |
634 | weight.push_back(track->Pt()); //! Weight : P_t trak / P_t jet ... the latter used at the end | |
635 | } | |
636 | }///END LOOP OVER REFERENCED TRACKS | |
637 | } | |
638 | //! Remove duplicate indices for counting | |
639 | sort( idx2.begin(), idx2.end() ); | |
640 | idx2.erase( unique( idx2.begin(), idx2.end() ), idx2.end() ); | |
641 | if (idx2.size() == 0) return 0x0; | |
7ae7961c | 642 | Double_t* arrayOfWeights = new Double_t[(UInt_t)idx2.size()]; |
643 | if(!arrayOfWeights){ | |
644 | return 0x0; | |
645 | } | |
646 | for(UInt_t ii=0;ii<(UInt_t)idx2.size();ii++)arrayOfWeights[ii]=0; | |
e8b1c5e0 | 647 | |
648 | for (unsigned int idxloop =0 ;idxloop<idx2.size();idxloop++){ | |
649 | for (unsigned int z=0; z< idx.size() ; ++z){ | |
650 | int a = idx.at(z); | |
651 | double w = weight.at(z); | |
652 | if(a == idx2.at(idxloop)) arrayOfWeights[idxloop] += w; | |
653 | } | |
654 | } | |
655 | ||
656 | int winner = -1; | |
657 | double c=-1.; | |
658 | for (unsigned int idxloop =0 ;idxloop<idx2.size();idxloop++){ | |
659 | if(c < arrayOfWeights[idxloop]){ | |
660 | winner =idxloop; | |
661 | c=arrayOfWeights[idxloop]; | |
662 | } | |
663 | } | |
664 | ||
665 | AliAODMCParticle *parton = 0x0; | |
7ae7961c | 666 | if(winner>0){ |
667 | parton=(AliAODMCParticle*)arrayMC->At(idx.at(winner)); | |
668 | contribution = arrayOfWeights[winner]/jet->Pt(); | |
669 | } | |
670 | else { | |
671 | ||
672 | if(arrayOfWeights) delete [] arrayOfWeights; | |
673 | return 0x0; | |
674 | } | |
675 | if(arrayOfWeights) delete [] arrayOfWeights; | |
e8b1c5e0 | 676 | |
677 | return parton; | |
678 | ||
679 | ||
680 | } | |
681 | ||
682 | //--------------------------------------------------------------- | |
683 | AliAODMCParticle *AliAnalysisTaskSEHFCJqa::GetMCPartonOrigin(TClonesArray *arrayMC,AliAODMCParticle *p, Int_t &idx) | |
684 | { //Follows chain of track mothers until q/g or idx = -1 | |
685 | AliAODMCParticle *p2=0x0; | |
686 | Int_t mlabel = TMath::Abs(p->GetMother()) ; | |
687 | Double_t pz=0.; | |
688 | while(mlabel > 1){ | |
689 | p2 = (AliAODMCParticle*)arrayMC->At(mlabel); | |
690 | pz=TMath::Abs(p2->Pz()); | |
691 | //printf("Mother label %d, pdg %d, pz %f\n",mlabel,p2->PdgCode(),pz); | |
692 | if( p2->PdgCode() == 21 || (p2->PdgCode() != 0 && abs(p2->PdgCode()) <6) ) | |
693 | { | |
694 | idx = mlabel; | |
695 | return p2; | |
696 | } | |
697 | mlabel = TMath::Abs(p2->GetMother()); | |
698 | } | |
699 | idx=-1; | |
700 | return 0x0; | |
701 | ||
702 | } | |
703 | ||
704 | //_____________________________________________________________ | |
705 | void AliAnalysisTaskSEHFCJqa::FillJetRecoHisto(const AliAODJet *jet,Int_t partonnat,Double_t contribution,Double_t ptpart){ | |
706 | //FIll sparse with reco jets properties | |
2942f542 | 707 | 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 | 708 | fSparseRecoJets->Fill(point); |
709 | } | |
710 | ||
711 | ||
712 | //_______________________________________________________________ | |
713 | //void AliAnalysisTaskSEHFCJqa::FillTrackHistosPID(AliAODTrack *aodtr){ | |
714 | // | |
715 | //} | |
716 | ||
717 | //_______________________________________________________________ | |
718 | void AliAnalysisTaskSEHFCJqa::Terminate(const Option_t*){ | |
719 | //TERMINATE METHOD: NOTHING TO DO | |
720 | ||
721 | ||
722 | ||
723 | } | |
724 |