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