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     fCuts(0),
  42     fRadiusCut(0.7),
  43     fTriggerPtCut(0.7),
  44     fAssociatePtCut(0.4),
  45     fLeadingOrRandom(1),
  46     fMode(0),
  47     fVertexZCut(10.),
  48     fESDEvent(0),
  49     fAODEvent(0),
  50     fHists(0),
  51     fHistPt(0),
  52     fHistPtMC(0)
  53
  54 {
  55   for(Int_t i = 0;i< 4;i++){
  56     fVertexZ[i]               = 0;
  57     fPt[i]                    = 0;
  58     fEta[i]                   = 0;
  59     fPhi[i]                   = 0;
  60
  61     fDPhiDEtaEventAxis[i]     = 0;
  62     fTriggerNch[i]            = 0;
  63     fTriggerTracklet[i]       = 0;
  64
  65     fNch07Nch[i]              = 0;
  66     pNch07Nch[i]              = 0;
  67     fNch07Tracklet[i]         = 0;
  68     fNchTracklet[i]         = 0;
  69     pNch07Tracklet[i]         = 0;
  70     for(Int_t j=0;j<150;j++){
  71       fDPhiEventAxisNchBin[i][j]   = 0;
  72       fDPhiEventAxisNchBinTrig[i][j]   = 0;
  73
  74       fDPhiEventAxisTrackletBin[i][j]   = 0;
  75       fDPhiEventAxisTrackletBinTrig[i][j]   = 0;
  76     }
  77   }
  78   DefineOutput(1, TList::Class());
  79 }
  80
  81 //________________________________________________________________________
  82 AliAnalysisTaskMinijet::~AliAnalysisTaskMinijet()
  83 {
  84   // Destructor
  85
  86   if (fHists && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHists;
  87 }
  88
  89 //________________________________________________________________________
  90 void AliAnalysisTaskMinijet::UserCreateOutputObjects()
  91 {
  92   // Create histograms
  93   // Called once
  94
  95
  96   fHistPt = new TH1F("fHistPt", "P_{T} distribution REC", 15, 0.1, 3.1);
  97   fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
  98   fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
  99   fHistPt->SetMarkerStyle(kFullCircle);
 100
 101
 102   if (fUseMC) {
 103     fHistPtMC = new TH1F("fHistPtMC", "P_{T} distribution MC", 15, 0.1, 3.1);
 104     fHistPtMC->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 105     fHistPtMC->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 106     fHistPtMC->SetMarkerStyle(kFullCircle);
 107   }
 108
 109   TString labels[4]={"ESDrec", "ESDmc", "AODrec", "AODmc"};
 110
 111   for(Int_t i=2*fMode;i<1+2*fMode+fUseMC;i++){
 112
 113     fVertexZ[i]                  = new TH1F(Form("fVertexZ%s",labels[i].Data()),
 114                                                 Form("fVertexZ%s",labels[i].Data()) ,
 115                                             200, -10., 10);
 116     fPt[i]                       = new TH1F(Form("fPt%s",labels[i].Data()),
 117                                             Form("fPt%s",labels[i].Data()) ,
 118                                             500, 0., 50);
 119     fEta[i]                      = new TH1F (Form("fEta%s",labels[i].Data()),
 120                                                  Form("fEta%s",labels[i].Data()) ,
 121                                                  100, -1., 1);
 122     fPhi[i]                      = new TH1F(Form("fPhi%s",labels[i].Data()),
 123                                                 Form("fPhi%s",labels[i].Data()) ,
 124                                                 360, 0.,2*TMath::Pi());
 125     fDPhiDEtaEventAxis[i]        = new TH2F(Form("fDPhiDEtaEventAxis%s",labels[i].Data()),
 126                                                 Form("fDPhiDEtaEventAxis%s",labels[i].Data()) ,
 127                                                 180, -1., 2*TMath::Pi()-1, 200, -2.,2.);
 128     fTriggerNch[i]               = new TH1F(Form("fTriggerNch%s",labels[i].Data()),
 129                                                 Form("fTriggerNch%s",labels[i].Data()) ,
 130                                                 250, -0.5, 249.5);
 131     fTriggerTracklet[i]          = new TH1F(Form("fTriggerTracklet%s",labels[i].Data()),
 132                                                 Form("fTriggerTracklet%s",labels[i].Data()) ,
 133                                                 250, -0.5, 249.5);
 134     fNch07Nch[i]                 = new TH2F(Form("fNch07Nch%s",labels[i].Data()),
 135                                                 Form("fNch07Nch%s",labels[i].Data()) ,
 136                                                 250, -2.5, 247.5,250, -2.5, 247.5);
 137     pNch07Nch[i]                 = new TProfile(Form("pNch07Nch%s",labels[i].Data()),
 138                                                     Form("pNch07Nch%s",labels[i].Data()) ,
 139                                                     250, -2.5, 247.5);
 140     fNch07Tracklet[i]            = new TH2F(Form("fNch07Tracklet%s",labels[i].Data()),
 141                                                 Form("fNch07Tracklet%s",labels[i].Data()) ,
 142                                                 250, -2.5, 247.5,250, -2.5, 247.5);
 143     fNchTracklet[i]              = new TH2F(Form("fNchTracklet%s",labels[i].Data()),
 144                                                 Form("fNchTracklet%s",labels[i].Data()) ,
 145                                                 250, -2.5, 247.5,250, -2.5, 247.5);
 146     pNch07Tracklet[i]            = new TProfile(Form("pNch07Tracklet%s",labels[i].Data()),
 147                                                     Form("pNch07Tracklet%s",labels[i].Data()) ,
 148                                                     250, -2.5, 247.5);
 149     for(Int_t j=0;j<150;j++){
 150       fDPhiEventAxisNchBin[i][j]          = new TH1F(Form("fDPhiEventAxisNchBin%s%02d",
 151                                                               labels[i].Data(),j),
 152                                                          Form("fDPhiEventAxisNchBin%s%02d",
 153                                                               labels[i].Data(),j) ,
 154                                                          180, 0., TMath::Pi());
 155       fDPhiEventAxisNchBinTrig[i][j]          = new TH1F(Form("fDPhiEventAxisNchBinTrig%s%02d",
 156                                                               labels[i].Data(),j),
 157                                                          Form("fDPhiEventAxisNchBinTrig%s%02d",
 158                                                               labels[i].Data(),j) ,
 159                                                          180, 0., TMath::Pi());
 160
 161       fDPhiEventAxisTrackletBin[i][j]     = new TH1F(Form("fDPhiEventAxisTrackletBin%s%02d",
 162                                                               labels[i].Data(),j),
 163                                                          Form("fDPhiEventAxisTrackletBin%s%02d",
 164                                                               labels[i].Data(),j) ,
 165                                                          180, 0., TMath::Pi());
 166       fDPhiEventAxisTrackletBinTrig[i][j]     = new TH1F(Form("fDPhiEventAxisTrackletBinTrig%s%02d",
 167                                                               labels[i].Data(),j),
 168                                                          Form("fDPhiEventAxisTrackletBinTrig%s%02d",
 169                                                               labels[i].Data(),j) ,
 170                                                          180, 0., TMath::Pi());
 171     }
 172   }
 173
 174   fHists = new TList();
 175   fHists->SetOwner();
 176
 177   fHists->Add(fHistPt);
 178   if(fUseMC)fHists->Add(fHistPtMC);
 179
 180   for(Int_t i=2*fMode;i<1+2*fMode+fUseMC;i++){
 181     fHists->Add(fVertexZ[i]);
 182     fHists->Add(fPt[i]);
 183     fHists->Add(fEta[i]);
 184     fHists->Add(fPhi[i]);
 185     fHists->Add(fDPhiDEtaEventAxis[i]);
 186     fHists->Add(fTriggerNch[i]);
 187     fHists->Add(fTriggerTracklet[i]);
 188     fHists->Add(fNch07Nch[i]);
 189     fHists->Add(pNch07Nch[i]);
 190     fHists->Add(fNch07Tracklet[i]);
 191     fHists->Add(fNchTracklet[i]);
 192     fHists->Add(pNch07Tracklet[i]);
 193     for(Int_t j=0;j<150;j++){
 194       fHists->Add(fDPhiEventAxisNchBin[i][j]);
 195       fHists->Add(fDPhiEventAxisNchBinTrig[i][j]);
 196
 197       fHists->Add(fDPhiEventAxisTrackletBin[i][j]);
 198       fHists->Add(fDPhiEventAxisTrackletBinTrig[i][j]);
 199     }
 200   }
 201   PostData(1, fHists);
 202
 203 }
 204
 205 //________________________________________________________________________
 206 void AliAnalysisTaskMinijet::UserExec(Option_t *)
 207 {
 208   // Main loop
 209   // Called for each event
 210
 211   AliVEvent *event = InputEvent();
 212   if (!event) {
 213      Error("UserExec", "Could not retrieve event");
 214      return;
 215   }
 216
 217   //get events, either ESD or AOD
 218   fESDEvent = dynamic_cast<AliESDEvent*> (InputEvent());
 219   fAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
 220
 221
 222   //arrays for track properties
 223   Float_t *pt  = 0;
 224   Float_t *eta = 0;
 225   Float_t *phi = 0;
 226   //number of accepted tracks and tracklets
 227   Int_t ntracks = 0;  //return value for reading functions for ESD and AOD
 228   Int_t *numbers = 0; // numbers[0]=nAcceptedTracks, numbers[1]=ntracklets,
 229
 230
 231   //read data and analyse. Implemented for ESD, ESDmc, AOD, AODmc
 232   //-------------------------------------------------------------
 233   if (fESDEvent && fMode==0) {//ESDs
 234     //ESD reading and analysis
 235     //------
 236     ntracks = LoopESD(&pt, &eta, &phi, &numbers); //read tracks
 237     if(ntracks>0){
 238       if(fDebug){
 239         Printf("ntracks=%d", numbers[0]);
 240         Printf("ntracklets=%d", numbers[1]);
 241       }
 242       Analyse(pt, eta, phi, ntracks, numbers[1], 0); //analyse tracks
 243     }
 244     CleanArrays(pt, eta, phi, numbers); // clean up array memory
 245     //------
 246
 247     //ESD MC reading and analysis
 248     //------
 249     if (fUseMC){
 250       ntracks = LoopESDMC(&pt, &eta, &phi); //read mc particles
 251       if(ntracks>0){
 252         Analyse(pt, eta, phi, ntracks, 0, 1);//analyse
 253       }
 254       CleanArrays(pt, eta, phi);// clean up array memory
 255     }
 256     //------
 257   }
 258
 259   if (fAODEvent && fMode==1) {//AODs
 260
 261     //AOD reading and analysis
 262     //------
 263     ntracks = LoopAOD(&pt, &eta, &phi, &numbers);//read tracks
 264     if(ntracks>0){
 265       if(fDebug){
 266         Printf("ntracks=%d", numbers[0]);
 267         Printf("ntracklets=%d", numbers[1]);
 268       }
 269       Analyse(pt, eta, phi, ntracks, numbers[1], 2);//analyse
 270     }
 271     CleanArrays(pt, eta, phi, numbers);// clean up array memory
 272     //------
 273
 274     //AOD MCreading and analysis
 275     //------
 276     if (fUseMC){
 277       ntracks = LoopAODMC(&pt, &eta, &phi);//read tracks
 278       if(ntracks>0){
 279         Analyse(pt, eta, phi, ntracks, 0, 3);//analyse
 280       }
 281       CleanArrays(pt, eta, phi);// clean up array memory
 282     }
 283     //------
 284   }
 285   //-------------------------------------------------------------
 286
 287
 288   // Post output data.
 289   //  PostData(1, fHists);
 290
 291 }
 292
 293
 294 //________________________________________________________________________
 295 Int_t AliAnalysisTaskMinijet::LoopESD(Float_t **ptArray, Float_t ** etaArray,
 296                                       Float_t ** phiArray, Int_t **numbers )
 297 {
 298   // gives back the number of esd tracks and pointer to arrays with track
 299   // properties (pt, eta, phi)
 300
 301   // Retreive the number of all tracks for this event.
 302   Int_t ntracks = fESDEvent->GetNumberOfTracks();
 303   if(fDebug)Printf("all ESD tracks: %d", ntracks);
 304
 305
 306   //first loop to check how many tracks are accepted
 307   //------------------
 308   Int_t nAcceptedTracks=0;
 309   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 310     AliESDtrack *track = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
 311     if (!track) {
 312       Error("UserExec", "Could not receive track %d", iTracks);
 313       continue;
 314     }
 315     if (fCuts->AcceptTrack(track)) ++nAcceptedTracks;
 316   }
 317   //------------------
 318
 319   //generate arrays
 320   *ptArray = new Float_t[nAcceptedTracks];
 321   *etaArray = new Float_t[nAcceptedTracks];
 322   *phiArray = new Float_t[nAcceptedTracks];
 323   *numbers = new Int_t[2]; //ntracksAccepted, ntracklets
 324
 325   //check if event is pile up or no tracks are accepted, return to user exec
 326   if(fESDEvent->IsPileupFromSPD(3,0,8)) return 0;
 327   if(nAcceptedTracks==0) return 0;
 328
 329   //accept event only, if vertex is good and is within fVertexZcut region
 330   const AliESDVertex*   vertexESD   = fESDEvent->GetPrimaryVertex();
 331   if(vertexESD->GetNContributors()==0)return 0;
 332   Float_t fVz= vertexESD->GetZ();
 333   if(TMath::Abs(fVz)>fVertexZCut) return 0;
 334   fVertexZ[0]->Fill(fVz);
 335
 336
 337   // Track loop
 338   Int_t iAcceptedTrack=0;
 339   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 340     AliESDtrack *track = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
 341     if (!track) {
 342       Error("UserExec", "Could not receive track %d", iTracks);
 343       continue;
 344     }
 345     if (fCuts->AcceptTrack(track)){
 346       (*ptArray)[iAcceptedTrack]  = track->Pt();
 347       (*etaArray)[iAcceptedTrack] = track->Eta();
 348       (*phiArray)[iAcceptedTrack++] = track->Phi();
 349       fHistPt->Fill(track->Pt());
 350     }
 351   }
 352
 353
 354   //tracklet loop
 355   Int_t ntrackletsAccept=0;
 356   AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
 357   Int_t ntracklets = mult->GetNumberOfTracklets();
 358   for(Int_t i=0;i< ntracklets;i++){
 359     if(mult->GetDeltaPhi(i)<0.05){
 360       ntrackletsAccept++;
 361     }
 362   }
 363
 364   (*numbers)[0] = nAcceptedTracks;
 365   (*numbers)[1] = ntrackletsAccept;
 366
 367   return nAcceptedTracks;
 368
 369 }
 370
 371 //________________________________________________________________________
 372 Int_t AliAnalysisTaskMinijet::LoopESDMC(Float_t **ptArray, Float_t ** etaArray,
 373                                         Float_t ** phiArray)
 374 {
 375   // gives back the number of charged prim MC particle and pointer to arrays
 376   // with particle properties (pt, eta, phi)
 377
 378   AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
 379   if (!mcEvent) {
 380     Error("UserExec", "Could not retrieve MC event");
 381     return 0;
 382   }
 383
 384   AliStack* stack = 0x0;
 385   if(fUseMC) stack = MCEvent()->Stack();
 386
 387   Int_t ntracks = mcEvent->GetNumberOfTracks();
 388   if(fDebug)Printf("MC particles: %d", ntracks);
 389
 390
 391   //----------------------------------
 392   //first track loop to check how many chared primary tracks are available
 393   Int_t nChargedPrimaries=0;
 394   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 395     AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
 396     if (!track) {
 397       Error("UserExec", "Could not receive track %d", iTracks);
 398       continue;
 399     }
 400     if (!(stack->IsPhysicalPrimary(track->Label()))) continue;
 401     if (track->Particle()->GetPDG()->Charge() == 0) continue;
 402     //same cuts as on ESDtracks
 403     if(TMath::Abs(track->Eta())>0.9 || track->Pt()<0.2 || track->Pt()>200) continue;
 404     ++nChargedPrimaries;
 405   }
 406   //----------------------------------
 407
 408
 409   *ptArray = new Float_t[nChargedPrimaries];
 410   *etaArray = new Float_t[nChargedPrimaries];
 411   *phiArray = new Float_t[nChargedPrimaries];
 412
 413   if(nChargedPrimaries==0) return 0;
 414
 415   AliGenEventHeader*  header = MCEvent()->GenEventHeader();
 416   TArrayF mcV;
 417   header->PrimaryVertex(mcV);
 418   Float_t vzMC = mcV[2];
 419   if(TMath::Abs(vzMC)>fVertexZCut) return 0;
 420   fVertexZ[1]->Fill(vzMC);
 421
 422
 423
 424   //track loop
 425   Int_t iChargedPrimaries=0;
 426   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 427     AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
 428     if (!track) {
 429       Error("UserExec", "Could not receive track %d", iTracks);
 430       continue;
 431     }
 432     if (!(stack->IsPhysicalPrimary(track->Label()))) continue;
 433     if (track->Particle()->GetPDG()->Charge() == 0) continue;
 434     //same cuts as on ESDtracks
 435     if(TMath::Abs(track->Eta())>0.9 || track->Pt()<0.2 || track->Pt()>200) continue;
 436
 437
 438     fHistPtMC->Fill(track->Pt());
 439     //fills arrays with track properties
 440     (*ptArray)[iChargedPrimaries]  = track->Pt();
 441     (*etaArray)[iChargedPrimaries] = track->Eta();
 442     (*phiArray)[iChargedPrimaries++] = track->Phi();
 443
 444   } //track loop
 445
 446   return nChargedPrimaries;
 447
 448 }
 449
 450 //________________________________________________________________________
 451 Int_t AliAnalysisTaskMinijet::LoopAOD(Float_t **ptArray, Float_t ** etaArray,
 452                                       Float_t ** phiArray, Int_t ** numbers)
 453 {
 454   // gives back the number of AOD tracks and pointer to arrays with track
 455   // properties (pt, eta, phi)
 456
 457
 458   // Retreive the number of tracks for this event.
 459   Int_t ntracks = fAODEvent->GetNumberOfTracks();
 460   if(fDebug) Printf("AOD tracks: %d", ntracks);
 461
 462
 463   Int_t nAcceptedTracks=0;
 464   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 465     AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
 466     if (!track) {
 467       Error("UserExec", "Could not receive track %d", iTracks);
 468       continue;
 469     }
 470     if(track->TestFilterBit(16))nAcceptedTracks++;
 471   }
 472
 473   *ptArray = new Float_t[nAcceptedTracks];
 474   *etaArray = new Float_t[nAcceptedTracks];
 475   *phiArray = new Float_t[nAcceptedTracks];
 476   *numbers = new Int_t[2];
 477
 478
 479   if(nAcceptedTracks==0) return 0;
 480   AliAODVertex* vertex= fAODEvent->GetPrimaryVertex();
 481   if(!vertex) return 0;
 482   Double_t vzAOD=vertex->GetZ();
 483   if(vertex->GetNContributors()==0) return 0;
 484   if(TMath::Abs(vzAOD)>fVertexZCut) return 0;
 485   fVertexZ[2]->Fill(vzAOD);
 486
 487   // Track loop to fill a pT spectrum
 488   Int_t iAcceptedTracks=0;
 489   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 490     AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
 491     if (!track) {
 492      Error("UserExec", "Could not receive track %d", iTracks);
 493      continue;
 494     }
 495     if(!track->TestFilterBit(16))continue;
 496     fHistPt->Fill(track->Pt());
 497
 498     //fills arrays with track properties
 499     (*ptArray)[iAcceptedTracks]  = track->Pt();
 500     (*etaArray)[iAcceptedTracks] = track->Eta();
 501     (*phiArray)[iAcceptedTracks++] = track->Phi();
 502
 503   } //track loop
 504
 505   //tracklet loop
 506   Int_t ntrackletsAccept=0;
 507   AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
 508   for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
 509     if(TMath::Abs(mult->GetDeltaPhi(i))<0.05){
 510       ++ntrackletsAccept;
 511     }
 512   }
 513
 514   (*numbers)[0] = nAcceptedTracks;
 515   (*numbers)[1] = ntrackletsAccept;
 516
 517   return nAcceptedTracks;
 518
 519 }
 520
 521 //________________________________________________________________________
 522 Int_t AliAnalysisTaskMinijet::LoopAODMC(Float_t **ptArray, Float_t ** etaArray,
 523                                         Float_t ** phiArray)
 524 {
 525   // gives back the number of AOD MC particles and pointer to arrays with particle
 526   // properties (pt, eta, phi)
 527
 528   //retreive MC particles from event
 529   TClonesArray *mcArray = (TClonesArray*)fAODEvent->
 530     FindListObject(AliAODMCParticle::StdBranchName());
 531   if(!mcArray){
 532     Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__);
 533     return kFALSE;
 534   }
 535
 536   Int_t ntracks = mcArray->GetEntriesFast();
 537   if(fDebug)Printf("MC particles: %d", ntracks);
 538
 539
 540   // Track loop: chek how many particles will be accepted
 541   Float_t vzMC=0.;
 542   Int_t nAcceptedTracks=0;
 543   for (Int_t it = 0; it < ntracks; it++) {
 544     AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
 545     if (!track) {
 546       Error("UserExec", "Could not receive track %d", it);
 547       continue;
 548     }
 549     if(!track->IsPhysicalPrimary())continue;
 550     if (track->Charge() == 0) continue;
 551     if(TMath::Abs(track->Eta())>1. || track->Pt()<0.2 || track->Pt()>200) continue; //same cuts as in ESD filter
 552     if(track->TestBit(16))nAcceptedTracks++;
 553     if(nAcceptedTracks==1) vzMC= track->Zv(); // check only one time. (only one vertex per event allowed)
 554   }
 555
 556   //generate array with size of number of accepted tracks
 557   *ptArray = new Float_t[nAcceptedTracks];
 558   *etaArray = new Float_t[nAcceptedTracks];
 559   *phiArray = new Float_t[nAcceptedTracks];
 560
 561   if(nAcceptedTracks==0) return 0;
 562
 563   //check vertex
 564   if(TMath::Abs(vzMC)>fVertexZCut) return 0;
 565   fVertexZ[3]->Fill(vzMC);
 566
 567
 568   // Track loop: fill arrays for accepted tracks
 569   Int_t iAcceptedTracks=0;
 570   for (Int_t it = 0; it < ntracks; it++) {
 571     AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
 572     if (!track) {
 573       Error("UserExec", "Could not receive track %d", it);
 574       continue;
 575     }
 576     if(!track->IsPhysicalPrimary())continue;
 577     if (track->Charge() == 0) continue;
 578     if(TMath::Abs(track->Eta())>0.9 || track->Pt()<0.2 || track->Pt()>200) continue;
 579
 580     if(track->TestBit(16)){
 581       fHistPtMC->Fill(track->Pt());
 582       (*ptArray)[iAcceptedTracks]  = track->Pt();
 583       (*etaArray)[iAcceptedTracks] = track->Eta();
 584       (*phiArray)[iAcceptedTracks++] = track->Phi();
 585     }
 586   }
 587
 588   return nAcceptedTracks;
 589
 590 }
 591
 592 //________________________________________________________________________
 593 void AliAnalysisTaskMinijet::Analyse(Float_t *pt, Float_t *eta, Float_t *phi, Int_t ntracks,
 594                                      Int_t ntracklets, Int_t mode)
 595 {
 596
 597   // analyse track properties (comming from either ESDs or AODs) in order to compute
 598   // mini jet activity (chared tracks) as function of charged multiplicity
 599
 600   // ntracks and ntracklets are already the number of accepted tracks and tracklets
 601
 602   if(fDebug) Printf("In Analysis\n");
 603
 604   Float_t ptEventAxis=0;  // pt leading
 605   Float_t etaEventAxis=0; // eta leading
 606   Float_t phiEventAxis=0; // phi leading
 607
 608   Float_t ptOthers  = 0; // pt others // for all other tracks around event axis -> see loop
 609   Float_t etaOthers = 0; // eta others
 610   Float_t phiOthers = 0; // phi others
 611
 612   Int_t *pindex  = new Int_t[ntracks];//index needed for sorting of pt array
 613
 614   for (Int_t i = 0; i < ntracks; i++) {
 615     //filling of simple check plots
 616     fPt[mode]  ->Fill( pt[i]);
 617     fEta[mode] ->Fill(eta[i]);
 618     fPhi[mode] ->Fill(phi[i]);
 619     pindex[i]=0; //set all values to zero
 620   }
 621
 622   //sort all tracks according to their pt
 623   if(ntracks) TMath::Sort(ntracks, pt, pindex, kTRUE);
 624
 625
 626   // define event axis: leading or random track (with pt>fTriggerPtCut)
 627   // ---------------------------------------
 628   Int_t highPtTracks=0;
 629   for(Int_t i=0;i<ntracks;i++){//show just the filled entries, skip empty ones.
 630     if(fDebug) Printf("%d:  pt = %f, number %i \n",mode, pt[pindex[i]],i);
 631     if(pt[pindex[i]]>fTriggerPtCut) highPtTracks++;
 632   }
 633
 634   //  plot of multiplicity distributions
 635   fNch07Nch[mode]->Fill(ntracks, highPtTracks);
 636   pNch07Nch[mode]->Fill(ntracks, highPtTracks);
 637   if(ntracklets){
 638     fNch07Tracklet[mode]->Fill(ntracklets, highPtTracks);
 639     fNchTracklet[mode]->Fill(ntracklets, ntracks);
 640     pNch07Tracklet[mode]->Fill(ntracklets, highPtTracks);
 641   }
 642
 643   //analysis can only be performed with event axis, defined by high pt track
 644   if(highPtTracks>0){
 645
 646     //check setter of event axis
 647     //default option: random=1,
 648     //second option:leading=0
 649     Int_t axis=-1;
 650     if(fLeadingOrRandom==0)axis=0;
 651     else if (fLeadingOrRandom==1)axis= Int_t((highPtTracks)*gRandom->Rndm());
 652     else Printf("Wrong settings for event axis.");
 653     if(fDebug)Printf("Axis tracks has pT=%f", pt[pindex[axis]]);
 654
 655     //---------------------------------------
 656
 657
 658     if(ntracks>1){ // require at least two tracks (leading and prob. accosicates)
 659
 660       //EventAxisRandom track properties
 661       ptEventAxis  = pt [pindex[axis]];
 662       etaEventAxis = eta[pindex[axis]];
 663       phiEventAxis = phi[pindex[axis]];
 664
 665
 666       //track loop for event propoerties around event axis with pt>triggerPtCut
 667       //loop only over already accepted tracks except event axis
 668       if(ptEventAxis>fTriggerPtCut){
 669
 670         for (Int_t iTrack = 0; iTrack < ntracks; iTrack++) {
 671
 672           if(axis==iTrack)continue; // no double counting
 673
 674           ptOthers   = pt [pindex[iTrack]];
 675           etaOthers  = eta[pindex[iTrack]];
 676           phiOthers  = phi[pindex[iTrack]];
 677
 678
 679           if(ptOthers>fAssociatePtCut){ // only tracks which fullfill associate pt cut
 680
 681             Float_t dPhi=TMath::Abs(phiOthers-phiEventAxis);
 682             if(dPhi>TMath::Pi())      dPhi=2*TMath::Pi()-dPhi;
 683             Float_t dEta=etaOthers-etaEventAxis;
 684
 685             Float_t dphiplot = phiOthers-phiEventAxis;
 686             if(dphiplot>2*TMath::Pi()-1) dphiplot = dphiplot-2*TMath::Pi();
 687             else if(dphiplot<-1)dphiplot=dphiplot+2*TMath::Pi();
 688             fDPhiDEtaEventAxis[mode]->Fill(dphiplot, dEta);
 689
 690             if(ntracks<150){
 691               fDPhiEventAxisNchBin[mode][ntracks]->Fill(dPhi);
 692               if(ptOthers>fTriggerPtCut)
 693                 fDPhiEventAxisNchBinTrig[mode][ntracks]->Fill(dPhi);
 694             }
 695
 696             if(ntracklets<150){
 697               fDPhiEventAxisTrackletBin[mode][ntracklets]->Fill(dPhi);
 698               if(ptOthers>fTriggerPtCut)
 699                 fDPhiEventAxisTrackletBinTrig[mode][ntracklets]->Fill(dPhi);
 700             }
 701
 702           }//tracks fulfill assoc track cut
 703
 704         }// end track loop
 705
 706
 707         // fill histogram with number of tracks (pt>fAssociatePtCut) around event axis
 708         // how often is there a trigger particle at a certain Nch bin
 709         fTriggerNch[mode]->Fill(ntracks);
 710         fTriggerTracklet[mode]->Fill(ntracklets);
 711
 712       }//if track pt is at least trigger pt
 713
 714     } //if there are more than 1 track
 715
 716   }//if there is at least one high pt track
 717
 718   if(pindex){// clean up array memory used for TMath::Sort
 719     delete[] pindex;
 720     pindex=0;
 721   }
 722
 723   PostData(1, fHists);
 724
 725 }
 726
 727 //________________________________________________________________________
 728 void AliAnalysisTaskMinijet::Terminate(Option_t*)
 729 {
 730   //terminate function is called at the end
 731   //can be used to draw histograms etc.
 732
 733 }
 734
 735 //________________________________________________________________________
 736 void AliAnalysisTaskMinijet::CleanArrays(Float_t* pt, Float_t* eta, Float_t* phi,Int_t* numbers)
 737 {
 738   //clean up of memory used for arrays of track properties
 739
 740   if(pt){
 741     delete[] pt;
 742     pt=0;
 743   }
 744   if(eta){
 745     delete[] eta;
 746     eta=0;
 747   }
 748   if(phi){
 749     delete[] phi;
 750     phi=0;
 751   }
 752   if(numbers){
 753     delete[] numbers;
 754     numbers=0;
 755   }
 756
 757 }