new functions, mc analysis can be performed also for pi0 and k0. eta of seed can...
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskMinijet.cxx
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