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