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