]>
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()); | |
282 | if (!aod) { | |
283 | Printf("ERROR: aod not available"); | |
284 | return; | |
285 | } | |
286 | fhEventCounter->Fill(0); | |
287 | ||
288 | if(!fCuts->IsEventSelected(aod)){ | |
289 | PostData(1,fhEventCounter); | |
290 | return; | |
291 | } | |
292 | fhEventCounter->Fill(1); | |
293 | ||
294 | TClonesArray *arrayJets; | |
295 | if(!aod && AODEvent() && IsStandardAOD()) { | |
296 | // In case there is an AOD handler writing a standard AOD, use the AOD | |
297 | // event in memory rather than the input (ESD) event. | |
298 | aod = dynamic_cast<AliAODEvent*> (AODEvent()); | |
299 | // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) | |
300 | // have to taken from the AOD event hold by the AliAODExtension | |
301 | AliAODHandler* aodHandler = (AliAODHandler*) | |
302 | ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); | |
303 | ||
304 | if(fLoadJet>=1&&aodHandler->GetExtensions()) { | |
305 | AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.Jets.root"); | |
306 | AliAODEvent* aodFromExt = ext->GetAOD(); | |
307 | ||
308 | // load jet array | |
309 | arrayJets = (TClonesArray*)aodFromExt->GetList()->FindObject(fJetArrayString.Data()); | |
310 | if(!arrayJets) { | |
311 | Printf("AliAnalysisTaskSEHFCJqa::UserExec: desired branch %s not found!",fJetArrayString.Data()); | |
312 | PostData(1,fhEventCounter); | |
313 | return; | |
314 | ||
315 | } | |
316 | ||
317 | } | |
318 | } else { | |
319 | // load jet array | |
320 | if(fLoadJet>=1){ | |
321 | arrayJets = (TClonesArray*)aod->GetList()->FindObject(fJetArrayString.Data()); | |
322 | if(!arrayJets) { | |
323 | Printf("AliAnalysisTaskSEHFCJqa::UserExec: desired Jet branch %s not found!",fJetArrayString.Data()); | |
324 | PostData(1,fhEventCounter); | |
325 | return; | |
326 | } | |
327 | ||
328 | } | |
329 | } | |
330 | ||
331 | if(fLoadJet!=0 && !arrayJets) { | |
332 | printf("AliAnalysisTaskSEHFCJqa::UserExec: desired jet input branch not found!\n"); | |
333 | PostData(1,fhEventCounter); | |
334 | return; | |
335 | } | |
336 | ||
337 | fhEventCounter->Fill(2); | |
338 | ||
339 | // AOD primary vertex | |
340 | AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex(); | |
341 | TString primTitle = vtx1->GetTitle(); | |
342 | if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) { | |
343 | fhEventCounter->Fill(3); | |
344 | ||
345 | } | |
346 | else { | |
347 | PostData(1,fhEventCounter); | |
348 | return; | |
349 | ||
350 | } | |
351 | // Convert primary vertex to esd vertex (needed for later usage of RelateToVertex) | |
352 | Double_t pos[3],cov[6]; | |
353 | vtx1->GetXYZ(pos); | |
354 | vtx1->GetCovarianceMatrix(cov); | |
355 | const AliESDVertex vESD(pos,cov,100.,100); | |
356 | Double_t magfield=aod->GetMagneticField(); | |
357 | ||
358 | TClonesArray *arrayMC=0x0; | |
359 | AliAODMCHeader *aodmcHeader=0x0; | |
360 | Double_t vtxTrue[3]; | |
361 | ||
362 | ||
363 | if(fReadMC){ | |
364 | // load MC particles | |
365 | arrayMC = | |
366 | (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()); | |
367 | if(!arrayMC) { | |
368 | Printf("AliAnalysisTaskSEHFCJqa::UserExec: MC particles branch not found!\n"); | |
369 | return; | |
370 | } | |
371 | // load MC header | |
372 | aodmcHeader = | |
373 | (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()); | |
374 | if(!aodmcHeader) { | |
375 | Printf("AliAnalysisTaskSEHFCJqa::UserExec: MC header branch not found!\n"); | |
376 | return; | |
377 | } | |
378 | // MC primary vertex | |
379 | aodmcHeader->GetVertex(vtxTrue); | |
380 | ||
381 | } | |
382 | ||
383 | // Starting the fun part | |
384 | SetupPIDresponse();// this sets the pid reponse to fPPIDRespons; could also get from the cut object (AliAODPidHF::GetPidResponse) | |
385 | ||
386 | // Looping over aod tracks | |
387 | for(Int_t j=0;j<aod->GetNTracks();j++){ | |
388 | ||
389 | AliAODTrack *aodtrack=aod->GetTrack(j); | |
390 | // CHECK FILTER MAPS | |
391 | if(!FillTrackHistosAndSelectTrack(aodtrack,vESD,magfield))continue; | |
392 | // if(j%100==0) | |
393 | ||
394 | Double_t p=aodtrack->P(); | |
395 | Double_t eta=aodtrack->Eta(); | |
396 | // START PID: TPC | |
397 | ||
398 | Double_t nsigmaEleTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kElectron); | |
399 | Double_t nsigmaPionTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kPion); | |
400 | Double_t nsigmaKaonTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kKaon); | |
401 | Double_t nsigmaProtonTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kProton); | |
402 | ||
403 | ||
404 | // TOF | |
405 | Double_t nsigmaEleTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kElectron); | |
406 | Double_t nsigmaPionTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kPion); | |
407 | Double_t nsigmaKaonTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kKaon); | |
408 | Double_t nsigmaProtonTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kProton); | |
409 | ||
410 | ||
411 | fhnSigmaTPCTOFEle->Fill(p,nsigmaEleTPC,nsigmaEleTOF); | |
412 | fhnSigmaTPCTOFPion->Fill(p,nsigmaPionTPC,nsigmaPionTOF); | |
413 | fhnSigmaTPCTOFKaon->Fill(p,nsigmaKaonTPC,nsigmaKaonTOF); | |
414 | fhnSigmaTPCTOFProton->Fill(p,nsigmaProtonTPC,nsigmaProtonTOF); | |
415 | // if(j%100==0) | |
416 | ||
417 | ||
418 | // NOW EMCAL | |
419 | // CHECK WHETHER THERE IS A EMCAL CLUSTER | |
420 | Int_t nClsId = aodtrack->GetEMCALcluster(); | |
421 | if(nClsId >0) { | |
422 | AliVCluster *cluster = aod->GetCaloCluster(nClsId); | |
423 | Double_t clsE = cluster->E(); | |
424 | if(!cluster->IsEMCAL()) continue; | |
425 | ||
426 | // Check whether it is an electron candidate: do not reject here hadrons to allow QA checks of (E/p,nsigmaTPC,pt) | |
427 | Double_t nEoverP = clsE/p; | |
428 | Double_t eOverPpidResp; | |
429 | Double_t showerShape[4]; | |
430 | Double_t nsigmaEleEMCal=fpidResp->NumberOfSigmasEMCAL(aodtrack,AliPID::kElectron,eOverPpidResp,showerShape); | |
431 | Double_t poix[4]={p, nsigmaEleTPC, nsigmaEleEMCal, nEoverP}; | |
432 | fhSparseEoverPeleTPC->Fill(poix); | |
433 | ||
434 | ||
435 | Double_t emcTrackDx=cluster->GetTrackDx(); | |
436 | Double_t emcTrackDz=cluster->GetTrackDz(); | |
437 | Double_t pointET[6]={p,eta,clsE,nEoverP,emcTrackDx,emcTrackDz}; | |
438 | fhTrackEMCal->Fill(pointET); | |
439 | ||
440 | Double_t pointEmShSh[7]={aodtrack->Pt(),nsigmaEleTPC,nEoverP,cluster->GetM02(),cluster->GetM20(),cluster->GetDispersion(),cluster->GetNCells()}; | |
441 | ||
442 | fhSparseShowShapeEleTPC->Fill(pointEmShSh); | |
443 | ||
444 | } | |
445 | } | |
446 | ||
447 | ||
448 | // NOW LOOP OVER JETS | |
449 | ||
450 | Int_t nJets=arrayJets->GetEntries();//calcolo numero di jet nell'event | |
451 | AliAODJet *jet; | |
452 | for(Int_t jetcand=0;jetcand<nJets;jetcand++){ | |
453 | // if(jetcand%100==0) | |
454 | ||
455 | jet=(AliAODJet*)arrayJets->UncheckedAt(jetcand); | |
456 | ||
457 | Double_t contribution=0,ptpart=-1; | |
458 | Int_t partonnat=0; | |
459 | if(fReadMC){ | |
460 | ||
461 | AliAODMCParticle *parton=IsMCJet(arrayMC,jet,contribution); | |
462 | if(parton){ | |
463 | Int_t pdg=TMath::Abs(parton->PdgCode()); | |
464 | //printf("pdg parton: %d \n",pdg); | |
465 | if(pdg==21)partonnat=1; | |
466 | else if(pdg<4)partonnat=2; | |
467 | else if(pdg==4)partonnat=3; | |
468 | else if(pdg==5)partonnat=4; | |
469 | ptpart=parton->Pt(); | |
470 | } | |
471 | ||
472 | } | |
473 | ||
474 | ||
475 | FillJetRecoHisto(jet,partonnat,contribution,ptpart); | |
476 | } | |
477 | ||
478 | ||
479 | ||
480 | PostData(1,fhEventCounter); | |
481 | PostData(2,fListTrackAndPID); | |
482 | PostData(3,fListJets); | |
483 | ||
484 | return; | |
485 | } | |
486 | ||
487 | ||
488 | ||
489 | void AliAnalysisTaskSEHFCJqa::SetupPIDresponse(){ | |
490 | ||
491 | AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); | |
492 | AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler(); | |
493 | fpidResp=inputHandler->GetPIDResponse(); | |
494 | if(!fpidResp)AliFatal("No PID response could be set"); | |
495 | } | |
496 | ||
497 | //_______________________________________________________________ | |
498 | Bool_t AliAnalysisTaskSEHFCJqa::FillTrackHistosAndSelectTrack(AliAODTrack *aodtr, const AliESDVertex primary, const Double_t magfield){ | |
499 | ||
500 | Bool_t retval=kTRUE; | |
501 | // THnSparse for filter bits | |
502 | Double_t point[16]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,999.,-999.,-999.,-999.,-999.,-999.}; | |
503 | Double_t pointAcc[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,999.,-999.,-999.,-999.,-999.}; | |
504 | Double_t pointImpPar[6]={999.,-999.,-999.,-999.,-999.,-999.}; | |
505 | ||
506 | for(Int_t j=0;j<10;j++){ | |
507 | if(aodtr->TestFilterBit(TMath::Power(2,j))){ | |
508 | point[j]=1; | |
509 | pointAcc[j]=1; | |
510 | } | |
511 | } | |
512 | // check ID | |
513 | Int_t trid=aodtr->GetID(); | |
514 | if(aodtr->GetID()>0)point[10]=1.; | |
515 | else if(aodtr->GetID()>0)point[10]=0.; | |
516 | if(aodtr->GetID()==0)point[10]=-1.; | |
517 | ||
518 | Float_t iparxy,iparz; | |
519 | Float_t pt=aodtr->Pt(); | |
520 | Int_t clustITS=aodtr->GetITSNcls();// for histo | |
521 | AliESDtrack esdtrack(aodtr); | |
522 | // needed to calculate impact parameters | |
523 | ||
524 | ||
525 | // check refit status | |
526 | Int_t refit=-1; | |
527 | UInt_t status = esdtrack.GetStatus(); | |
528 | if(status&AliESDtrack::kTPCrefit)refit+=1; | |
529 | if(status&AliESDtrack::kITSrefit)refit+=2; | |
530 | point[11]=refit; | |
531 | // CHECK SPD | |
532 | point[12]=-1; | |
533 | if(aodtr->HasPointOnITSLayer(0)){ | |
534 | clustITS+=10; | |
535 | point[12]+=1; | |
536 | } | |
537 | if(aodtr->HasPointOnITSLayer(1)){ | |
538 | point[12]+=2; | |
539 | clustITS+=20; | |
540 | } | |
541 | point[13]=aodtr->GetTPCNcls(); | |
542 | point[14]=aodtr->GetTPCNCrossedRows(); | |
543 | esdtrack.RelateToVertex(&primary,magfield,4.);// CHECK THIS : I put 4.. usually we set it to 3 | |
544 | esdtrack.GetImpactParameters(iparxy,iparz); | |
545 | point[15]=iparxy; | |
546 | ||
547 | ||
548 | ||
549 | ||
550 | fhSparseFilterMask->Fill(point); | |
551 | if(!aodtr->TestBit(ffilterbit))retval =kFALSE; | |
552 | ||
553 | if(aodtr->GetID()<0&&!fKeepTrackNegID)retval = kFALSE; | |
554 | if(retval) fhImpParResolITSsel->Fill(pt,iparxy*10000.,clustITS); | |
555 | ||
556 | AliESDtrackCuts *cuts=fCuts->GetTrackCuts(); | |
557 | if(!cuts->IsSelected(&esdtrack)) retval = kFALSE; | |
558 | ||
559 | if(fCuts->GetUseKinkRejection()){ | |
560 | AliAODVertex *maybeKink=aodtr->GetProdVertex(); | |
561 | if(maybeKink->GetType()==AliAODVertex::kKink) retval=kFALSE; | |
562 | } | |
563 | ||
564 | if(retval){ | |
565 | pointAcc[10]=1*trid; | |
566 | } | |
567 | else { | |
568 | if(trid!=0) | |
569 | pointAcc[10]=2*trid; | |
570 | else pointAcc[10]=-3; | |
571 | } | |
572 | ||
573 | pointAcc[11]=point[12]; | |
574 | pointAcc[12]=aodtr->Phi(); | |
575 | if(pointAcc[12]<0.)pointAcc[12]+=2.*TMath::Pi(); | |
576 | pointAcc[13]=aodtr->Eta(); | |
577 | pointAcc[14]=pt; | |
578 | fhSparseFilterMaskTrackAcc->Fill(pointAcc); | |
579 | ||
580 | pointImpPar[0]=pointAcc[10]; | |
581 | pointImpPar[1]=pointAcc[11]; | |
582 | pointImpPar[2]=pointAcc[12]; | |
583 | pointImpPar[3]=pointAcc[13]; | |
584 | pointImpPar[4]=pointAcc[14]; | |
585 | pointImpPar[5]=iparxy; | |
586 | fhSparseFilterMaskImpPar->Fill(pointImpPar); | |
587 | ||
588 | if(retval) { | |
589 | fhImpParResolITSselGoodTracks->Fill(pt,iparxy*10000.,clustITS); | |
590 | } | |
591 | ||
592 | return retval; | |
593 | ||
594 | } | |
595 | ||
596 | ||
597 | //--------------------------------------------------------------- | |
598 | AliAODMCParticle* AliAnalysisTaskSEHFCJqa::IsMCJet(TClonesArray *arrayMC,const AliAODJet *jet, Double_t &contribution){// assignment of parton ID to jet | |
599 | // method by L. Feldkamp | |
600 | std::vector< int > idx; | |
601 | std::vector< int > idx2; | |
602 | std::vector< double > weight; | |
603 | ||
604 | int counter =0; | |
605 | int num = jet->GetRefTracks()->GetEntries(); | |
606 | ||
607 | for(int k=0;k<num;++k){ | |
608 | ||
609 | AliAODTrack * track = (AliAODTrack*)(jet->GetRefTracks()->At(k)); | |
610 | if (track->GetLabel() >=0){ | |
611 | AliAODMCParticle* part = (AliAODMCParticle*) arrayMC->At(track->GetLabel()); | |
612 | if(!part)continue; | |
613 | ||
614 | int label =0 ; | |
615 | AliAODMCParticle* motherParton=GetMCPartonOrigin(arrayMC,part, label); | |
616 | if (!motherParton) Printf("no mother"); | |
617 | else { | |
618 | counter++; | |
619 | idx.push_back(label); //! Label of Mother | |
620 | idx2.push_back(label); | |
621 | weight.push_back(track->Pt()); //! Weight : P_t trak / P_t jet ... the latter used at the end | |
622 | } | |
623 | }///END LOOP OVER REFERENCED TRACKS | |
624 | } | |
625 | //! Remove duplicate indices for counting | |
626 | sort( idx2.begin(), idx2.end() ); | |
627 | idx2.erase( unique( idx2.begin(), idx2.end() ), idx2.end() ); | |
628 | if (idx2.size() == 0) return 0x0; | |
629 | Double_t* arrayOfWeights = new Double_t [(int)idx2.size()]; | |
630 | for(Int_t ii=0;ii<(Int_t)idx2.size();ii++)arrayOfWeights[ii]=0; | |
631 | ||
632 | for (unsigned int idxloop =0 ;idxloop<idx2.size();idxloop++){ | |
633 | for (unsigned int z=0; z< idx.size() ; ++z){ | |
634 | int a = idx.at(z); | |
635 | double w = weight.at(z); | |
636 | if(a == idx2.at(idxloop)) arrayOfWeights[idxloop] += w; | |
637 | } | |
638 | } | |
639 | ||
640 | int winner = -1; | |
641 | double c=-1.; | |
642 | for (unsigned int idxloop =0 ;idxloop<idx2.size();idxloop++){ | |
643 | if(c < arrayOfWeights[idxloop]){ | |
644 | winner =idxloop; | |
645 | c=arrayOfWeights[idxloop]; | |
646 | } | |
647 | } | |
648 | ||
649 | AliAODMCParticle *parton = 0x0; | |
650 | parton=(AliAODMCParticle*)arrayMC->At(idx.at(winner)); | |
651 | contribution = arrayOfWeights[winner]/jet->Pt(); | |
652 | ||
653 | if(arrayOfWeights) delete arrayOfWeights; | |
654 | ||
655 | ||
656 | return parton; | |
657 | ||
658 | ||
659 | } | |
660 | ||
661 | //--------------------------------------------------------------- | |
662 | AliAODMCParticle *AliAnalysisTaskSEHFCJqa::GetMCPartonOrigin(TClonesArray *arrayMC,AliAODMCParticle *p, Int_t &idx) | |
663 | { //Follows chain of track mothers until q/g or idx = -1 | |
664 | AliAODMCParticle *p2=0x0; | |
665 | Int_t mlabel = TMath::Abs(p->GetMother()) ; | |
666 | Double_t pz=0.; | |
667 | while(mlabel > 1){ | |
668 | p2 = (AliAODMCParticle*)arrayMC->At(mlabel); | |
669 | pz=TMath::Abs(p2->Pz()); | |
670 | //printf("Mother label %d, pdg %d, pz %f\n",mlabel,p2->PdgCode(),pz); | |
671 | if( p2->PdgCode() == 21 || (p2->PdgCode() != 0 && abs(p2->PdgCode()) <6) ) | |
672 | { | |
673 | idx = mlabel; | |
674 | return p2; | |
675 | } | |
676 | mlabel = TMath::Abs(p2->GetMother()); | |
677 | } | |
678 | idx=-1; | |
679 | return 0x0; | |
680 | ||
681 | } | |
682 | ||
683 | //_____________________________________________________________ | |
684 | void AliAnalysisTaskSEHFCJqa::FillJetRecoHisto(const AliAODJet *jet,Int_t partonnat,Double_t contribution,Double_t ptpart){ | |
685 | //FIll sparse with reco jets properties | |
686 | Double_t point[8]={jet->Pt(),jet->Eta(),jet->Phi()-TMath::Pi(),jet->GetRefTracks()->GetEntriesFast(),0,partonnat,contribution,ptpart}; | |
687 | fSparseRecoJets->Fill(point); | |
688 | } | |
689 | ||
690 | ||
691 | //_______________________________________________________________ | |
692 | //void AliAnalysisTaskSEHFCJqa::FillTrackHistosPID(AliAODTrack *aodtr){ | |
693 | // | |
694 | //} | |
695 | ||
696 | //_______________________________________________________________ | |
697 | void AliAnalysisTaskSEHFCJqa::Terminate(const Option_t*){ | |
698 | //TERMINATE METHOD: NOTHING TO DO | |
699 | ||
700 | ||
701 | ||
702 | } | |
703 |