]>
Commit | Line | Data |
---|---|---|
1 | #include <TChain.h> | |
2 | #include <TList.h> | |
3 | ||
4 | #include <TTree.h> | |
5 | #include <TH1F.h> | |
6 | #include <TH2F.h> | |
7 | #include <TProfile.h> | |
8 | #include <TCanvas.h> | |
9 | #include "TRandom.h" | |
10 | ||
11 | #include "AliVEvent.h" | |
12 | #include "AliVParticle.h" | |
13 | ||
14 | #include "AliESDEvent.h" | |
15 | #include "AliESDtrack.h" | |
16 | #include "AliMultiplicity.h" | |
17 | #include "AliESDtrackCuts.h" | |
18 | ||
19 | #include "AliAODEvent.h" | |
20 | #include "AliAODTrack.h" | |
21 | #include "AliAODMCParticle.h" | |
22 | ||
23 | #include "AliStack.h" | |
24 | #include "AliMCEvent.h" | |
25 | #include "AliMCParticle.h" | |
26 | #include "AliGenEventHeader.h" | |
27 | ||
28 | #include "AliAnalysisManager.h" | |
29 | ||
30 | #include "AliAnalysisTaskMinijet.h" | |
31 | ||
32 | // Analysis Task for mini jet activity analysis | |
33 | // Authors: Eva Sicking | |
34 | ||
35 | ClassImp(AliAnalysisTaskMinijet) | |
36 | ||
37 | //________________________________________________________________________ | |
38 | AliAnalysisTaskMinijet::AliAnalysisTaskMinijet(const char *name) | |
39 | : AliAnalysisTaskSE(name), | |
40 | fUseMC(kFALSE), | |
41 | fMcOnly(kFALSE), | |
42 | fCuts(0), | |
43 | fRadiusCut(0.7), | |
44 | fTriggerPtCut(0.7), | |
45 | fAssociatePtCut(0.4), | |
46 | fLeadingOrRandom(1), | |
47 | fMode(0), | |
48 | fVertexZCut(10.), | |
49 | fEtaCut(0.9), | |
50 | fEtaCutSeed(0.2), | |
51 | fSelectParticles(1), | |
52 | fESDEvent(0), | |
53 | fAODEvent(0), | |
54 | fHists(0), | |
55 | fHistPt(0), | |
56 | fHistPtMC(0), | |
57 | fChargedPi0(0) | |
58 | ||
59 | { | |
60 | for(Int_t i = 0;i< 4;i++){ | |
61 | fVertexZ[i] = 0; | |
62 | ||
63 | fPt[i] = 0; | |
64 | fEta[i] = 0; | |
65 | fPhi[i] = 0; | |
66 | ||
67 | fPtSeed[i] = 0; | |
68 | fEtaSeed[i] = 0; | |
69 | fPhiSeed[i] = 0; | |
70 | ||
71 | fPhiEta[i] = 0; | |
72 | ||
73 | fDPhiDEtaEventAxis[i] = 0; | |
74 | fTriggerNch[i] = 0; | |
75 | fTriggerTracklet[i] = 0; | |
76 | ||
77 | fNch07Nch[i] = 0; | |
78 | pNch07Nch[i] = 0; | |
79 | fNch07Tracklet[i] = 0; | |
80 | fNchTracklet[i] = 0; | |
81 | pNch07Tracklet[i] = 0; | |
82 | for(Int_t j=0;j<150;j++){ | |
83 | fDPhiEventAxisNchBin[i][j] = 0; | |
84 | fDPhiEventAxisNchBinTrig[i][j] = 0; | |
85 | ||
86 | fDPhiEventAxisTrackletBin[i][j] = 0; | |
87 | fDPhiEventAxisTrackletBinTrig[i][j] = 0; | |
88 | } | |
89 | } | |
90 | DefineOutput(1, TList::Class()); | |
91 | } | |
92 | ||
93 | //________________________________________________________________________ | |
94 | AliAnalysisTaskMinijet::~AliAnalysisTaskMinijet() | |
95 | { | |
96 | // Destructor | |
97 | ||
98 | if (fHists && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHists; | |
99 | } | |
100 | ||
101 | //________________________________________________________________________ | |
102 | void AliAnalysisTaskMinijet::UserCreateOutputObjects() | |
103 | { | |
104 | // Create histograms | |
105 | // Called once | |
106 | if(fDebug) Printf("In User Create Output Objects."); | |
107 | ||
108 | fHistPt = new TH1F("fHistPt", "P_{T} distribution REC", 15, 0.1, 3.1); | |
109 | fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); | |
110 | fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); | |
111 | fHistPt->SetMarkerStyle(kFullCircle); | |
112 | ||
113 | ||
114 | if (fUseMC) { | |
115 | fHistPtMC = new TH1F("fHistPtMC", "P_{T} distribution MC", 15, 0.1, 3.1); | |
116 | fHistPtMC->GetXaxis()->SetTitle("P_{T} (GeV/c)"); | |
117 | fHistPtMC->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); | |
118 | fHistPtMC->SetMarkerStyle(kFullCircle); | |
119 | } | |
120 | ||
121 | fChargedPi0 = new TH2F("fChargedPi0", "fChargedPi0", 200, -0.5, 199.5, 200, -0.5, 199.5); | |
122 | ||
123 | TString labels[4]={"ESDrec", "ESDmc", "AODrec", "AODmc"}; | |
124 | ||
125 | for(Int_t i=2*fMode+fMcOnly;i<1+2*fMode+fUseMC;i++){ | |
126 | ||
127 | fVertexZ[i] = new TH1F(Form("fVertexZ%s",labels[i].Data()), | |
128 | Form("fVertexZ%s",labels[i].Data()) , | |
129 | 200, -10., 10.); | |
130 | fPt[i] = new TH1F(Form("fPt%s",labels[i].Data()), | |
131 | Form("fPt%s",labels[i].Data()) , | |
132 | 500, 0., 50); | |
133 | fEta[i] = new TH1F (Form("fEta%s",labels[i].Data()), | |
134 | Form("fEta%s",labels[i].Data()) , | |
135 | 100, -1., 1); | |
136 | fPhi[i] = new TH1F(Form("fPhi%s",labels[i].Data()), | |
137 | Form("fPhi%s",labels[i].Data()) , | |
138 | 360, 0.,2*TMath::Pi()); | |
139 | ||
140 | fPtSeed[i] = new TH1F(Form("fPSeedt%s",labels[i].Data()), | |
141 | Form("fPtSeed%s",labels[i].Data()) , | |
142 | 500, 0., 50); | |
143 | fEtaSeed[i] = new TH1F (Form("fEtaSeed%s",labels[i].Data()), | |
144 | Form("fEtaSeed%s",labels[i].Data()) , | |
145 | 100, -1., 1); | |
146 | fPhiSeed[i] = new TH1F(Form("fPhiSeed%s",labels[i].Data()), | |
147 | Form("fPhiSeed%s",labels[i].Data()) , | |
148 | 360, 0.,2*TMath::Pi()); | |
149 | ||
150 | fPhiEta[i] = new TH2F(Form("fPhiEta%s",labels[i].Data()), | |
151 | Form("fPhiEta%s",labels[i].Data()) , | |
152 | 180, 0., 2*TMath::Pi(), 100, -1.,1.); | |
153 | fDPhiDEtaEventAxis[i] = new TH2F(Form("fDPhiDEtaEventAxis%s",labels[i].Data()), | |
154 | Form("fDPhiDEtaEventAxis%s",labels[i].Data()) , | |
155 | 180, -1., 2*TMath::Pi()-1, 200, -2.,2.); | |
156 | fTriggerNch[i] = new TH1F(Form("fTriggerNch%s",labels[i].Data()), | |
157 | Form("fTriggerNch%s",labels[i].Data()) , | |
158 | 250, -0.5, 249.5); | |
159 | fTriggerTracklet[i] = new TH1F(Form("fTriggerTracklet%s",labels[i].Data()), | |
160 | Form("fTriggerTracklet%s",labels[i].Data()) , | |
161 | 250, -0.5, 249.5); | |
162 | fNch07Nch[i] = new TH2F(Form("fNch07Nch%s",labels[i].Data()), | |
163 | Form("fNch07Nch%s",labels[i].Data()) , | |
164 | 250, -2.5, 247.5,250, -2.5, 247.5); | |
165 | pNch07Nch[i] = new TProfile(Form("pNch07Nch%s",labels[i].Data()), | |
166 | Form("pNch07Nch%s",labels[i].Data()) , | |
167 | 250, -2.5, 247.5); | |
168 | fNch07Tracklet[i] = new TH2F(Form("fNch07Tracklet%s",labels[i].Data()), | |
169 | Form("fNch07Tracklet%s",labels[i].Data()) , | |
170 | 250, -2.5, 247.5,250, -2.5, 247.5); | |
171 | fNchTracklet[i] = new TH2F(Form("fNchTracklet%s",labels[i].Data()), | |
172 | Form("fNchTracklet%s",labels[i].Data()) , | |
173 | 250, -2.5, 247.5,250, -2.5, 247.5); | |
174 | pNch07Tracklet[i] = new TProfile(Form("pNch07Tracklet%s",labels[i].Data()), | |
175 | Form("pNch07Tracklet%s",labels[i].Data()) , | |
176 | 250, -2.5, 247.5); | |
177 | for(Int_t j=0;j<150;j++){ | |
178 | fDPhiEventAxisNchBin[i][j] = new TH1F(Form("fDPhiEventAxisNchBin%s%02d", | |
179 | labels[i].Data(),j), | |
180 | Form("fDPhiEventAxisNchBin%s%02d", | |
181 | labels[i].Data(),j) , | |
182 | 180, 0., TMath::Pi()); | |
183 | fDPhiEventAxisNchBinTrig[i][j] = new TH1F(Form("fDPhiEventAxisNchBinTrig%s%02d", | |
184 | labels[i].Data(),j), | |
185 | Form("fDPhiEventAxisNchBinTrig%s%02d", | |
186 | labels[i].Data(),j) , | |
187 | 180, 0., TMath::Pi()); | |
188 | ||
189 | fDPhiEventAxisTrackletBin[i][j] = new TH1F(Form("fDPhiEventAxisTrackletBin%s%02d", | |
190 | labels[i].Data(),j), | |
191 | Form("fDPhiEventAxisTrackletBin%s%02d", | |
192 | labels[i].Data(),j) , | |
193 | 180, 0., TMath::Pi()); | |
194 | fDPhiEventAxisTrackletBinTrig[i][j] = new TH1F(Form("fDPhiEventAxisTrackletBinTrig%s%02d", | |
195 | labels[i].Data(),j), | |
196 | Form("fDPhiEventAxisTrackletBinTrig%s%02d", | |
197 | labels[i].Data(),j) , | |
198 | 180, 0., TMath::Pi()); | |
199 | } | |
200 | } | |
201 | ||
202 | fHists = new TList(); | |
203 | fHists->SetOwner(); | |
204 | ||
205 | fHists->Add(fHistPt); | |
206 | if(fUseMC)fHists->Add(fHistPtMC); | |
207 | fHists->Add(fChargedPi0); | |
208 | ||
209 | for(Int_t i=2*fMode+fMcOnly;i<1+2*fMode+fUseMC;i++){ | |
210 | fHists->Add(fVertexZ[i]); | |
211 | fHists->Add(fPt[i]); | |
212 | fHists->Add(fEta[i]); | |
213 | fHists->Add(fPhi[i]); | |
214 | fHists->Add(fPtSeed[i]); | |
215 | fHists->Add(fEtaSeed[i]); | |
216 | fHists->Add(fPhiSeed[i]); | |
217 | fHists->Add(fPhiEta[i]); | |
218 | fHists->Add(fDPhiDEtaEventAxis[i]); | |
219 | fHists->Add(fTriggerNch[i]); | |
220 | fHists->Add(fTriggerTracklet[i]); | |
221 | fHists->Add(fNch07Nch[i]); | |
222 | fHists->Add(pNch07Nch[i]); | |
223 | fHists->Add(fNch07Tracklet[i]); | |
224 | fHists->Add(fNchTracklet[i]); | |
225 | fHists->Add(pNch07Tracklet[i]); | |
226 | for(Int_t j=0;j<150;j++){ | |
227 | fHists->Add(fDPhiEventAxisNchBin[i][j]); | |
228 | fHists->Add(fDPhiEventAxisNchBinTrig[i][j]); | |
229 | ||
230 | fHists->Add(fDPhiEventAxisTrackletBin[i][j]); | |
231 | fHists->Add(fDPhiEventAxisTrackletBinTrig[i][j]); | |
232 | } | |
233 | } | |
234 | PostData(1, fHists); | |
235 | ||
236 | } | |
237 | ||
238 | //________________________________________________________________________ | |
239 | void AliAnalysisTaskMinijet::UserExec(Option_t *) | |
240 | { | |
241 | // Main loop | |
242 | // Called for each event | |
243 | ||
244 | AliVEvent *event = InputEvent(); | |
245 | if (!event) { | |
246 | Error("UserExec", "Could not retrieve event"); | |
247 | return; | |
248 | } | |
249 | ||
250 | //get events, either ESD or AOD | |
251 | fESDEvent = dynamic_cast<AliESDEvent*> (InputEvent()); | |
252 | fAODEvent = dynamic_cast<AliAODEvent*> (InputEvent()); | |
253 | ||
254 | ||
255 | //arrays for track properties | |
256 | Float_t *pt = 0; | |
257 | Float_t *eta = 0; | |
258 | Float_t *phi = 0; | |
259 | //number of accepted tracks and tracklets | |
260 | Int_t ntracks = 0; //return value for reading functions for ESD and AOD | |
261 | Int_t *nTracksTracklets = 0; // [0]=nAccepted,1=ntracklets,2=nall(also neutral in case of mc, | |
262 | //for real nall=ncharged) | |
263 | ||
264 | ||
265 | //read data and analyse. Implemented for ESD, ESDmc, AOD, AODmc | |
266 | //------------------------------------------------------------- | |
267 | if (fESDEvent && fMode==0) {//ESDs | |
268 | //ESD reading and analysis | |
269 | //------ | |
270 | if(!fMcOnly){ | |
271 | ntracks = LoopESD(&pt, &eta, &phi, &nTracksTracklets); //read tracks | |
272 | if(ntracks>0){ | |
273 | if(fDebug){ | |
274 | Printf("ntracks=%d", nTracksTracklets[0]); | |
275 | Printf("ntracklets=%d", nTracksTracklets[1]); | |
276 | } | |
277 | Analyse(pt, eta, phi, ntracks, nTracksTracklets[1], nTracksTracklets[2], 0); //analyse tracks | |
278 | } | |
279 | CleanArrays(pt, eta, phi, nTracksTracklets); // clean up array memory | |
280 | } | |
281 | //------ | |
282 | ||
283 | //ESD MC reading and analysis | |
284 | //------ | |
285 | if (fUseMC){ | |
286 | ntracks = LoopESDMC(&pt, &eta, &phi, &nTracksTracklets); //read mc particles | |
287 | if(ntracks>0){ | |
288 | Analyse(pt, eta, phi, ntracks, nTracksTracklets[1], nTracksTracklets[2], 1);//analyse | |
289 | } | |
290 | CleanArrays(pt, eta, phi, nTracksTracklets);// clean up array memory | |
291 | } | |
292 | //------ | |
293 | } | |
294 | ||
295 | if (fAODEvent && fMode==1) {//AODs | |
296 | ||
297 | //AOD reading and analysis | |
298 | // //------ | |
299 | if(!fMcOnly){ | |
300 | ntracks = LoopAOD(&pt, &eta, &phi, &nTracksTracklets);//read tracks | |
301 | if(ntracks>0){ | |
302 | if(fDebug){ | |
303 | Printf("ntracks=%d", nTracksTracklets[0]); | |
304 | Printf("ntracklets=%d", nTracksTracklets[1]); | |
305 | } | |
306 | Analyse(pt, eta, phi, ntracks, nTracksTracklets[1], nTracksTracklets[2], 2);//analyse | |
307 | } | |
308 | CleanArrays(pt, eta, phi, nTracksTracklets);// clean up array memory | |
309 | } | |
310 | //------ | |
311 | ||
312 | //AOD MCreading and analysis | |
313 | //------ | |
314 | if (fUseMC){ | |
315 | ntracks = LoopAODMC(&pt, &eta, &phi, &nTracksTracklets);//read tracks | |
316 | if(ntracks>0){ | |
317 | Analyse(pt, eta, phi, ntracks, nTracksTracklets[1],nTracksTracklets[2], 3);//analyse | |
318 | } | |
319 | CleanArrays(pt, eta, phi, nTracksTracklets);// clean up array memory | |
320 | } | |
321 | //------ | |
322 | } | |
323 | //------------------------------------------------------------- | |
324 | ||
325 | ||
326 | // Post output data. | |
327 | // PostData(1, fHists); | |
328 | ||
329 | } | |
330 | ||
331 | ||
332 | //________________________________________________________________________ | |
333 | Int_t AliAnalysisTaskMinijet::LoopESD(Float_t **ptArray, Float_t ** etaArray, | |
334 | Float_t ** phiArray, Int_t **nTracksTracklets ) | |
335 | { | |
336 | // gives back the number of esd tracks and pointer to arrays with track | |
337 | // properties (pt, eta, phi) | |
338 | ||
339 | // Retreive the number of all tracks for this event. | |
340 | Int_t ntracks = fESDEvent->GetNumberOfTracks(); | |
341 | if(fDebug)Printf("all ESD tracks: %d", ntracks); | |
342 | ||
343 | ||
344 | //first loop to check how many tracks are accepted | |
345 | //------------------ | |
346 | Int_t nAcceptedTracks=0; | |
347 | for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { | |
348 | AliESDtrack *track = (AliESDtrack *)fESDEvent->GetTrack(iTracks); | |
349 | if (!track) { | |
350 | Error("UserExec", "Could not receive track %d", iTracks); | |
351 | continue; | |
352 | } | |
353 | if (fCuts->AcceptTrack(track) && TMath::Abs(track->Eta())<fEtaCut) ++nAcceptedTracks; | |
354 | } | |
355 | //------------------ | |
356 | ||
357 | //generate arrays | |
358 | *ptArray = new Float_t[nAcceptedTracks]; | |
359 | *etaArray = new Float_t[nAcceptedTracks]; | |
360 | *phiArray = new Float_t[nAcceptedTracks]; | |
361 | *nTracksTracklets = new Int_t[2]; //ntracksAccepted, ntracklets | |
362 | ||
363 | //check if event is pile up or no tracks are accepted, return to user exec | |
364 | // if(fESDEvent->IsPileupFromSPD(3,0,8)) return 0; | |
365 | if(nAcceptedTracks==0) return 0; | |
366 | ||
367 | //accept event only, if vertex is good and is within fVertexZcut region | |
368 | const AliESDVertex* vertexESD = fESDEvent->GetPrimaryVertex(); | |
369 | if(vertexESD->GetNContributors()==0)return 0; | |
370 | Float_t fVz= vertexESD->GetZ(); | |
371 | if(TMath::Abs(fVz)>fVertexZCut) return 0; | |
372 | fVertexZ[0]->Fill(fVz); | |
373 | ||
374 | ||
375 | // Track loop | |
376 | Int_t iAcceptedTrack=0; | |
377 | for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { | |
378 | AliESDtrack *track = (AliESDtrack *)fESDEvent->GetTrack(iTracks); | |
379 | if (!track) { | |
380 | Error("UserExec", "Could not receive track %d", iTracks); | |
381 | continue; | |
382 | } | |
383 | if (fCuts->AcceptTrack(track) && TMath::Abs(track->Eta())<fEtaCut){ | |
384 | (*ptArray)[iAcceptedTrack] = track->Pt(); | |
385 | (*etaArray)[iAcceptedTrack] = track->Eta(); | |
386 | (*phiArray)[iAcceptedTrack++] = track->Phi(); | |
387 | fHistPt->Fill(track->Pt()); | |
388 | } | |
389 | } | |
390 | ||
391 | ||
392 | //tracklet loop | |
393 | Int_t ntrackletsAccept=0; | |
394 | AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity()); | |
395 | Int_t ntracklets = mult->GetNumberOfTracklets(); | |
396 | for(Int_t i=0;i< ntracklets;i++){ | |
397 | if(mult->GetDeltaPhi(i)<0.05){ | |
398 | ntrackletsAccept++; | |
399 | } | |
400 | } | |
401 | ||
402 | (*nTracksTracklets)[0] = nAcceptedTracks; | |
403 | (*nTracksTracklets)[1] = ntrackletsAccept; | |
404 | (*nTracksTracklets)[2] = nAcceptedTracks;//in order to be compatible with mc analysis | |
405 | //where here also neutral particles are counted. | |
406 | ||
407 | return nAcceptedTracks; | |
408 | ||
409 | } | |
410 | ||
411 | //________________________________________________________________________ | |
412 | Int_t AliAnalysisTaskMinijet::LoopESDMC(Float_t **ptArray, Float_t ** etaArray, | |
413 | Float_t ** phiArray, Int_t ** nTracksTracklets) | |
414 | { | |
415 | // gives back the number of charged prim MC particle and pointer to arrays | |
416 | // with particle properties (pt, eta, phi) | |
417 | ||
418 | // Printf("Event starts: Loop ESD MC"); | |
419 | ||
420 | AliMCEvent *mcEvent = (AliMCEvent*) MCEvent(); | |
421 | if (!mcEvent) { | |
422 | Error("UserExec", "Could not retrieve MC event"); | |
423 | return 0; | |
424 | } | |
425 | ||
426 | AliStack* stack = 0x0; | |
427 | if(fUseMC) stack = MCEvent()->Stack(); | |
428 | ||
429 | Int_t ntracks = mcEvent->GetNumberOfTracks(); | |
430 | if(fDebug)Printf("MC particles: %d", ntracks); | |
431 | ||
432 | ||
433 | //---------------------------------- | |
434 | //first track loop to check how many chared primary tracks are available | |
435 | Int_t nChargedPrimaries=0; | |
436 | Int_t nAllPrimaries=0; | |
437 | ||
438 | Int_t nPseudoTracklets=0; | |
439 | for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { | |
440 | AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks)); | |
441 | if (!track) { | |
442 | Error("UserExec", "Could not receive track %d", iTracks); | |
443 | continue; | |
444 | } | |
445 | ||
446 | ||
447 | if(//count also charged particles in case of fSelectParticles==2 (only neutral) | |
448 | !SelectParticlePlusCharged( | |
449 | track->Charge(), | |
450 | track->PdgCode(), | |
451 | stack->IsPhysicalPrimary(track->Label()) | |
452 | ) | |
453 | ) | |
454 | continue; | |
455 | ||
456 | ||
457 | // Printf("fSelectParticles= %d\n", fSelectParticles); | |
458 | //count number of pseudo tracklets | |
459 | if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.035) ++nPseudoTracklets; | |
460 | //same cuts as on ESDtracks | |
461 | if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 || track->Pt()>200) continue; | |
462 | ||
463 | //count all primaries | |
464 | ++nAllPrimaries; | |
465 | //count charged primaries | |
466 | if (track->Charge() != 0) ++nChargedPrimaries; | |
467 | ||
468 | // Printf("PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) ); | |
469 | ||
470 | ||
471 | } | |
472 | //---------------------------------- | |
473 | ||
474 | // Printf("All in acceptance=%d",nAllPrimaries); | |
475 | // Printf("Charged in acceptance =%d",nChargedPrimaries); | |
476 | ||
477 | fChargedPi0->Fill(nAllPrimaries,nChargedPrimaries); | |
478 | ||
479 | if(fSelectParticles!=2){ | |
480 | *ptArray = new Float_t[nAllPrimaries]; | |
481 | *etaArray = new Float_t[nAllPrimaries]; | |
482 | *phiArray = new Float_t[nAllPrimaries]; | |
483 | } | |
484 | else{ | |
485 | *ptArray = new Float_t[nAllPrimaries-nChargedPrimaries]; | |
486 | *etaArray = new Float_t[nAllPrimaries-nChargedPrimaries]; | |
487 | *phiArray = new Float_t[nAllPrimaries-nChargedPrimaries]; | |
488 | } | |
489 | ||
490 | *nTracksTracklets = new Int_t[3]; | |
491 | ||
492 | if(nAllPrimaries==0) return 0; | |
493 | if(nChargedPrimaries==0) return 0; | |
494 | ||
495 | ||
496 | //vertex cut | |
497 | AliGenEventHeader* header = MCEvent()->GenEventHeader(); | |
498 | TArrayF mcV; | |
499 | header->PrimaryVertex(mcV); | |
500 | Float_t vzMC = mcV[2]; | |
501 | if(TMath::Abs(vzMC)>fVertexZCut) return 0; | |
502 | fVertexZ[1]->Fill(vzMC); | |
503 | ||
504 | ||
505 | //track loop | |
506 | Int_t iChargedPiK=0; | |
507 | for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { | |
508 | AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks)); | |
509 | if (!track) { | |
510 | Error("UserExec", "Could not receive track %d", iTracks); | |
511 | continue; | |
512 | } | |
513 | ||
514 | if(!SelectParticle( | |
515 | track->Charge(), | |
516 | track->PdgCode(), | |
517 | stack->IsPhysicalPrimary(track->Label()) | |
518 | ) | |
519 | ) | |
520 | continue; | |
521 | ||
522 | ||
523 | //same cuts as on ESDtracks | |
524 | if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 || track->Pt()>200) continue; | |
525 | ||
526 | // Printf("After: PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) ); | |
527 | ||
528 | ||
529 | fHistPtMC->Fill(track->Pt()); | |
530 | //fills arrays with track properties | |
531 | (*ptArray)[iChargedPiK] = track->Pt(); | |
532 | (*etaArray)[iChargedPiK] = track->Eta(); | |
533 | (*phiArray)[iChargedPiK++] = track->Phi(); | |
534 | ||
535 | } //track loop | |
536 | ||
537 | (*nTracksTracklets)[0] = nChargedPrimaries; | |
538 | (*nTracksTracklets)[1] = nPseudoTracklets; | |
539 | if(fSelectParticles!=2){ | |
540 | (*nTracksTracklets)[2] = nAllPrimaries; | |
541 | } | |
542 | else{ | |
543 | (*nTracksTracklets)[2] = nAllPrimaries-nChargedPrimaries; // only neutral | |
544 | } | |
545 | return nChargedPrimaries; | |
546 | ||
547 | } | |
548 | ||
549 | //________________________________________________________________________ | |
550 | Int_t AliAnalysisTaskMinijet::LoopAOD(Float_t **ptArray, Float_t ** etaArray, | |
551 | Float_t ** phiArray, Int_t ** nTracksTracklets) | |
552 | { | |
553 | // gives back the number of AOD tracks and pointer to arrays with track | |
554 | // properties (pt, eta, phi) | |
555 | ||
556 | ||
557 | // Retreive the number of tracks for this event. | |
558 | Int_t ntracks = fAODEvent->GetNumberOfTracks(); | |
559 | if(fDebug) Printf("AOD tracks: %d", ntracks); | |
560 | ||
561 | ||
562 | Int_t nAcceptedTracks=0; | |
563 | for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { | |
564 | AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks); | |
565 | if (!track) { | |
566 | Error("UserExec", "Could not receive track %d", iTracks); | |
567 | continue; | |
568 | } | |
569 | if(track->TestFilterBit(16) && TMath::Abs(track->Eta())<fEtaCut) nAcceptedTracks++; | |
570 | } | |
571 | ||
572 | *ptArray = new Float_t[nAcceptedTracks]; | |
573 | *etaArray = new Float_t[nAcceptedTracks]; | |
574 | *phiArray = new Float_t[nAcceptedTracks]; | |
575 | *nTracksTracklets = new Int_t[3]; //here, third entry only copy of first (compatibility for MC) | |
576 | ||
577 | ||
578 | if(nAcceptedTracks==0) return 0; | |
579 | AliAODVertex* vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertex(); | |
580 | ||
581 | // TODO: need to check how to implement IsPileupFromSPD(3,0.8) | |
582 | // function of esd event | |
583 | // first solution: exclude this call in esd loop for comparison (QA) | |
584 | ||
585 | if(!vertex) return 0; | |
586 | Double_t vzAOD=vertex->GetZ(); | |
587 | //if(vertex->GetNContributors()==0) return 0; | |
588 | if(TMath::Abs(vzAOD)<1e-9) return 0; | |
589 | ||
590 | if(TMath::Abs(vzAOD)>fVertexZCut) return 0; | |
591 | fVertexZ[2]->Fill(vzAOD); | |
592 | ||
593 | // Track loop to fill a pT spectrum | |
594 | Int_t iAcceptedTracks=0; | |
595 | for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) { | |
596 | AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks); | |
597 | if (!track) { | |
598 | Error("UserExec", "Could not receive track %d", iTracks); | |
599 | continue; | |
600 | } | |
601 | if(!track->TestFilterBit(16) || TMath::Abs(track->Eta())>fEtaCut) continue; | |
602 | fHistPt->Fill(track->Pt()); | |
603 | ||
604 | //fills arrays with track properties | |
605 | (*ptArray)[iAcceptedTracks] = track->Pt(); | |
606 | (*etaArray)[iAcceptedTracks] = track->Eta(); | |
607 | (*phiArray)[iAcceptedTracks++] = track->Phi(); | |
608 | ||
609 | } //track loop | |
610 | ||
611 | //tracklet loop | |
612 | Int_t ntrackletsAccept=0; | |
613 | AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets(); | |
614 | for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){ | |
615 | if(TMath::Abs(mult->GetDeltaPhi(i))<0.05){ | |
616 | ++ntrackletsAccept; | |
617 | } | |
618 | } | |
619 | ||
620 | (*nTracksTracklets)[0] = nAcceptedTracks; | |
621 | (*nTracksTracklets)[1] = ntrackletsAccept; | |
622 | (*nTracksTracklets)[2] = nAcceptedTracks; | |
623 | ||
624 | return nAcceptedTracks; | |
625 | ||
626 | } | |
627 | ||
628 | //________________________________________________________________________ | |
629 | Int_t AliAnalysisTaskMinijet::LoopAODMC(Float_t **ptArray, Float_t ** etaArray, | |
630 | Float_t ** phiArray, Int_t ** nTracksTracklets) | |
631 | { | |
632 | // gives back the number of AOD MC particles and pointer to arrays with particle | |
633 | // properties (pt, eta, phi) | |
634 | ||
635 | //retreive MC particles from event | |
636 | TClonesArray *mcArray = (TClonesArray*)fAODEvent-> | |
637 | FindListObject(AliAODMCParticle::StdBranchName()); | |
638 | if(!mcArray){ | |
639 | Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__); | |
640 | return kFALSE; | |
641 | } | |
642 | ||
643 | Int_t ntracks = mcArray->GetEntriesFast(); | |
644 | if(fDebug)Printf("MC particles: %d", ntracks); | |
645 | ||
646 | ||
647 | // Track loop: chek how many particles will be accepted | |
648 | Float_t vzMC=0.; | |
649 | Int_t nChargedPrim=0; | |
650 | Int_t nAllPrim=0; | |
651 | Int_t nPseudoTracklets=0; | |
652 | for (Int_t it = 0; it < ntracks; it++) { | |
653 | AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it); | |
654 | if (!track) { | |
655 | Error("UserExec", "Could not receive track %d", it); | |
656 | continue; | |
657 | } | |
658 | ||
659 | if(//count also charged particles in case of fSelectParticles==2 (only neutral) | |
660 | !SelectParticlePlusCharged( | |
661 | track->Charge(), | |
662 | track->PdgCode(), | |
663 | track->IsPhysicalPrimary() | |
664 | ) | |
665 | ) | |
666 | continue; | |
667 | ||
668 | ||
669 | if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.035)++nPseudoTracklets; | |
670 | if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 || track->Pt()>200) continue; | |
671 | // if(TMath::Abs(track->Eta())<1e-9 && TMath::Abs(track->Phi())<1e-9)continue; | |
672 | ||
673 | //same cuts as in ESD filter | |
674 | if(!track->TestBit(16))continue; //cuts set in ESD filter during AOD generation | |
675 | ||
676 | nAllPrim++; | |
677 | if(track->Charge()!=0) nChargedPrim++; | |
678 | Printf("eta=%f,phi=%f,pt=%f",track->Eta(),track->Phi(),track->Pt()); | |
679 | Printf("prim=%d,pdg=%d,charge=%d",track->IsPhysicalPrimary(),track->PdgCode(),track->Charge()); | |
680 | ||
681 | ||
682 | if(nAllPrim==1) vzMC= track->Zv(); // check only one time. (only one vertex per event allowed) | |
683 | } | |
684 | ||
685 | //generate array with size of number of accepted tracks | |
686 | fChargedPi0->Fill(nAllPrim,nChargedPrim); | |
687 | ||
688 | if(fSelectParticles!=2){ | |
689 | *ptArray = new Float_t[nAllPrim]; | |
690 | *etaArray = new Float_t[nAllPrim]; | |
691 | *phiArray = new Float_t[nAllPrim]; | |
692 | } | |
693 | else{ | |
694 | *ptArray = new Float_t[nAllPrim-nChargedPrim]; | |
695 | *etaArray = new Float_t[nAllPrim-nChargedPrim]; | |
696 | *phiArray = new Float_t[nAllPrim-nChargedPrim]; | |
697 | } | |
698 | ||
699 | ||
700 | *nTracksTracklets = new Int_t[3]; | |
701 | ||
702 | // Printf("nAllPrim=%d", nAllPrim); | |
703 | // Printf("nChargedPrim=%d", nChargedPrim); | |
704 | ||
705 | if(nAllPrim==0) return 0; | |
706 | if(nChargedPrim==0) return 0; | |
707 | ||
708 | ||
709 | if(TMath::Abs(vzMC)>fVertexZCut) return 0; | |
710 | fVertexZ[3]->Fill(vzMC); | |
711 | ||
712 | ||
713 | // Track loop: fill arrays for accepted tracks | |
714 | Int_t iChargedPrim=0; | |
715 | for (Int_t it = 0; it < ntracks; it++) { | |
716 | AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it); | |
717 | if (!track) { | |
718 | Error("UserExec", "Could not receive track %d", it); | |
719 | continue; | |
720 | } | |
721 | ||
722 | if(!SelectParticle( | |
723 | track->Charge(), | |
724 | track->PdgCode(), | |
725 | track->IsPhysicalPrimary() | |
726 | ) | |
727 | ) | |
728 | continue; | |
729 | ||
730 | ||
731 | if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 || track->Pt()>200) continue; | |
732 | //if(TMath::Abs(track->Eta())<1e-8 && TMath::Abs(track->Phi())<1e-8)continue; | |
733 | ||
734 | //if(track->TestBit(16))continue; | |
735 | ||
736 | fHistPtMC->Fill(track->Pt()); | |
737 | (*ptArray)[iChargedPrim] = track->Pt(); | |
738 | (*etaArray)[iChargedPrim] = track->Eta(); | |
739 | (*phiArray)[iChargedPrim++] = track->Phi(); | |
740 | ||
741 | } | |
742 | ||
743 | (*nTracksTracklets)[0] = nChargedPrim; | |
744 | (*nTracksTracklets)[1] = nPseudoTracklets; | |
745 | (*nTracksTracklets)[2] = nAllPrim; | |
746 | ||
747 | return nChargedPrim; | |
748 | // Printf("ntracks=%d", nChargedPrim); | |
749 | ||
750 | } | |
751 | ||
752 | //________________________________________________________________________ | |
753 | void AliAnalysisTaskMinijet::Analyse(Float_t *pt, Float_t *eta, Float_t *phi, Int_t ntracksCharged, | |
754 | Int_t ntracklets, Int_t nAll, Int_t mode) | |
755 | { | |
756 | ||
757 | // analyse track properties (comming from either ESDs or AODs) in order to compute | |
758 | // mini jet activity (chared tracks) as function of charged multiplicity | |
759 | ||
760 | // ntracks and ntracklets are already the number of accepted tracks and tracklets | |
761 | ||
762 | if(fDebug) Printf("In Analysis\n"); | |
763 | ||
764 | Float_t ptEventAxis=0; // pt leading | |
765 | Float_t etaEventAxis=0; // eta leading | |
766 | Float_t phiEventAxis=0; // phi leading | |
767 | ||
768 | Float_t ptOthers = 0; // pt others // for all other tracks around event axis -> see loop | |
769 | Float_t etaOthers = 0; // eta others | |
770 | Float_t phiOthers = 0; // phi others | |
771 | ||
772 | Int_t *pindexInnerEta = new Int_t[nAll]; | |
773 | Float_t *ptInnerEta = new Float_t[nAll]; | |
774 | ||
775 | for (Int_t i = 0; i < nAll; i++) { | |
776 | //filling of simple check plots | |
777 | fPt[mode] -> Fill( pt[i]); | |
778 | fEta[mode] -> Fill(eta[i]); | |
779 | fPhi[mode] -> Fill(phi[i]); | |
780 | fPhiEta[mode]-> Fill(phi[i], eta[i]); | |
781 | ||
782 | pindexInnerEta[i]=0; //set all values to zero | |
783 | //fill new array for eta check | |
784 | ptInnerEta[i]= pt[i]; | |
785 | ||
786 | } | |
787 | ||
788 | ||
789 | ||
790 | // define event axis: leading or random track (with pt>fTriggerPtCut) | |
791 | // --------------------------------------- | |
792 | Int_t highPtTracks=0; | |
793 | Int_t highPtTracksInnerEta=0; | |
794 | ||
795 | ||
796 | //count high pt tracks and high pt tracks in acceptance for seeds | |
797 | for(Int_t i=0;i<nAll;i++){ | |
798 | ||
799 | if(pt[i]>fTriggerPtCut) { | |
800 | highPtTracks++; | |
801 | } | |
802 | ||
803 | // seed should be place in middle of acceptance, that complete cone is in acceptance | |
804 | if(pt[i]>fTriggerPtCut && TMath::Abs(eta[i])<fEtaCutSeed){ | |
805 | ||
806 | // Printf("eta=%f", eta[i]); | |
807 | highPtTracksInnerEta++; | |
808 | } | |
809 | else{ | |
810 | ptInnerEta[i]=0; | |
811 | } | |
812 | } | |
813 | ||
814 | ||
815 | //sort array in order to get highest pt tracks first | |
816 | //index can be computed with pindexInnerEta[number] | |
817 | if(nAll) TMath::Sort(nAll, ptInnerEta, pindexInnerEta, kTRUE); | |
818 | ||
819 | ||
820 | // plot of multiplicity distributions | |
821 | fNch07Nch[mode]->Fill(ntracksCharged, highPtTracks); | |
822 | pNch07Nch[mode]->Fill(ntracksCharged, highPtTracks); | |
823 | if(ntracklets){ | |
824 | fNch07Tracklet[mode]->Fill(ntracklets, highPtTracks); | |
825 | fNchTracklet[mode]->Fill(ntracklets, ntracksCharged); | |
826 | pNch07Tracklet[mode]->Fill(ntracklets, highPtTracks); | |
827 | } | |
828 | ||
829 | //analysis can only be performed with event axis, defined by high pt track | |
830 | ||
831 | ||
832 | if(highPtTracks>0 && highPtTracksInnerEta>0){ | |
833 | ||
834 | //Printf("%s:%d",(char*)__FILE__,__LINE__); | |
835 | //check setter of event axis | |
836 | //default option: random=1, | |
837 | //second option:leading=0 | |
838 | Int_t axis=-1; | |
839 | if(fLeadingOrRandom==0)axis=0; | |
840 | else if (fLeadingOrRandom==1)axis= Int_t((highPtTracksInnerEta)*gRandom->Rndm()); | |
841 | else Printf("Wrong settings for event axis."); | |
842 | ||
843 | if(fDebug){ | |
844 | Printf("Axis tracks has pT=%f, test=%f", ptInnerEta[pindexInnerEta[axis]], pt[pindexInnerEta[axis]]); | |
845 | Printf("Axis tracks has eta=%f", eta[pindexInnerEta[axis]]); | |
846 | } | |
847 | ||
848 | //--------------------------------------- | |
849 | ||
850 | ||
851 | if(ntracksCharged>1){ // require at least two tracks (leading and prob. accosicates) | |
852 | ||
853 | //EventAxisRandom track properties | |
854 | ptEventAxis = pt [pindexInnerEta[axis]]; | |
855 | etaEventAxis = eta[pindexInnerEta[axis]]; | |
856 | phiEventAxis = phi[pindexInnerEta[axis]]; | |
857 | fPtSeed[mode] -> Fill( ptEventAxis); | |
858 | fEtaSeed[mode] -> Fill(etaEventAxis); | |
859 | fPhiSeed[mode] -> Fill(phiEventAxis); | |
860 | ||
861 | //track loop for event propoerties around event axis with pt>triggerPtCut | |
862 | //loop only over already accepted tracks except event axis | |
863 | if(ptEventAxis>fTriggerPtCut){ | |
864 | ||
865 | for (Int_t iTrack = 0; iTrack < nAll; iTrack++) { | |
866 | ||
867 | if(pindexInnerEta[axis]==iTrack)continue; // no double counting | |
868 | ||
869 | ptOthers = pt [iTrack]; | |
870 | etaOthers = eta[iTrack]; | |
871 | phiOthers = phi[iTrack]; | |
872 | ||
873 | ||
874 | if(ptOthers>fAssociatePtCut){ // only tracks which fullfill associate pt cut | |
875 | ||
876 | Float_t dPhi=TMath::Abs(phiOthers-phiEventAxis); | |
877 | if(dPhi>TMath::Pi()) dPhi=2*TMath::Pi()-dPhi; | |
878 | Float_t dEta=etaOthers-etaEventAxis; | |
879 | ||
880 | Float_t dphiplot = phiOthers-phiEventAxis; | |
881 | if(dphiplot>2*TMath::Pi()-1) dphiplot = dphiplot-2*TMath::Pi(); | |
882 | else if(dphiplot<-1)dphiplot=dphiplot+2*TMath::Pi(); | |
883 | fDPhiDEtaEventAxis[mode]->Fill(dphiplot, dEta); | |
884 | ||
885 | if(ntracksCharged<150){ | |
886 | fDPhiEventAxisNchBin[mode][ntracksCharged]->Fill(dPhi); | |
887 | if(ptOthers>fTriggerPtCut) | |
888 | fDPhiEventAxisNchBinTrig[mode][ntracksCharged]->Fill(dPhi); | |
889 | } | |
890 | ||
891 | if(ntracklets<150){ | |
892 | fDPhiEventAxisTrackletBin[mode][ntracklets]->Fill(dPhi); | |
893 | if(ptOthers>fTriggerPtCut) | |
894 | fDPhiEventAxisTrackletBinTrig[mode][ntracklets]->Fill(dPhi); | |
895 | } | |
896 | ||
897 | }//tracks fulfill assoc track cut | |
898 | ||
899 | }// end track loop | |
900 | ||
901 | ||
902 | // fill histogram with number of tracks (pt>fAssociatePtCut) around event axis | |
903 | // how often is there a trigger particle at a certain Nch bin | |
904 | fTriggerNch[mode]->Fill(ntracksCharged); | |
905 | fTriggerTracklet[mode]->Fill(ntracklets); | |
906 | ||
907 | }//if track pt is at least trigger pt | |
908 | ||
909 | } //if there are more than 1 track | |
910 | ||
911 | }//if there is at least one high pt track | |
912 | ||
913 | ||
914 | if(pindexInnerEta){// clean up array memory used for TMath::Sort | |
915 | delete[] pindexInnerEta; | |
916 | pindexInnerEta=0; | |
917 | } | |
918 | ||
919 | if(ptInnerEta){// clean up array memory used for TMath::Sort | |
920 | delete[] ptInnerEta; | |
921 | ptInnerEta=0; | |
922 | } | |
923 | ||
924 | PostData(1, fHists); | |
925 | ||
926 | } | |
927 | ||
928 | //________________________________________________________________________ | |
929 | void AliAnalysisTaskMinijet::Terminate(Option_t*) | |
930 | { | |
931 | //terminate function is called at the end | |
932 | //can be used to draw histograms etc. | |
933 | ||
934 | } | |
935 | ||
936 | //________________________________________________________________________ | |
937 | Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(Int_t charge, Int_t pdg, Bool_t prim) | |
938 | { | |
939 | //selection of mc particle | |
940 | //fSelectParticles=0: use charged primaries and pi0 and k0 | |
941 | //fSelectParticles=1: use only charged primaries | |
942 | //fSelectParticles=2: use only pi0 and k0 | |
943 | ||
944 | if(fSelectParticles==0 || fSelectParticles==2){ // in case of 2: need to count also charged particles | |
945 | if(charge==0){ | |
946 | if(!(pdg==111||pdg==130||pdg==310)) | |
947 | return false; | |
948 | } | |
949 | else{// charge !=0 | |
950 | if(!prim) | |
951 | return false; | |
952 | } | |
953 | } | |
954 | ||
955 | else if(fSelectParticles==1){ | |
956 | if (charge==0 || !prim){ | |
957 | return false; | |
958 | } | |
959 | } | |
960 | ||
961 | else{ | |
962 | Printf("Error: wrong selection of charged/pi0/k0"); | |
963 | return 0; | |
964 | } | |
965 | ||
966 | return true; | |
967 | ||
968 | } | |
969 | ||
970 | //________________________________________________________________________ | |
971 | Bool_t AliAnalysisTaskMinijet::SelectParticle(Int_t charge, Int_t pdg, Bool_t prim) | |
972 | { | |
973 | //selection of mc particle | |
974 | //fSelectParticles=0: use charged primaries and pi0 and k0 | |
975 | //fSelectParticles=1: use only charged primaries | |
976 | //fSelectParticles=2: use only pi0 and k0 | |
977 | ||
978 | if(fSelectParticles==0){ | |
979 | if(charge==0){ | |
980 | if(!(pdg==111||pdg==130||pdg==310)) | |
981 | return false; | |
982 | } | |
983 | else{// charge !=0 | |
984 | if(!prim) | |
985 | return false; | |
986 | } | |
987 | } | |
988 | ||
989 | else if (fSelectParticles==1){ | |
990 | if (charge==0 || !prim){ | |
991 | return false; | |
992 | } | |
993 | } | |
994 | else if(fSelectParticles==2){ | |
995 | if(!(pdg==111||pdg==130||pdg==310)) | |
996 | return false; | |
997 | } | |
998 | ||
999 | return true; | |
1000 | ||
1001 | } | |
1002 | ||
1003 | //________________________________________________________________________ | |
1004 | void AliAnalysisTaskMinijet::CleanArrays(Float_t* pt, Float_t* eta, Float_t* phi,Int_t* nTracksTracklets) | |
1005 | { | |
1006 | //clean up of memory used for arrays of track properties | |
1007 | ||
1008 | if(pt){ | |
1009 | delete[] pt; | |
1010 | pt=0; | |
1011 | } | |
1012 | if(eta){ | |
1013 | delete[] eta; | |
1014 | eta=0; | |
1015 | } | |
1016 | if(phi){ | |
1017 | delete[] phi; | |
1018 | phi=0; | |
1019 | } | |
1020 | if(nTracksTracklets){ | |
1021 | delete[] nTracksTracklets; | |
1022 | nTracksTracklets=0; | |
1023 | } | |
1024 | ||
1025 | } | |
1026 |