new functions, mc analysis can be performed also for pi0 and k0. eta of seed can...
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskMinijet.cxx
CommitLineData
117f99e3 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"
b7c0438e 21#include "AliAODMCParticle.h"
117f99e3 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
35ClassImp(AliAnalysisTaskMinijet)
36
37//________________________________________________________________________
b9ad6f04 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
117f99e3 59{
60 for(Int_t i = 0;i< 4;i++){
b9ad6f04 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;
117f99e3 76
b9ad6f04 77 fNch07Nch[i] = 0;
78 pNch07Nch[i] = 0;
79 fNch07Tracklet[i] = 0;
80 fNchTracklet[i] = 0;
81 pNch07Tracklet[i] = 0;
117f99e3 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//________________________________________________________________________
94AliAnalysisTaskMinijet::~AliAnalysisTaskMinijet()
95{
96 // Destructor
97
98 if (fHists && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHists;
99}
100
101//________________________________________________________________________
102void AliAnalysisTaskMinijet::UserCreateOutputObjects()
103{
104 // Create histograms
105 // Called once
b9ad6f04 106 if(fDebug) Printf("In User Create Output Objects.");
117f99e3 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
b9ad6f04 121 fChargedPi0 = new TH2F("fChargedPi0", "fChargedPi0", 200, -0.5, 199.5, 200, -0.5, 199.5);
122
117f99e3 123 TString labels[4]={"ESDrec", "ESDmc", "AODrec", "AODmc"};
124
b9ad6f04 125 for(Int_t i=2*fMode+fMcOnly;i<1+2*fMode+fUseMC;i++){
117f99e3 126
127 fVertexZ[i] = new TH1F(Form("fVertexZ%s",labels[i].Data()),
b9ad6f04 128 Form("fVertexZ%s",labels[i].Data()) ,
129 200, -10., 10.);
117f99e3 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()),
b9ad6f04 134 Form("fEta%s",labels[i].Data()) ,
135 100, -1., 1);
117f99e3 136 fPhi[i] = new TH1F(Form("fPhi%s",labels[i].Data()),
b9ad6f04 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.);
117f99e3 153 fDPhiDEtaEventAxis[i] = new TH2F(Form("fDPhiDEtaEventAxis%s",labels[i].Data()),
b9ad6f04 154 Form("fDPhiDEtaEventAxis%s",labels[i].Data()) ,
155 180, -1., 2*TMath::Pi()-1, 200, -2.,2.);
117f99e3 156 fTriggerNch[i] = new TH1F(Form("fTriggerNch%s",labels[i].Data()),
b9ad6f04 157 Form("fTriggerNch%s",labels[i].Data()) ,
158 250, -0.5, 249.5);
117f99e3 159 fTriggerTracklet[i] = new TH1F(Form("fTriggerTracklet%s",labels[i].Data()),
b9ad6f04 160 Form("fTriggerTracklet%s",labels[i].Data()) ,
161 250, -0.5, 249.5);
117f99e3 162 fNch07Nch[i] = new TH2F(Form("fNch07Nch%s",labels[i].Data()),
b9ad6f04 163 Form("fNch07Nch%s",labels[i].Data()) ,
164 250, -2.5, 247.5,250, -2.5, 247.5);
117f99e3 165 pNch07Nch[i] = new TProfile(Form("pNch07Nch%s",labels[i].Data()),
b9ad6f04 166 Form("pNch07Nch%s",labels[i].Data()) ,
167 250, -2.5, 247.5);
117f99e3 168 fNch07Tracklet[i] = new TH2F(Form("fNch07Tracklet%s",labels[i].Data()),
b9ad6f04 169 Form("fNch07Tracklet%s",labels[i].Data()) ,
170 250, -2.5, 247.5,250, -2.5, 247.5);
117f99e3 171 fNchTracklet[i] = new TH2F(Form("fNchTracklet%s",labels[i].Data()),
b9ad6f04 172 Form("fNchTracklet%s",labels[i].Data()) ,
173 250, -2.5, 247.5,250, -2.5, 247.5);
117f99e3 174 pNch07Tracklet[i] = new TProfile(Form("pNch07Tracklet%s",labels[i].Data()),
b9ad6f04 175 Form("pNch07Tracklet%s",labels[i].Data()) ,
176 250, -2.5, 247.5);
117f99e3 177 for(Int_t j=0;j<150;j++){
178 fDPhiEventAxisNchBin[i][j] = new TH1F(Form("fDPhiEventAxisNchBin%s%02d",
b9ad6f04 179 labels[i].Data(),j),
180 Form("fDPhiEventAxisNchBin%s%02d",
181 labels[i].Data(),j) ,
182 180, 0., TMath::Pi());
117f99e3 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",
b9ad6f04 190 labels[i].Data(),j),
191 Form("fDPhiEventAxisTrackletBin%s%02d",
192 labels[i].Data(),j) ,
193 180, 0., TMath::Pi());
117f99e3 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);
b9ad6f04 207 fHists->Add(fChargedPi0);
117f99e3 208
b9ad6f04 209 for(Int_t i=2*fMode+fMcOnly;i<1+2*fMode+fUseMC;i++){
117f99e3 210 fHists->Add(fVertexZ[i]);
211 fHists->Add(fPt[i]);
212 fHists->Add(fEta[i]);
213 fHists->Add(fPhi[i]);
b9ad6f04 214 fHists->Add(fPtSeed[i]);
215 fHists->Add(fEtaSeed[i]);
216 fHists->Add(fPhiSeed[i]);
217 fHists->Add(fPhiEta[i]);
117f99e3 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//________________________________________________________________________
239void AliAnalysisTaskMinijet::UserExec(Option_t *)
240{
241 // Main loop
242 // Called for each event
243
244 AliVEvent *event = InputEvent();
245 if (!event) {
b9ad6f04 246 Error("UserExec", "Could not retrieve event");
247 return;
117f99e3 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
b9ad6f04 261 Int_t *nTracksTracklets = 0; // [0]=nAccepted,1=ntracklets,2=nall(also neutral in case of mc,
262 //for real nall=ncharged)
117f99e3 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 //------
b9ad6f04 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
117f99e3 278 }
b9ad6f04 279 CleanArrays(pt, eta, phi, nTracksTracklets); // clean up array memory
117f99e3 280 }
117f99e3 281 //------
282
283 //ESD MC reading and analysis
284 //------
285 if (fUseMC){
b9ad6f04 286 ntracks = LoopESDMC(&pt, &eta, &phi, &nTracksTracklets); //read mc particles
117f99e3 287 if(ntracks>0){
b9ad6f04 288 Analyse(pt, eta, phi, ntracks, nTracksTracklets[1], nTracksTracklets[2], 1);//analyse
117f99e3 289 }
b9ad6f04 290 CleanArrays(pt, eta, phi, nTracksTracklets);// clean up array memory
117f99e3 291 }
292 //------
293 }
294
295 if (fAODEvent && fMode==1) {//AODs
296
297 //AOD reading and analysis
b9ad6f04 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
117f99e3 307 }
b9ad6f04 308 CleanArrays(pt, eta, phi, nTracksTracklets);// clean up array memory
117f99e3 309 }
b9ad6f04 310 //------
311
117f99e3 312 //AOD MCreading and analysis
313 //------
314 if (fUseMC){
b9ad6f04 315 ntracks = LoopAODMC(&pt, &eta, &phi, &nTracksTracklets);//read tracks
117f99e3 316 if(ntracks>0){
b9ad6f04 317 Analyse(pt, eta, phi, ntracks, nTracksTracklets[1],nTracksTracklets[2], 3);//analyse
117f99e3 318 }
b9ad6f04 319 CleanArrays(pt, eta, phi, nTracksTracklets);// clean up array memory
117f99e3 320 }
321 //------
322 }
323 //-------------------------------------------------------------
324
325
326 // Post output data.
327 // PostData(1, fHists);
328
329}
330
331
332//________________________________________________________________________
333Int_t AliAnalysisTaskMinijet::LoopESD(Float_t **ptArray, Float_t ** etaArray,
b9ad6f04 334 Float_t ** phiArray, Int_t **nTracksTracklets )
117f99e3 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 }
b9ad6f04 353 if (fCuts->AcceptTrack(track) && TMath::Abs(track->Eta())<fEtaCut) ++nAcceptedTracks;
117f99e3 354 }
355 //------------------
356
357 //generate arrays
358 *ptArray = new Float_t[nAcceptedTracks];
359 *etaArray = new Float_t[nAcceptedTracks];
360 *phiArray = new Float_t[nAcceptedTracks];
b9ad6f04 361 *nTracksTracklets = new Int_t[2]; //ntracksAccepted, ntracklets
117f99e3 362
363 //check if event is pile up or no tracks are accepted, return to user exec
b9ad6f04 364 // if(fESDEvent->IsPileupFromSPD(3,0,8)) return 0;
117f99e3 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 }
b9ad6f04 383 if (fCuts->AcceptTrack(track) && TMath::Abs(track->Eta())<fEtaCut){
117f99e3 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
b9ad6f04 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.
117f99e3 406
407 return nAcceptedTracks;
408
409}
410
411//________________________________________________________________________
412Int_t AliAnalysisTaskMinijet::LoopESDMC(Float_t **ptArray, Float_t ** etaArray,
b9ad6f04 413 Float_t ** phiArray, Int_t ** nTracksTracklets)
117f99e3 414{
415 // gives back the number of charged prim MC particle and pointer to arrays
416 // with particle properties (pt, eta, phi)
417
b9ad6f04 418 // Printf("Event starts: Loop ESD MC");
419
117f99e3 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;
b9ad6f04 436 Int_t nAllPrimaries=0;
437
438 Int_t nPseudoTracklets=0;
117f99e3 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 }
b9ad6f04 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;
117f99e3 460 //same cuts as on ESDtracks
b9ad6f04 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
117f99e3 471 }
472 //----------------------------------
473
b9ad6f04 474 // Printf("All in acceptance=%d",nAllPrimaries);
475 // Printf("Charged in acceptance =%d",nChargedPrimaries);
476
477 fChargedPi0->Fill(nAllPrimaries,nChargedPrimaries);
117f99e3 478
b9ad6f04 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];
117f99e3 491
b9ad6f04 492 if(nAllPrimaries==0) return 0;
117f99e3 493 if(nChargedPrimaries==0) return 0;
494
b9ad6f04 495
496 //vertex cut
117f99e3 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);
117f99e3 503
504
505 //track loop
b9ad6f04 506 Int_t iChargedPiK=0;
117f99e3 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 }
b9ad6f04 513
514 if(!SelectParticle(
515 track->Charge(),
516 track->PdgCode(),
517 stack->IsPhysicalPrimary(track->Label())
518 )
519 )
520 continue;
521
522
117f99e3 523 //same cuts as on ESDtracks
b9ad6f04 524 if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<0.2 || track->Pt()>200) continue;
117f99e3 525
b9ad6f04 526 // Printf("After: PDG=%d, IsPrim=%d", track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
527
117f99e3 528
529 fHistPtMC->Fill(track->Pt());
530 //fills arrays with track properties
b9ad6f04 531 (*ptArray)[iChargedPiK] = track->Pt();
532 (*etaArray)[iChargedPiK] = track->Eta();
533 (*phiArray)[iChargedPiK++] = track->Phi();
117f99e3 534
535 } //track loop
536
b9ad6f04 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 }
117f99e3 545 return nChargedPrimaries;
546
547}
548
549//________________________________________________________________________
550Int_t AliAnalysisTaskMinijet::LoopAOD(Float_t **ptArray, Float_t ** etaArray,
b9ad6f04 551 Float_t ** phiArray, Int_t ** nTracksTracklets)
117f99e3 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 }
b9ad6f04 569 if(track->TestFilterBit(16) && TMath::Abs(track->Eta())<fEtaCut) nAcceptedTracks++;
117f99e3 570 }
571
572 *ptArray = new Float_t[nAcceptedTracks];
573 *etaArray = new Float_t[nAcceptedTracks];
574 *phiArray = new Float_t[nAcceptedTracks];
b9ad6f04 575 *nTracksTracklets = new Int_t[3]; //here, third entry only copy of first (compatibility for MC)
117f99e3 576
577
578 if(nAcceptedTracks==0) return 0;
b9ad6f04 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
117f99e3 585 if(!vertex) return 0;
586 Double_t vzAOD=vertex->GetZ();
b9ad6f04 587 //if(vertex->GetNContributors()==0) return 0;
588 if(TMath::Abs(vzAOD)<1e-9) return 0;
589
117f99e3 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) {
b9ad6f04 598 Error("UserExec", "Could not receive track %d", iTracks);
599 continue;
117f99e3 600 }
b9ad6f04 601 if(!track->TestFilterBit(16) || TMath::Abs(track->Eta())>fEtaCut) continue;
117f99e3 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
b9ad6f04 620 (*nTracksTracklets)[0] = nAcceptedTracks;
621 (*nTracksTracklets)[1] = ntrackletsAccept;
622 (*nTracksTracklets)[2] = nAcceptedTracks;
117f99e3 623
624 return nAcceptedTracks;
625
626}
627
628//________________________________________________________________________
629Int_t AliAnalysisTaskMinijet::LoopAODMC(Float_t **ptArray, Float_t ** etaArray,
b9ad6f04 630 Float_t ** phiArray, Int_t ** nTracksTracklets)
117f99e3 631{
b7c0438e 632 // gives back the number of AOD MC particles and pointer to arrays with particle
633 // properties (pt, eta, phi)
117f99e3 634
b7c0438e 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;
117f99e3 641 }
642
b7c0438e 643 Int_t ntracks = mcArray->GetEntriesFast();
117f99e3 644 if(fDebug)Printf("MC particles: %d", ntracks);
645
b7c0438e 646
647 // Track loop: chek how many particles will be accepted
648 Float_t vzMC=0.;
b9ad6f04 649 Int_t nChargedPrim=0;
650 Int_t nAllPrim=0;
651 Int_t nPseudoTracklets=0;
b7c0438e 652 for (Int_t it = 0; it < ntracks; it++) {
653 AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
117f99e3 654 if (!track) {
b7c0438e 655 Error("UserExec", "Could not receive track %d", it);
117f99e3 656 continue;
657 }
b9ad6f04 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)
b7c0438e 683 }
117f99e3 684
b7c0438e 685 //generate array with size of number of accepted tracks
b9ad6f04 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];
b7c0438e 701
b9ad6f04 702 // Printf("nAllPrim=%d", nAllPrim);
703 // Printf("nChargedPrim=%d", nChargedPrim);
704
705 if(nAllPrim==0) return 0;
706 if(nChargedPrim==0) return 0;
117f99e3 707
b9ad6f04 708
b7c0438e 709 if(TMath::Abs(vzMC)>fVertexZCut) return 0;
710 fVertexZ[3]->Fill(vzMC);
711
117f99e3 712
b7c0438e 713 // Track loop: fill arrays for accepted tracks
b9ad6f04 714 Int_t iChargedPrim=0;
b7c0438e 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 }
b9ad6f04 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
b7c0438e 741 }
b9ad6f04 742
743 (*nTracksTracklets)[0] = nChargedPrim;
744 (*nTracksTracklets)[1] = nPseudoTracklets;
745 (*nTracksTracklets)[2] = nAllPrim;
117f99e3 746
b9ad6f04 747 return nChargedPrim;
748 // Printf("ntracks=%d", nChargedPrim);
117f99e3 749
750}
751
752//________________________________________________________________________
b9ad6f04 753void AliAnalysisTaskMinijet::Analyse(Float_t *pt, Float_t *eta, Float_t *phi, Int_t ntracksCharged,
754 Int_t ntracklets, Int_t nAll, Int_t mode)
117f99e3 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
b9ad6f04 772 Int_t *pindexInnerEta = new Int_t[nAll];
773 Float_t *ptInnerEta = new Float_t[nAll];
117f99e3 774
b9ad6f04 775 for (Int_t i = 0; i < nAll; i++) {
117f99e3 776 //filling of simple check plots
b9ad6f04 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
117f99e3 786 }
787
117f99e3 788
789
790 // define event axis: leading or random track (with pt>fTriggerPtCut)
791 // ---------------------------------------
792 Int_t highPtTracks=0;
b9ad6f04 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 }
117f99e3 812 }
813
b9ad6f04 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
117f99e3 820 // plot of multiplicity distributions
b9ad6f04 821 fNch07Nch[mode]->Fill(ntracksCharged, highPtTracks);
822 pNch07Nch[mode]->Fill(ntracksCharged, highPtTracks);
117f99e3 823 if(ntracklets){
824 fNch07Tracklet[mode]->Fill(ntracklets, highPtTracks);
b9ad6f04 825 fNchTracklet[mode]->Fill(ntracklets, ntracksCharged);
117f99e3 826 pNch07Tracklet[mode]->Fill(ntracklets, highPtTracks);
827 }
828
829 //analysis can only be performed with event axis, defined by high pt track
117f99e3 830
b9ad6f04 831
832 if(highPtTracks>0 && highPtTracksInnerEta>0){
833
834 //Printf("%s:%d",(char*)__FILE__,__LINE__);
117f99e3 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;
b9ad6f04 840 else if (fLeadingOrRandom==1)axis= Int_t((highPtTracksInnerEta)*gRandom->Rndm());
117f99e3 841 else Printf("Wrong settings for event axis.");
b9ad6f04 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 }
117f99e3 847
848 //---------------------------------------
849
850
b9ad6f04 851 if(ntracksCharged>1){ // require at least two tracks (leading and prob. accosicates)
117f99e3 852
853 //EventAxisRandom track properties
b9ad6f04 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);
117f99e3 860
117f99e3 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
b9ad6f04 865 for (Int_t iTrack = 0; iTrack < nAll; iTrack++) {
117f99e3 866
b9ad6f04 867 if(pindexInnerEta[axis]==iTrack)continue; // no double counting
117f99e3 868
b9ad6f04 869 ptOthers = pt [iTrack];
870 etaOthers = eta[iTrack];
871 phiOthers = phi[iTrack];
117f99e3 872
b9ad6f04 873
117f99e3 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;
b9ad6f04 879
117f99e3 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
b9ad6f04 885 if(ntracksCharged<150){
886 fDPhiEventAxisNchBin[mode][ntracksCharged]->Fill(dPhi);
117f99e3 887 if(ptOthers>fTriggerPtCut)
b9ad6f04 888 fDPhiEventAxisNchBinTrig[mode][ntracksCharged]->Fill(dPhi);
117f99e3 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
b9ad6f04 904 fTriggerNch[mode]->Fill(ntracksCharged);
117f99e3 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
b9ad6f04 913
914 if(pindexInnerEta){// clean up array memory used for TMath::Sort
915 delete[] pindexInnerEta;
916 pindexInnerEta=0;
117f99e3 917 }
918
b9ad6f04 919 if(ptInnerEta){// clean up array memory used for TMath::Sort
920 delete[] ptInnerEta;
921 ptInnerEta=0;
922 }
923
117f99e3 924 PostData(1, fHists);
925
926}
927
928//________________________________________________________________________
929void 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//________________________________________________________________________
b9ad6f04 937Bool_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//________________________________________________________________________
971Bool_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//________________________________________________________________________
1004void AliAnalysisTaskMinijet::CleanArrays(Float_t* pt, Float_t* eta, Float_t* phi,Int_t* nTracksTracklets)
117f99e3 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 }
b9ad6f04 1020 if(nTracksTracklets){
1021 delete[] nTracksTracklets;
1022 nTracksTracklets=0;
117f99e3 1023 }
1024
1025}
b9ad6f04 1026