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