PWG4/JetTasks/AliAnalysisTaskMinijet.cxx

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