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 <THnSparse.h>
   8 #include <TProfile.h>
   9 #include <TCanvas.h>
  10 #include "TRandom.h"
  11
  12 #include "AliVEvent.h"
  13 #include "AliVParticle.h"
  14
  15 #include "AliESDEvent.h"
  16 #include "AliESDtrack.h"
  17 #include "AliMultiplicity.h"
  18 #include "AliESDtrackCuts.h"
  19
  20 #include "AliAODEvent.h"
  21 #include "AliAODTrack.h"
  22 #include "AliAODMCParticle.h"
  23 #include "AliAODMCHeader.h"
  24
  25 #include "AliStack.h"
  26 #include "AliMCEvent.h"
  27 #include "AliMCParticle.h"
  28 #include "AliGenEventHeader.h"
  29
  30 #include "AliAnalysisManager.h"
  31 #include "AliInputEventHandler.h"
  32
  33 #include <vector>
  34 #include <algorithm>
  35 #include <functional>
  36 using namespace std;
  37
  38 #include "AliAnalysisTaskMinijet.h"
  39
  40 // Analysis task for two-particle correlations using all particles over pt threshold
  41 // pt_trig threshold for trigger particle (event axis) and pt_assoc for possible associated particles.
  42 // Extract mini-jet yield and fragmentation properties via Delta-Phi histograms of these correlations
  43 // post processing of analysis output via macro plot3and2Gaus.C
  44 // Can use ESD or AOD, reconstructed and Monte Carlo data as input
  45 // Author: Eva Sicking
  46
  47
  48 ClassImp(AliAnalysisTaskMinijet)
  49
  50 //________________________________________________________________________
  51   AliAnalysisTaskMinijet::AliAnalysisTaskMinijet(const char *name)
  52     : AliAnalysisTaskSE(name),
  53       fUseMC(kFALSE),
  54       fMcOnly(kFALSE),
  55       fBSign(0),
  56       fAnalysePrimOnly(kFALSE),// not used
  57       fPtMin(0.2),
  58       fPtMax(50.0),
  59       fCuts(0),
  60       fTriggerPtCut(0.7),
  61       fAssociatePtCut(0.4),
  62       fMode(0),
  63       fTriggerType(1),
  64       fFilterBit(128),
  65       fVertexZCut(10.),
  66       fEtaCut(0.9),
  67       fEtaCutSeed(0.9),
  68       fSelectParticles(1),
  69       fSelectParticlesAssoc(1),
  70       fESDEvent(0),
  71       fAODEvent(0),
  72       fNMcPrimAccept(0),
  73       fNRecAccept(0),
  74       fVzEvent(0),
  75       fMeanPtRec(0),
  76       fLeadingPtRec(0),
  77       fHists(0),
  78       fStep(0),
  79       fHistPt(0),
  80       fHistPtMC(0),
  81       fNmcNch(0),
  82       fPNmcNch(0),
  83       fChargedPi0(0)
  84 {
  85
  86   //Constructor
  87
  88   for(Int_t i = 0;i< 6;i++){
  89     fMapSingleTrig[i]         =  0;
  90     fMapPair[i]               =  0;
  91     fMapEvent[i]              =  0;
  92     fMapAll[i]                =  0;
  93
  94     fVertexZ[i]               =  0;
  95
  96     fNcharge[i]               =  0;
  97     fPt[i]                    =  0;
  98     fEta[i]                   =  0;
  99     fPhi[i]                   =  0;
 100     fDcaXY[i]                 =  0;
 101     fDcaZ[i]                  =  0;
 102
 103     fPtSeed[i]                =  0;
 104     fEtaSeed[i]               =  0;
 105     fPhiSeed[i]               =  0;
 106
 107     fPtOthers[i]              =  0;
 108     fEtaOthers[i]             =  0;
 109     fPhiOthers[i]             =  0;
 110     fPtEtaOthers[i]           =  0;
 111
 112     fPhiEta[i]                =  0;
 113
 114     fDPhiDEtaEventAxis[i]     =  0;
 115     fDPhiDEtaEventAxisSeeds[i]=  0;
 116     fTriggerNch[i]            =  0;
 117
 118     fTriggerNchSeeds[i]       =  0;
 119     fTriggerTracklet[i]       =  0;
 120
 121     fNch07Nch[i]              =  0;
 122     fPNch07Nch[i]             =  0;
 123
 124     fNch07Tracklet[i]         =  0;
 125     fNchTracklet[i]           =  0;
 126     fPNch07Tracklet[i]        =  0;
 127     fDPhiEventAxis[i]         =  0;
 128   }
 129   DefineOutput(1, TList::Class());
 130 }
 131
 132 //________________________________________________________________________
 133 AliAnalysisTaskMinijet::~AliAnalysisTaskMinijet()
 134 {
 135   // Destructor
 136
 137   if (fHists && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHists;
 138 }
 139
 140 //________________________________________________________________________
 141 void AliAnalysisTaskMinijet::UserCreateOutputObjects()
 142 {
 143   // Create histograms
 144   // Called once
 145   if(fDebug) Printf("In User Create Output Objects.");
 146
 147   fStep = new TH1F("fStep", "fStep", 10, -0.5, 9.5);
 148   fHistPt = new TH1F("fHistPt", "P_{T} distribution REC", 150, 0.1, 3.1);
 149   fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 150   fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 151   fHistPt->SetMarkerStyle(kFullCircle);
 152
 153   if (fUseMC) {
 154     fHistPtMC = new TH1F("fHistPtMC", "P_{T} distribution MC", 150, 0.1, 3.1);
 155     fHistPtMC->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 156     fHistPtMC->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 157     fHistPtMC->SetMarkerStyle(kFullCircle);
 158
 159     fNmcNch = new TH2F("fNmcNch", "fNmcNch", 100,-0.5,99.5,100,-0.5,99.5);
 160     fPNmcNch = new TProfile("pNmcNch", "pNmcNch", 100,-0.5,99.5);
 161
 162   }
 163
 164   fChargedPi0 = new TH2F("fChargedPi0", "fChargedPi0", 200, -0.5, 199.5, 200, -0.5, 199.5);
 165
 166
 167   //----------------------
 168   //bins for pt in THnSpare
 169   Double_t ptMin = 0.0, ptMax = 100.;
 170   Int_t nPtBins = 39;
 171   Double_t binsPt[]  = {0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8,
 172                         0.9, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0,
 173                         10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};
 174
 175   //  Int_t nPtBins = 100;
 176   //   Double_t ptMin = 1.e-2, ptMax = 100.;
 177   //   Double_t *binsPt = 0;
 178   //   binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 179
 180
 181   Double_t ptMin2 = 0.0, ptMax2 = 100.;
 182   Int_t nPtBins2 = 9;
 183   Double_t binsPt2[]  = {0.1, 0.4, 0.7, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 100.0};
 184
 185   //  Int_t nPtBins2 = 10;
 186   //   Double_t ptMin2 = 0.4, ptMax2 = 100.;
 187   //   Double_t *binsPt2 = 0;
 188   //   binsPt2 = CreateLogAxis(nPtBins2,ptMin2,ptMax2);
 189
 190   //3 dim matrix
 191   Int_t binsEffHisto[3]   = {Int_t(fEtaCut*20),  nPtBins,    150   };
 192   Double_t minEffHisto[3] = {-fEtaCut,           ptMin,        -0.5 };
 193   Double_t maxEffHisto[3] = {fEtaCut,            ptMax,        149.5 };
 194
 195   //5 dim matrix
 196   Int_t binsEffHisto5[6]   = {  nPtBins2,  nPtBins2,  Int_t(fEtaCut*10),    180,                   150 ,      2 };
 197   Double_t minEffHisto5[6] = {  ptMin2,    ptMin2,    -2*fEtaCut,          -0.5*TMath::Pi(),      -0.5 ,   -0.5 };
 198   Double_t maxEffHisto5[6] = {  ptMax2,    ptMax2,     2*fEtaCut,           1.5*TMath::Pi(),     149.5 ,    1.5 };
 199
 200
 201   //4 dim matrix
 202   Int_t binsEvent[4]   = {   150,        20,   50,  nPtBins };
 203   Double_t minEvent[4] = {  -0.5,       -10,    0,    ptMin };
 204   Double_t maxEvent[4] = { 149.5,        10,   10,    ptMax };
 205
 206   //3 dim matrix
 207   Int_t binsAll[3]   = {Int_t(fEtaCut*20),  nPtBins2,       150   };
 208   Double_t minAll[3] = {-fEtaCut,           ptMin2,        -0.5 };
 209   Double_t maxAll[3] = {fEtaCut,            ptMax2,        149.5 };
 210
 211   //--------------------
 212   TString dataType[2] ={"ESD", "AOD"};
 213   TString labels[6]={Form("%sAllAllMcNmc",dataType[fMode].Data()),
 214                      Form("%sTrigAllMcNmc",dataType[fMode].Data()),
 215                      Form("%sTrigAllMcNrec",dataType[fMode].Data()),
 216                      Form("%sTrigVtxMcNrec",dataType[fMode].Data()),
 217                      Form("%sTrigVtxRecMcPropNrec",dataType[fMode].Data()),
 218                      Form("%sTrigVtxRecNrec",dataType[fMode].Data())};
 219
 220
 221   for(Int_t i=0;i<6;i++){
 222
 223     fMapSingleTrig[i] = new THnSparseD(Form("fMapSingleTrig%s", labels[i].Data()),"eta:pt:Nrec",
 224                                        3,binsEffHisto,minEffHisto,maxEffHisto);
 225     fMapSingleTrig[i]->SetBinEdges(1,binsPt);
 226     fMapSingleTrig[i]->GetAxis(0)->SetTitle("#eta");
 227     fMapSingleTrig[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
 228     fMapSingleTrig[i]->GetAxis(2)->SetTitle("N_{rec}");
 229     fMapSingleTrig[i]->Sumw2();
 230
 231     fMapPair[i] = new THnSparseD(Form("fMapPair%s", labels[i].Data()),"pt_trig:pt_assoc:DeltaEta:DeltaPhi:Nrec:likesign",
 232                                  6,binsEffHisto5,minEffHisto5,maxEffHisto5);
 233     fMapPair[i]->SetBinEdges(0,binsPt2);
 234     fMapPair[i]->SetBinEdges(1,binsPt2);
 235     fMapPair[i]->GetAxis(0)->SetTitle("p_{T, trig} (GeV/c)");
 236     fMapPair[i]->GetAxis(1)->SetTitle("p_{T, assoc} (GeV/c)");
 237     fMapPair[i]->GetAxis(2)->SetTitle("#Delta #eta");
 238     fMapPair[i]->GetAxis(3)->SetTitle("#Delta #phi");
 239     fMapPair[i]->GetAxis(4)->SetTitle("N_{rec}");
 240     fMapPair[i]->GetAxis(5)->SetTitle("Like-sign or Unlike-sign");
 241     fMapPair[i]->Sumw2();
 242
 243
 244     fMapEvent[i] = new THnSparseD(Form("fMapEvent%s", labels[i].Data()),"Nrec:vertexZ:meanPt:leadingPt",
 245                                   4,binsEvent,minEvent,maxEvent);
 246     fMapEvent[i]->GetAxis(0)->SetTitle("N_{rec}");
 247     fMapEvent[i]->GetAxis(1)->SetTitle("z_{vertex} (cm)");
 248     fMapEvent[i]->GetAxis(2)->SetTitle("meanPt");
 249     fMapEvent[i]->SetBinEdges(3,binsPt);
 250     fMapEvent[i]->GetAxis(3)->SetTitle("leadingPt");
 251     fMapEvent[i]->Sumw2();
 252
 253     fMapAll[i] = new THnSparseD(Form("fMapAll%s", labels[i].Data()),"eta:pt:Nrec",
 254                                 3,binsAll,minAll,maxAll);
 255     fMapAll[i]->SetBinEdges(1,binsPt2);
 256     fMapAll[i]->GetAxis(0)->SetTitle("#eta");
 257     fMapAll[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
 258     fMapAll[i]->GetAxis(2)->SetTitle("N_{rec}");
 259     fMapAll[i]->Sumw2();
 260
 261
 262     fVertexZ[i]                  = new TH1F(Form("fVertexZ%s",labels[i].Data()),
 263                                             Form("fVertexZ%s",labels[i].Data()) ,
 264                                             200, -10., 10.);
 265     fPt[i]                       = new TH1F(Form("fPt%s",labels[i].Data()),
 266                                             Form("fPt%s",labels[i].Data()) ,
 267                                             200, 0., 50);
 268     fNcharge[i]                  = new TH1F(Form("fNcharge%s",labels[i].Data()),
 269                                             Form("fNcharge%s",labels[i].Data()) ,
 270                                             250, -0.5, 249.5);
 271     fEta[i]                      = new TH1F (Form("fEta%s",labels[i].Data()),
 272                                              Form("fEta%s",labels[i].Data()) ,
 273                                              100, -2., 2);
 274     fPhi[i]                      = new TH1F(Form("fPhi%s",labels[i].Data()),
 275                                             Form("fPhi%s",labels[i].Data()) ,
 276                                             360, 0.,2*TMath::Pi());
 277     fDcaXY[i]                    = new TH1F(Form("fDcaXY%s",labels[i].Data()),
 278                                             Form("fDcaXY%s",labels[i].Data()) ,
 279                                             200, -3,3);
 280     fDcaZ[i]                     = new TH1F(Form("fDcaZ%s",labels[i].Data()),
 281                                             Form("fDcaZ%s",labels[i].Data()) ,
 282                                             200, -10,10);
 283     fPtSeed[i]                       = new TH1F(Form("fPSeedt%s",labels[i].Data()),
 284                                                 Form("fPtSeed%s",labels[i].Data()) ,
 285                                                 500, 0., 50);
 286     fEtaSeed[i]                      = new TH1F (Form("fEtaSeed%s",labels[i].Data()),
 287                                                  Form("fEtaSeed%s",labels[i].Data()) ,
 288                                                  100, -2., 2);
 289     fPhiSeed[i]                      = new TH1F(Form("fPhiSeed%s",labels[i].Data()),
 290                                                 Form("fPhiSeed%s",labels[i].Data()) ,
 291                                                 360, 0.,2*TMath::Pi());
 292
 293     fPtOthers[i]                       = new TH1F(Form("fPOtherst%s",labels[i].Data()),
 294                                                   Form("fPtOthers%s",labels[i].Data()) ,
 295                                                   500, 0., 50);
 296     fEtaOthers[i]                      = new TH1F (Form("fEtaOthers%s",labels[i].Data()),
 297                                                    Form("fEtaOthers%s",labels[i].Data()) ,
 298                                                    100, -2., 2);
 299     fPhiOthers[i]                      = new TH1F(Form("fPhiOthers%s",labels[i].Data()),
 300                                                   Form("fPhiOthers%s",labels[i].Data()) ,
 301                                                   360, 0.,2*TMath::Pi());
 302     fPtEtaOthers[i]                      = new TH2F(Form("fPtEtaOthers%s",labels[i].Data()),
 303                                                     Form("fPtEtaOthers%s",labels[i].Data()) ,
 304                                                     500, 0., 50, 100, -1., 1);
 305     fPhiEta[i]                   = new TH2F(Form("fPhiEta%s",labels[i].Data()),
 306                                             Form("fPhiEta%s",labels[i].Data()) ,
 307                                             180, 0., 2*TMath::Pi(), 100, -1.,1.);
 308     fDPhiDEtaEventAxis[i]        = new TH2F(Form("fDPhiDEtaEventAxis%s",labels[i].Data()),
 309                                             Form("fDPhiDEtaEventAxis%s",labels[i].Data()) ,
 310                                             180, -0.5* TMath::Pi(), 1.5*TMath::Pi(), 9, -1.8,1.8);
 311     fTriggerNch[i]               = new TH1F(Form("fTriggerNch%s",labels[i].Data()),
 312                                             Form("fTriggerNch%s",labels[i].Data()) ,
 313                                             250, -0.5, 249.5);
 314     fTriggerNchSeeds[i]          = new TH2F(Form("fTriggerNchSeeds%s",labels[i].Data()),
 315                                             Form("fTriggerNchSeeds%s",labels[i].Data()) ,
 316                                             250, -0.5, 249.5, 100, -0.5, 99.5);
 317     fTriggerTracklet[i]          = new TH1F(Form("fTriggerTracklet%s",labels[i].Data()),
 318                                             Form("fTriggerTracklet%s",labels[i].Data()) ,
 319                                             250, -0.5, 249.5);
 320     fNch07Nch[i]                 = new TH2F(Form("fNch07Nch%s",labels[i].Data()),
 321                                             Form("fNch07Nch%s",labels[i].Data()) ,
 322                                             250, -2.5, 247.5,250, -2.5, 247.5);
 323     fPNch07Nch[i]                 = new TProfile(Form("pNch07Nch%s",labels[i].Data()),
 324                                                  Form("pNch07Nch%s",labels[i].Data()) ,
 325                                                  250, -2.5, 247.5);
 326     fNch07Tracklet[i]            = new TH2F(Form("fNch07Tracklet%s",labels[i].Data()),
 327                                             Form("fNch07Tracklet%s",labels[i].Data()) ,
 328                                             250, -2.5, 247.5,250, -2.5, 247.5);
 329     fNchTracklet[i]              = new TH2F(Form("fNchTracklet%s",labels[i].Data()),
 330                                             Form("fNchTracklet%s",labels[i].Data()) ,
 331                                             250, -2.5, 247.5,250, -2.5, 247.5);
 332     fPNch07Tracklet[i]            = new TProfile(Form("pNch07Tracklet%s",labels[i].Data()),
 333                                                  Form("pNch07Tracklet%s",labels[i].Data()) ,
 334                                                  250, -2.5, 247.5);
 335     fDPhiEventAxis[i]          = new TH1F(Form("fDPhiEventAxis%s",
 336                                                labels[i].Data()),
 337                                           Form("fDPhiEventAxis%s",
 338                                                labels[i].Data()) ,
 339                                           180, -0.5*TMath::Pi(), 1.5*TMath::Pi());
 340
 341   }
 342
 343   fHists = new TList();
 344   fHists->SetOwner();
 345
 346   fHists->Add(fStep);
 347   fHists->Add(fHistPt);
 348
 349   if(fUseMC){
 350     fHists->Add(fHistPtMC);
 351     fHists->Add(fNmcNch);
 352     fHists->Add(fPNmcNch);
 353   }
 354   fHists->Add(fChargedPi0);
 355
 356   for(Int_t i=0;i<6;i++){
 357     fHists->Add(fMapSingleTrig[i]);
 358     fHists->Add(fMapPair[i]);
 359     fHists->Add(fMapEvent[i]);
 360     fHists->Add(fMapAll[i]);
 361     fHists->Add(fVertexZ[i]);
 362     fHists->Add(fPt[i]);
 363     fHists->Add(fNcharge[i]);
 364     fHists->Add(fEta[i]);
 365     fHists->Add(fPhi[i]);
 366     fHists->Add(fDcaXY[i]);
 367     fHists->Add(fDcaZ[i]);
 368     fHists->Add(fPtSeed[i]);
 369     fHists->Add(fEtaSeed[i]);
 370     fHists->Add(fPhiSeed[i]);
 371     fHists->Add(fPtOthers[i]);
 372     fHists->Add(fEtaOthers[i]);
 373     fHists->Add(fPhiOthers[i]);
 374     fHists->Add(fPtEtaOthers[i]);
 375     fHists->Add(fPhiEta[i]);
 376     fHists->Add(fDPhiDEtaEventAxis[i]);
 377     fHists->Add(fTriggerNch[i]);
 378     fHists->Add(fTriggerNchSeeds[i]);
 379     fHists->Add(fTriggerTracklet[i]);
 380     fHists->Add(fNch07Nch[i]);
 381     fHists->Add(fPNch07Nch[i]);
 382     fHists->Add(fNch07Tracklet[i]);
 383     fHists->Add(fNchTracklet[i]);
 384     fHists->Add(fPNch07Tracklet[i]);
 385     fHists->Add(fDPhiEventAxis[i]);
 386   }
 387
 388   PostData(1, fHists);
 389
 390 }
 391
 392 //________________________________________________________________________
 393 void AliAnalysisTaskMinijet::UserExec(Option_t *)
 394 {
 395   // Main function, called for each event
 396   // Kinematics-only, ESD and AOD can be processed.
 397   // Data is read (ReadEventESD, ReadEventAOD...) and then analysed (Analyse).
 398   //  - in case of MC with full detector simulation, all correction steps(0-5) can be processed
 399   //  - for Data, only step 5 is performed
 400   //  - for kinematics-only, only step 0 is processed
 401   // step 5 =  Triggered events, reconstructed accepted vertex, reconstructed tracks,                    reconstructed multiplicity,
 402   // step 4 =  Triggered events, reconstructed accepted vertex, reconstructed tracks with MC properties, reconstructed multiplicity
 403   // step 3 =  Triggered events, reconstructed accepted vertex, mc primary particles,                    reconstructed multiplicity,
 404   // step 2 =  Triggered events, all                            mc primary particles,                    reconstructed multiplicity
 405   // step 1 =  Triggered events, all                            mc primary particles,                    true multiplicity
 406   // step 0 =  All events,       all                            mc primary particles,                    true multiplicity
 407
 408   // can be changed with new AOD definition -> AOD076
 409
 410   if(fDebug) Printf("UserExec: Event starts");
 411
 412   AliVEvent *event = InputEvent();
 413   if (!event) {
 414     Error("UserExec", "Could not retrieve event");
 415     return;
 416   }
 417   fBSign= event->GetMagneticField();
 418
 419   //get events, either ESD or AOD
 420   fESDEvent = dynamic_cast<AliESDEvent*> (InputEvent());
 421   fAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
 422
 423   vector<Float_t> pt;
 424   vector<Float_t> eta;
 425   vector<Float_t> phi;
 426   vector<Short_t> charge;
 427   vector<Float_t> px;
 428   vector<Float_t> py;
 429   vector<Float_t> pz;
 430   vector<Float_t> theta;
 431
 432
 433   //number of accepted tracks and tracklets
 434   Int_t ntracks = 0;  //return value for reading functions for ESD and AOD
 435   //Int_t ntracksRemove = 0;  //return value for reading functions for ESD and AOD
 436   vector<Int_t> nTracksTracklets; // [0]=nAccepted,1=ntracklets,2=nall(also neutral in case of mc,
 437                                   //for real nall=ncharged)
 438
 439
 440   if(!fAODEvent && !fESDEvent)return;
 441
 442   //=================== AOD ===============
 443
 444
 445
 446
 447   //reset values
 448   fNMcPrimAccept=0;// number of accepted primaries
 449   fNRecAccept=0;   // number of accepted tracks
 450
 451   // instead of task->SelectCollisionCandidate(mask) in AddTask macro
 452   Bool_t isSelected = (((AliInputEventHandler*)
 453                         (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 454                        ->IsEventSelected() &  fTriggerType);
 455
 456
 457   if(fDebug){
 458     Printf("IsSelected = %d", isSelected);
 459     Printf("CheckEvent(true)= %d", CheckEvent(true));
 460     Printf("CheckEvent(false)= %d", CheckEvent(false));
 461   }
 462
 463   //check trigger
 464   if(isSelected){ // has offline trigger
 465
 466     if(CheckEvent(true)){//step 5 = TrigVtxRecNrec, step 4 = TrigVtxRecMcPropNrec ,step 3 = TrigVtxMcNrec
 467
 468       if(!fMcOnly){
 469         //step 5 = TrigVtxRecNrec
 470
 471         // read tracks
 472         if(fESDEvent)     ntracks = ReadEventESD(pt, eta, phi, charge,nTracksTracklets, 5);
 473         else if(fAODEvent)ntracks = ReadEventAOD(pt, eta, phi, charge,nTracksTracklets, 5);
 474         else Printf("Fatal Error");
 475
 476         // analyse
 477         if(pt.size()){ //(internally ntracks=fNRecAccept)
 478           Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 5);
 479         }
 480
 481         if(fUseMC){
 482           // step 4 = TrigVtxRecMcPropNrec
 483
 484           // read tracks
 485           if(fESDEvent)       ntracks = ReadEventESDRecMcProp(pt, eta, phi, charge, nTracksTracklets, 4);
 486           else if(fAODEvent)  ntracks = ReadEventAODRecMcProp(pt, eta, phi, charge, nTracksTracklets, 4);
 487           else Printf("Fatal Error");
 488
 489           //analyse
 490           if(pt.size()){//(internally ntracks=fNRecAccept)
 491             Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 4);
 492           }
 493         }
 494
 495
 496       if(fUseMC){
 497         // step 3 = TrigVtxMcNrec
 498
 499         // read tracks
 500         if(fESDEvent)       ntracks = ReadEventESDMC(pt, eta, phi, charge, nTracksTracklets, 3);
 501         else if(fAODEvent)  ntracks = ReadEventAODMC(pt, eta, phi, charge, nTracksTracklets, 3);
 502         else Printf("Fatal Error");
 503
 504         // analyse
 505         if(pt.size()){//(internally ntracks=fNRecAccept)
 506           Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1],nTracksTracklets[2], 3);
 507         }
 508       }
 509       }
 510
 511     }//check event (true)
 512     if(fUseMC && !fMcOnly){
 513       //reset values
 514       fNMcPrimAccept=0;// number of accepted primaries
 515       fNRecAccept=0;   // number of accepted tracks
 516
 517       if(CheckEvent(false)){// all events, with and without reconstucted vertex
 518         if(fESDEvent){
 519           ntracks  = ReadEventESD  (pt, eta, phi, charge, nTracksTracklets, 2);//need to compute Nrec once more for all events
 520           ntracks  = ReadEventESDMC(pt, eta, phi, charge, nTracksTracklets, 1);//read tracks
 521         }
 522         else if(fAODEvent){
 523           ntracks  = ReadEventAOD  (pt, eta, phi, charge, nTracksTracklets, 2);//need to compute Nrec once more for all events
 524           ntracks  = ReadEventAODMC(pt, eta, phi, charge, nTracksTracklets, 1);//read tracks
 525         }
 526         else Printf("Fatal Error");
 527
 528         // analyse
 529         if(pt.size()){
 530           Analyse(pt, eta, phi, charge, fNRecAccept, nTracksTracklets[1],nTracksTracklets[2], 2); // step 2 = TrigAllMcNrec
 531
 532           Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 1);  // step 1 = TrigAllMcNmc
 533
 534           Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);  //first part of step 0 // step 0 = AllAllMcNmc
 535         }
 536
 537
 538       }
 539     }
 540   }// triggered event
 541
 542   else { // not selected by physics selection task = not triggered
 543     if(fUseMC && !fMcOnly){
 544       if(CheckEvent(false)){
 545
 546         //read tracks
 547         if(fESDEvent)      ntracks  = ReadEventESDMC(pt, eta, phi, charge, nTracksTracklets, 0);
 548         else if(fAODEvent) ntracks  = ReadEventAODMC(pt, eta, phi, charge, nTracksTracklets, 0);
 549         else Printf("Fatal Error");
 550
 551         //analyse
 552         if(pt.size()){
 553           Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);  //second part of step 0 // step 0 = AllAllMcNmc
 554         }
 555       }
 556     }
 557   }
 558
 559   if(fMcOnly){
 560     // read event
 561     if(fMode==0)       ntracks  = ReadEventESDMC(pt, eta, phi, charge, nTracksTracklets, 0);
 562     else if (fMode==1) ntracks  = ReadEventAODMC(pt, eta, phi, charge, nTracksTracklets, 0);
 563
 564     // analyse
 565     if(pt.size()){
 566       Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);
 567     }
 568   }
 569
 570
 571 }
 572
 573
 574 //________________________________________________________________________
 575 Int_t AliAnalysisTaskMinijet::ReadEventESD( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 576                                             vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
 577                                             vector<Int_t> &nTracksTracklets, const Int_t step)
 578 {
 579   // gives back the number of esd tracks and pointer to arrays with track
 580   // properties (pt, eta, phi)
 581   // Uses TPC tracks with SPD vertex from now on
 582
 583   ptArray.clear();
 584   etaArray.clear();
 585   phiArray.clear();
 586   chargeArray.clear();
 587   nTracksTracklets.clear();
 588
 589   const AliESDVertex*   vtxSPD   = fESDEvent->GetPrimaryVertexSPD(); // uses track or SPD vertexer
 590   fVertexZ[step]->Fill(vtxSPD->GetZ());
 591
 592   // Retreive the number of all tracks for this event.
 593   Int_t ntracks = fESDEvent->GetNumberOfTracks();
 594   if(fDebug>1)  Printf("all ESD tracks: %d", ntracks);
 595
 596   //first loop to check how many tracks are accepted
 597   //------------------
 598   Int_t nAcceptedTracks=0;
 599   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 600
 601     AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
 602     if (!esdTrack) {
 603       Error("ReadEventESD", "Could not receive track %d", iTracks);
 604       continue;
 605     }
 606
 607     // use TPC only tracks with non default SPD vertex
 608     AliESDtrack *track = AliESDtrackCuts::
 609       GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
 610     if(!track) continue;
 611     if(!fCuts->AcceptTrack(track)) {
 612       delete track;
 613       continue;// apply TPC track cuts before constrain to SPD vertex
 614     }
 615     if(track->Pt()>0.){
 616       // only constrain tracks above threshold
 617       AliExternalTrackParam exParam;
 618       // take the B-field from the ESD, no 3D fieldMap available at this point
 619       Bool_t relate = false;
 620       relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
 621                                         kVeryBig,&exParam);
 622       if(!relate){
 623         delete track;
 624         continue;
 625       }
 626       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
 627                  exParam.GetCovariance());
 628     }
 629     else{
 630       delete track;
 631       continue;// only if tracks have pt<=0
 632     }
 633
 634     if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax) {
 635       ptArray.push_back(track->Pt());
 636       etaArray.push_back(track->Eta());
 637       phiArray.push_back(track->Phi());
 638       chargeArray.push_back(track->Charge());
 639       ++nAcceptedTracks;
 640       fHistPt->Fill(track->Pt());
 641     }
 642
 643     // TPC only track memory needs to be cleaned up
 644     if(track)
 645       delete track;
 646
 647   }
 648
 649   //need to be checked
 650   if(nAcceptedTracks==0) return 0;
 651
 652   //tracklet loop
 653   Int_t ntrackletsAccept=0;
 654   AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
 655   Int_t ntracklets = mult->GetNumberOfTracklets();
 656   for(Int_t i=0;i< ntracklets;i++){
 657     if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
 658       ntrackletsAccept++;
 659     }
 660   }
 661   nTracksTracklets.push_back(nAcceptedTracks);
 662   nTracksTracklets.push_back(ntrackletsAccept);
 663   nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
 664                                               //where here also neutral particles are counted.
 665
 666
 667   fVzEvent=vtxSPD->GetZ(); // needed for correction map
 668   if(step==5 || step ==2) fNRecAccept=nAcceptedTracks;
 669   return fNRecAccept;
 670
 671
 672 }
 673
 674 //________________________________________________________________________
 675 Int_t AliAnalysisTaskMinijet::ReadEventESDRecMcProp( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 676                                                 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
 677                                                 vector<Int_t> &nTracksTracklets, const Int_t step)
 678 {
 679   // gives back the number of esd tracks and pointer to arrays with track
 680   // properties (pt, eta, phi) of mc particles if available
 681   // Uses TPC tracks with SPD vertex from now on
 682
 683   ptArray.clear();
 684   etaArray.clear();
 685   phiArray.clear();
 686   chargeArray.clear();
 687   nTracksTracklets.clear();
 688
 689
 690   AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
 691   if (!mcEvent) {
 692     Error("ReadEventESDRecMcProp", "Could not retrieve MC event");
 693     return 0;
 694   }
 695   AliStack* stack = MCEvent()->Stack();
 696   if(!stack) return 0;
 697
 698
 699   // Retreive the number of all tracks for this event.
 700   Int_t ntracks = fESDEvent->GetNumberOfTracks();
 701   if(fDebug>1)Printf("all ESD tracks: %d", ntracks);
 702
 703   const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
 704   fVertexZ[step]->Fill(vtxSPD->GetZ());
 705
 706   //track loop
 707   Int_t nAcceptedTracks=0;
 708   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 709
 710     AliVParticle *vtrack = fESDEvent->GetTrack(iTracks);
 711     AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
 712     if (!esdTrack) {
 713       Error("ReadEventESDRecMcProp", "Could not receive track %d", iTracks);
 714       continue;
 715     }
 716
 717     // use TPC only tracks with non default SPD vertex
 718     AliESDtrack *track = AliESDtrackCuts::
 719       GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
 720     if(!track) continue;
 721     if(!fCuts->AcceptTrack(track)) {
 722       delete track;
 723       continue;// apply TPC track cuts before constrain to SPD vertex
 724     }
 725     if(track->Pt()>0.){
 726       // only constrain tracks above threshold
 727       AliExternalTrackParam exParam;
 728       // take the B-field from the ESD, no 3D fieldMap available at this point
 729       Bool_t relate = false;
 730       relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
 731                                         kVeryBig,&exParam);
 732       if(!relate){
 733         delete track;
 734         continue;
 735       }
 736       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
 737                  exParam.GetCovariance());
 738     }
 739     else{
 740       delete track;
 741       continue;// only if tracks have pt<=0
 742     }
 743
 744     //count tracks, if available, use mc particle properties
 745     if(vtrack->GetLabel()<0){
 746       if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax){
 747         ptArray.push_back(track->Pt());
 748         etaArray.push_back(track->Eta());
 749         phiArray.push_back(track->Phi());
 750         chargeArray.push_back(track->Charge());
 751         ++nAcceptedTracks;
 752       }
 753     }
 754     else{
 755       TParticle *partOfTrack = stack->Particle(vtrack->GetLabel());
 756       if (TMath::Abs(partOfTrack->Eta())<fEtaCut && partOfTrack->Pt()>fPtMin && partOfTrack->Pt()<fPtMax) {
 757         ptArray.push_back(partOfTrack->Pt());
 758         etaArray.push_back(partOfTrack->Eta());
 759         phiArray.push_back(partOfTrack->Phi());
 760         chargeArray.push_back(vtrack->Charge());
 761         ++nAcceptedTracks;
 762       }
 763     }
 764
 765     // TPC only track memory needs to be cleaned up
 766     if(track)
 767       delete track;
 768
 769   }
 770
 771   if(nAcceptedTracks==0) return 0;
 772
 773   //tracklet loop
 774   Int_t ntrackletsAccept=0;
 775   AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
 776   Int_t ntracklets = mult->GetNumberOfTracklets();
 777   for(Int_t i=0;i< ntracklets;i++){
 778     if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
 779       ntrackletsAccept++;
 780     }
 781   }
 782
 783   nTracksTracklets.push_back(nAcceptedTracks);
 784   nTracksTracklets.push_back(ntrackletsAccept);
 785   nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
 786                                               //where here also neutral particles are counted.
 787
 788
 789   //get mc vertex for correction maps
 790   AliGenEventHeader*  header = MCEvent()->GenEventHeader();
 791   TArrayF mcV;
 792   header->PrimaryVertex(mcV);
 793   fVzEvent= mcV[2];
 794
 795   return fNRecAccept; // give back reconstructed value
 796
 797
 798 }
 799
 800
 801
 802
 803 //________________________________________________________________________
 804 Int_t AliAnalysisTaskMinijet::ReadEventESDMC(vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 805                                         vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
 806                                         vector<Int_t> &nTracksTracklets, const Int_t step)
 807 {
 808   // gives back the number of charged prim MC particle and pointer to arrays
 809   // with particle properties (pt, eta, phi)
 810
 811   ptArray.clear();
 812   etaArray.clear();
 813   phiArray.clear();
 814   chargeArray.clear();
 815   nTracksTracklets.clear();
 816
 817   fNMcPrimAccept=0;
 818
 819   AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
 820   if (!mcEvent) {
 821     Error("ReadEventESDMC", "Could not retrieve MC event");
 822     return 0;
 823   }
 824
 825   AliStack* stack = MCEvent()->Stack();
 826   if(!stack) return 0;
 827
 828   Int_t ntracks = mcEvent->GetNumberOfTracks();
 829   if(fDebug>1)Printf("MC particles: %d", ntracks);
 830
 831   //vertex
 832   AliGenEventHeader*  header = MCEvent()->GenEventHeader();
 833   TArrayF mcV;
 834   Float_t vzMC=0.;
 835   if(header){
 836     header->PrimaryVertex(mcV);
 837     vzMC = mcV[2];
 838     fVertexZ[step]->Fill(vzMC);
 839   }
 840
 841   //----------------------------------
 842   //first track loop to check how many chared primary tracks are available
 843   Int_t nChargedPrimaries=0;
 844   Int_t nAllPrimaries=0;
 845
 846   Int_t nPseudoTracklets=0;
 847   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 848     AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
 849     if (!track) {
 850       Error("ReadEventESDMC", "Could not receive track %d", iTracks);
 851       continue;
 852     }
 853
 854
 855     if(//count also charged particles in case of fSelectParticles==2 (only neutral)
 856        !SelectParticlePlusCharged(
 857                                   track->Charge(),
 858                                   track->PdgCode(),
 859                                   stack->IsPhysicalPrimary(track->Label())
 860                                   )
 861        )
 862       continue;
 863
 864     //count number of pseudo tracklets
 865     if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.0) ++nPseudoTracklets; //0.035
 866     //same cuts as on ESDtracks
 867     if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
 868
 869     //count all primaries
 870     ++nAllPrimaries;
 871     //count charged primaries
 872     if (track->Charge() != 0) ++nChargedPrimaries;
 873
 874     if(fDebug>2) Printf("PDG=%d, IsPrim=%d",  track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
 875
 876
 877   }
 878   //----------------------------------
 879
 880   if(fDebug>2){
 881     Printf("All in acceptance=%d",nAllPrimaries);
 882     Printf("Charged in acceptance =%d",nChargedPrimaries);
 883   }
 884
 885   fChargedPi0->Fill(nAllPrimaries,nChargedPrimaries);
 886
 887   if(nAllPrimaries==0) return 0;
 888   if(nChargedPrimaries==0) return 0;
 889
 890   //track loop
 891   //Int_t iChargedPiK=0;
 892   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 893     AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
 894     if (!track) {
 895       Error("ReadEventESDMC", "Could not receive track %d", iTracks);
 896       continue;
 897     }
 898
 899     if(!SelectParticle(
 900                        track->Charge(),
 901                        track->PdgCode(),
 902                        stack->IsPhysicalPrimary(track->Label())
 903                        )
 904        )
 905       continue;
 906
 907
 908     //same cuts as on ESDtracks
 909     if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin  || track->Pt()>fPtMax) continue;
 910
 911     if(fDebug>2) Printf("After: PDG=%d, IsPrim=%d",  track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
 912
 913
 914     fHistPtMC->Fill(track->Pt());
 915     //fills arrays with track properties
 916     ptArray.push_back(track->Pt());
 917     etaArray.push_back(track->Eta());
 918     phiArray.push_back(track->Phi());
 919     chargeArray.push_back(track->Charge());
 920
 921
 922   } //track loop
 923
 924   nTracksTracklets.push_back(nChargedPrimaries);
 925   nTracksTracklets.push_back(nPseudoTracklets);
 926   if(fSelectParticles!=2){
 927     nTracksTracklets.push_back(nAllPrimaries);
 928   }
 929   else{
 930     nTracksTracklets.push_back(nAllPrimaries-nChargedPrimaries); // only neutral
 931   }
 932
 933
 934   fNMcPrimAccept=nChargedPrimaries;
 935
 936   if(step==1){
 937     fNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
 938     fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
 939   }
 940
 941   fVzEvent= vzMC;
 942   return fNRecAccept;
 943
 944 }
 945
 946 //________________________________________________________________________
 947 Int_t AliAnalysisTaskMinijet::ReadEventAOD( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 948                                        vector<Float_t> &phiArray,  vector<Short_t> &chargeArray,
 949                                        vector<Int_t> &nTracksTracklets, const Int_t step)
 950 {
 951   // gives back the number of AOD tracks and pointer to arrays with track
 952   // properties (pt, eta, phi)
 953
 954   ptArray.clear();
 955   etaArray.clear();
 956   phiArray.clear();
 957   chargeArray.clear();
 958   nTracksTracklets.clear();
 959
 960   TClonesArray *mcArray=0x0;
 961   if(fAnalysePrimOnly){
 962     mcArray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
 963   }
 964
 965
 966   AliAODVertex* vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
 967   Double_t vzAOD=vertex->GetZ();
 968   fVertexZ[step]->Fill(vzAOD);
 969
 970   // Retreive the number of tracks for this event.
 971   Int_t ntracks = fAODEvent->GetNumberOfTracks();
 972   if(fDebug>1) Printf("AOD tracks: %d", ntracks);
 973
 974
 975   Int_t nAcceptedTracks=0;
 976   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 977     AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
 978     if (!track) {
 979       Error("ReadEventAOD", "Could not receive track %d", iTracks);
 980       continue;
 981     }
 982
 983     AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
 984
 985     //use only tracks from primaries
 986     if(fAnalysePrimOnly){
 987       if(vtrack->GetLabel()<0)continue;
 988       if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
 989     }
 990
 991     if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut
 992        && track->Pt()>fPtMin && track->Pt()<fPtMax){
 993
 994       nAcceptedTracks++;
 995
 996       // Printf("dca= %f", track->DCA());
 997       //save track properties in vector
 998       ptArray.push_back(track->Pt());
 999       etaArray.push_back(track->Eta());
1000       phiArray.push_back(track->Phi());
1001       chargeArray.push_back(track->Charge());
1002       fHistPt->Fill(track->Pt());
1003
1004     }
1005   }
1006   //need to check this option for MC
1007   if(nAcceptedTracks==0) return 0;
1008
1009
1010   //tracklet loop
1011   Int_t ntrackletsAccept=0;
1012   AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1013   for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1014     if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut){
1015       ++ntrackletsAccept;
1016     }
1017   }
1018
1019
1020   nTracksTracklets.push_back(nAcceptedTracks);
1021   nTracksTracklets.push_back(ntrackletsAccept);
1022   nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1023                                               //where here also neutral particles are counted.
1024
1025
1026   fVzEvent= vzAOD;
1027   if(step==5 || step==2)fNRecAccept = nAcceptedTracks; // needed for MC case //step5 = TrigVtxRecNrec
1028   return fNRecAccept; // at the moment, always return reconstructed multiplicity
1029
1030 }
1031
1032 //________________________________________________________________________
1033 Int_t AliAnalysisTaskMinijet::ReadEventAODRecMcProp( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1034                                                 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1035                                                 vector<Int_t> &nTracksTracklets, const Int_t step)
1036 {
1037   // gives back the number of AOD tracks and pointer to arrays with track
1038   // properties (pt, eta, phi)
1039
1040
1041   ptArray.clear();
1042   etaArray.clear();
1043   phiArray.clear();
1044   chargeArray.clear();
1045   nTracksTracklets.clear();
1046
1047
1048   // Retreive the number of tracks for this event.
1049   Int_t ntracks = fAODEvent->GetNumberOfTracks();
1050   if(fDebug>1) Printf("AOD tracks: %d", ntracks);
1051
1052
1053   //get array of mc particles
1054   TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1055     FindListObject(AliAODMCParticle::StdBranchName());
1056   if(!mcArray){
1057     Printf("No MC particle branch found");
1058     return kFALSE;
1059   }
1060
1061   AliAODVertex* vtx= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1062   Double_t vzAOD=vtx->GetZ();
1063   fVertexZ[step]->Fill(vzAOD);
1064
1065   Int_t nAcceptedTracks=0;
1066   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1067     AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
1068
1069     AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
1070
1071     if (!track) {
1072       Error("ReadEventAODRecMcProp", "Could not receive track %d", iTracks);
1073       continue;
1074     }
1075
1076     //use only tracks from primaries
1077     if(fAnalysePrimOnly){
1078       if(vtrack->GetLabel()<0)continue;
1079       if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
1080     }
1081
1082     if(track->TestFilterBit(fFilterBit) &&  TMath::Abs(track->Eta())<fEtaCut &&
1083        track->Pt()>fPtMin && track->Pt()<fPtMax){
1084
1085       nAcceptedTracks++;
1086
1087       //save track properties in vector
1088       if(vtrack->GetLabel()<0){ //fake tracks
1089         //      Printf("Fake track");
1090         //      continue;
1091         ptArray.push_back(track->Pt());
1092         etaArray.push_back(track->Eta());
1093         phiArray.push_back(track->Phi());
1094         chargeArray.push_back(track->Charge());
1095       }
1096       else{//mc properties
1097         AliAODMCParticle *partOfTrack = (AliAODMCParticle*)mcArray->At(vtrack->GetLabel());
1098
1099         ptArray.push_back(partOfTrack->Pt());
1100         etaArray.push_back(partOfTrack->Eta());
1101         phiArray.push_back(partOfTrack->Phi());
1102         chargeArray.push_back(vtrack->Charge());//partOfTrack?
1103       }
1104
1105     }
1106   }
1107   //need to check this option for MC
1108   if(nAcceptedTracks==0) return 0;
1109
1110   //tracklet loop
1111   Int_t ntrackletsAccept=0;
1112   AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1113   for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1114     if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut ){
1115       ++ntrackletsAccept;
1116     }
1117   }
1118
1119
1120   nTracksTracklets.push_back(nAcceptedTracks);
1121   nTracksTracklets.push_back(ntrackletsAccept);
1122   nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1123                                               //where here also neutral particles are counted.
1124
1125
1126   //check vertex (mc)
1127   AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1128     FindListObject(AliAODMCHeader::StdBranchName());
1129   Float_t vzMC = aodMCheader->GetVtxZ();
1130
1131   fVzEvent= vzMC;
1132   return fNRecAccept;//this is the rec value from step 5
1133
1134 }
1135
1136
1137
1138 //________________________________________________________________________
1139 Int_t AliAnalysisTaskMinijet::ReadEventAODMC( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1140                                          vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1141                                          vector<Int_t> &nTracksTracklets, const Int_t step)
1142 {
1143   // gives back the number of AOD MC particles and pointer to arrays with particle
1144   // properties (pt, eta, phi)
1145
1146   ptArray.clear();
1147   etaArray.clear();
1148   phiArray.clear();
1149   chargeArray.clear();
1150   nTracksTracklets.clear();
1151
1152   //check vertex
1153   AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1154     FindListObject(AliAODMCHeader::StdBranchName());
1155   Float_t vzMC = aodMCheader->GetVtxZ();
1156   fVertexZ[step]->Fill(vzMC);
1157
1158
1159   //retreive MC particles from event
1160   TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1161     FindListObject(AliAODMCParticle::StdBranchName());
1162   if(!mcArray){
1163     Printf("No MC particle branch found");
1164     return kFALSE;
1165   }
1166
1167   Int_t ntracks = mcArray->GetEntriesFast();
1168   if(fDebug>1)Printf("MC particles: %d", ntracks);
1169
1170
1171   // Track loop: chek how many particles will be accepted
1172   //Float_t vzMC=0.;
1173   Int_t nChargedPrim=0;
1174   Int_t nAllPrim=0;
1175   Int_t nPseudoTracklets=0;
1176   for (Int_t it = 0; it < ntracks; it++) {
1177     AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1178     if (!track) {
1179       Error("ReadEventAODMC", "Could not receive particle %d", it);
1180       continue;
1181     }
1182
1183     if(!SelectParticlePlusCharged(
1184                                   track->Charge(),
1185                                   track->PdgCode(),
1186                                   track->IsPhysicalPrimary()
1187                                   )
1188        )
1189       continue;
1190
1191     if(TMath::Abs(track->Eta())<fEtaCut && track->Pt()>0.0)++nPseudoTracklets; //0.035
1192     if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1193
1194     nAllPrim++;
1195     if(track->Charge()!=0) nChargedPrim++;
1196
1197   }
1198
1199
1200   if(nAllPrim==0) return 0;
1201   if(nChargedPrim==0) return 0;
1202
1203   //generate array with size of number of accepted tracks
1204   fChargedPi0->Fill(nAllPrim,nChargedPrim);
1205
1206
1207   // Track loop: fill arrays for accepted tracks
1208   // Int_t iChargedPrim=0;
1209   for (Int_t it = 0; it < ntracks; it++) {
1210     AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1211     if (!track) {
1212       Error("ReadEventAODMC", "Could not receive particle %d", it);
1213       continue;
1214     }
1215
1216     if(!SelectParticle(
1217                        track->Charge(),
1218                        track->PdgCode(),
1219                        track->IsPhysicalPrimary()
1220                        )
1221        )
1222       continue;
1223
1224     if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1225
1226     fHistPtMC->Fill(track->Pt());
1227     ptArray.push_back(track->Pt());
1228     etaArray.push_back(track->Eta());
1229     phiArray.push_back(track->Phi());
1230     chargeArray.push_back(track->Charge());
1231   }
1232
1233   nTracksTracklets.push_back(nChargedPrim);
1234   nTracksTracklets.push_back(nPseudoTracklets);
1235   if(fSelectParticles!=2){
1236     nTracksTracklets.push_back(nAllPrim);
1237   }
1238   else{
1239     nTracksTracklets.push_back(nAllPrim-nChargedPrim); // only neutral
1240   }
1241
1242
1243
1244   fVzEvent= vzMC;
1245   fNMcPrimAccept = nChargedPrim;
1246   if(step==1){ // step 1 = Trig All Mc Nmc
1247     fNmcNch->Fill( fNMcPrimAccept,fNRecAccept);
1248     fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1249   }
1250   return fNRecAccept; // rec value from step 5 or step 2
1251
1252
1253 }
1254
1255 //________________________________________________________________________
1256 void AliAnalysisTaskMinijet::Analyse(const vector<Float_t> &pt,
1257                                      const vector<Float_t> &eta,
1258                                      const vector<Float_t> &phi,
1259                                      const vector<Short_t> &charge,
1260                                      const Int_t ntracksCharged,
1261                                      const Int_t ntracklets,
1262                                      const Int_t nAll,
1263                                      const Int_t step)
1264 {
1265
1266   // analyse track properties (comming from either ESDs or AODs) in order to compute
1267   // mini jet activity (chared tracks) as function of charged multiplicity
1268
1269   fStep->Fill(step);
1270
1271   if(fDebug){
1272     Printf("Analysis Step=%d", step);
1273     if(fDebug>2){
1274       Printf("nAll=%d",nAll);
1275       Printf("nCharged=%d",ntracksCharged);
1276       for (Int_t i = 0; i < nAll; i++) {
1277         Printf("pt[%d]=%f",i,pt[i]);
1278       }
1279     }
1280   }
1281
1282   //calculation of mean pt for all tracks in case of step==0
1283   if(step==5 || step ==2){ // rec step
1284     Double_t meanPt=0.;
1285     Double_t leadingPt=0.;
1286     for (Int_t i = 0; i < nAll; i++) {
1287       if(pt[i]<0.01)continue;
1288       meanPt+=pt[i];
1289       if(leadingPt<pt[i])leadingPt=pt[i];
1290     }
1291     meanPt=meanPt/nAll;
1292     fMeanPtRec=meanPt;
1293     fLeadingPtRec=leadingPt;
1294   }
1295
1296   Double_t propEvent[4] = {ntracksCharged,fVzEvent,fMeanPtRec, fLeadingPtRec}; //vz: {rec, mc, mc}, meanPt and Nrec is always rec value
1297   fMapEvent[step]->Fill(propEvent);
1298
1299   fNcharge[step]->Fill(ntracksCharged);
1300
1301   Float_t ptEventAxis=0;  // pt event axis
1302   Float_t etaEventAxis=0; // eta event axis
1303   Float_t phiEventAxis=0; // phi event axis
1304   Short_t chargeEventAxis=0; // charge event axis
1305
1306   Float_t ptOthers  = 0; // pt others // for all other tracks around event axis -> see loop
1307   Float_t etaOthers = 0; // eta others
1308   Float_t phiOthers = 0; // phi others
1309   Short_t chargeOthers = 0; // charge others
1310
1311   Int_t   *pindexInnerEta  = new Int_t  [nAll+1];
1312   Float_t *ptInnerEta      = new Float_t[nAll+1];
1313
1314
1315
1316   for (Int_t i = 0; i < nAll; i++) {
1317
1318     if(pt[i]<0.01)continue;
1319
1320     //fill single particle correction for first step of pair correction
1321     Double_t propAll[3] = {eta[i],pt[i],ntracksCharged };
1322     fMapAll[step]->Fill(propAll);
1323
1324
1325     //filling of simple check plots
1326     if(pt[i]<0.7) continue;
1327     fPt[step]    -> Fill( pt[i]);
1328     fEta[step]   -> Fill(eta[i]);
1329     fPhi[step]   -> Fill(phi[i]);
1330     fPhiEta[step]-> Fill(phi[i], eta[i]);
1331
1332     pindexInnerEta[i]=0; //set all values to zero
1333     //fill new array for eta check
1334     ptInnerEta[i]= pt[i];
1335
1336   }
1337
1338
1339
1340   // define event axis: leading or random track (with pt>fTriggerPtCut)
1341   // ---------------------------------------
1342   Int_t highPtTracks=0;
1343   Int_t highPtTracksInnerEta=0;
1344   // Double_t highPtTracksInnerEtaCorr=0;
1345   Int_t mult09=0;
1346
1347   //count high pt tracks and high pt tracks in acceptance for seeds
1348   for(Int_t i=0;i<nAll;i++){
1349
1350     if(pt[i]<0.01)continue;
1351
1352     if(TMath::Abs(eta[i])<0.9){
1353       mult09++;
1354     }
1355
1356     if(pt[i]>fTriggerPtCut) {
1357       highPtTracks++;
1358     }
1359
1360     // seed should be place in middle of acceptance, that complete cone is in acceptance
1361     if(pt[i]>fTriggerPtCut && TMath::Abs(eta[i])<fEtaCutSeed && charge[i]!=0){
1362
1363       // Printf("eta=%f", eta[i]);
1364       highPtTracksInnerEta++;
1365
1366     }
1367     else{
1368       ptInnerEta[i]=0;
1369     }
1370   }
1371
1372
1373   //sort array in order to get highest pt tracks first
1374   //index can be computed with pindexInnerEta[number]
1375   if(nAll) TMath::Sort(nAll, ptInnerEta, pindexInnerEta, kTRUE);
1376
1377   //     plot of multiplicity distributions
1378   fNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1379   fPNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1380
1381   if(ntracklets){
1382     fNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1383     fNchTracklet[step]->Fill(ntracklets, ntracksCharged);
1384     fPNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1385   }
1386
1387   //analysis can only be performed with event axis, defined by high pt track
1388
1389
1390   if(highPtTracks>0 && highPtTracksInnerEta>0){
1391
1392     // build pairs in two track loops
1393     // loop over all possible trigger particles (defined by pt_trig and eta_acceptance)
1394     for(Int_t axis=0;(axis<nAll) && (pt[pindexInnerEta[axis]]>=fTriggerPtCut); axis++){
1395
1396       //EventAxisRandom track properties
1397       ptEventAxis  = pt [pindexInnerEta[axis]];
1398       etaEventAxis = eta[pindexInnerEta[axis]];
1399       phiEventAxis = phi[pindexInnerEta[axis]];
1400       chargeEventAxis = charge[pindexInnerEta[axis]];
1401       fPtSeed[step]    -> Fill( ptEventAxis);
1402       fEtaSeed[step]   -> Fill(etaEventAxis);
1403       fPhiSeed[step]   -> Fill(phiEventAxis);
1404
1405
1406       Double_t prop[3] = {etaEventAxis,ptEventAxis,ntracksCharged };
1407       fMapSingleTrig[step]->Fill(prop);
1408
1409       //associated tracks
1410       for (Int_t iTrack = axis+1; iTrack < nAll; iTrack++) {
1411
1412         if(pt[pindexInnerEta[iTrack]]<fAssociatePtCut) continue;
1413
1414         if(fSelectParticlesAssoc==1){
1415           if(charge[pindexInnerEta[iTrack]]==0)continue;
1416         }
1417         if(fSelectParticlesAssoc==2){
1418           if(charge[pindexInnerEta[iTrack]]!=0)continue;
1419         }
1420
1421
1422         ptOthers   = pt [pindexInnerEta[iTrack]];
1423         etaOthers  = eta[pindexInnerEta[iTrack]];
1424         phiOthers  = phi[pindexInnerEta[iTrack]];
1425         chargeOthers = charge[pindexInnerEta[iTrack]];
1426
1427
1428         //plot only properties of tracks with pt>ptassoc
1429         fPtOthers[step]    -> Fill( ptOthers);
1430         fEtaOthers[step]   -> Fill(etaOthers);
1431         fPhiOthers[step]   -> Fill(phiOthers);
1432         fPtEtaOthers[step]   -> Fill(ptOthers, etaOthers);
1433
1434         //      if(fDebug>2)Printf("%f, %f", pt[pindexInnerEta[axis]], pt[pindexInnerEta[iTrack]]);
1435
1436         Float_t dPhi = phiOthers-phiEventAxis;
1437         if(dPhi>1.5*TMath::Pi()) dPhi = dPhi-2*TMath::Pi();
1438         else if(dPhi<-0.5*TMath::Pi())dPhi=dPhi+2*TMath::Pi();
1439         Float_t dEta=etaOthers-etaEventAxis;
1440
1441
1442         fDPhiDEtaEventAxis[step]->Fill(dPhi, dEta);
1443         fDPhiEventAxis[step]->Fill(dPhi);
1444
1445         //check outliers
1446         if(ptEventAxis< 0.4 || ptEventAxis > 100) Printf("particles out of range pt");
1447         if(ntracksCharged<0 || ntracksCharged>150) Printf("particles out of range ncharge");
1448         if(TMath::Abs(dEta)>2*fEtaCut) {
1449           Printf("particles out of range dEta");
1450           Printf("eta1=%f, eta2=%f", etaOthers, etaEventAxis);
1451           Printf("step=%d",step);
1452         }
1453         if(dPhi<-0.5*TMath::Pi() || dPhi>1.5*TMath::Pi()){
1454           Printf("particles out of range dPhi");
1455           Printf("phi1=%f, phi2=%f", phiOthers, phiEventAxis);
1456         }
1457
1458         Bool_t isLikeSign = CheckLikeSign(chargeEventAxis, chargeOthers);
1459
1460         Double_t prop6[6] = {ptEventAxis,ptOthers,dEta,dPhi,ntracksCharged, isLikeSign };
1461         fMapPair[step]->Fill(prop6);
1462
1463       }// end of inner track loop
1464
1465     } //end of outer track loop
1466
1467     fTriggerNch[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1468     fTriggerNchSeeds[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1469     fTriggerTracklet[step]->Fill(ntracklets);
1470
1471
1472   }//if there is at least one high pt track
1473
1474
1475   if(pindexInnerEta){// clean up array memory used for TMath::Sort
1476     delete[] pindexInnerEta;
1477     pindexInnerEta=0;
1478   }
1479
1480   if(ptInnerEta){// clean up array memory used for TMath::Sort
1481     delete[] ptInnerEta;
1482     ptInnerEta=0;
1483   }
1484
1485 }
1486
1487
1488
1489 //________________________________________________________________________
1490 void AliAnalysisTaskMinijet::Terminate(Option_t*)
1491 {
1492   //terminate function is called at the end
1493   //can be used to draw histograms etc.
1494
1495
1496 }
1497
1498 //________________________________________________________________________
1499 Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(const Short_t charge, const Int_t pdg, Bool_t prim)
1500 {
1501   //selection of mc particle
1502   //fSelectParticles=0: use charged primaries and pi0 and k0
1503   //fSelectParticles=1: use only charged primaries
1504   //fSelectParticles=2: use only pi0 and k0
1505
1506   if(fSelectParticles==0 || fSelectParticles==2){ // in case of 2: need to count also charged particles
1507     if(charge==0){
1508       if(!(pdg==111||pdg==130||pdg==310))
1509         return false;
1510     }
1511     else{// charge !=0
1512       if(!prim)
1513         return false;
1514     }
1515   }
1516
1517   else if(fSelectParticles==1){
1518     if (charge==0 || !prim){
1519       return false;
1520     }
1521   }
1522
1523   else{
1524     Printf("Error: wrong selection of charged/pi0/k0");
1525     return 0;
1526   }
1527
1528   return true;
1529
1530 }
1531
1532 //________________________________________________________________________
1533 Bool_t AliAnalysisTaskMinijet::SelectParticle(const Short_t charge, const Int_t pdg, const Bool_t prim)
1534 {
1535   //selection of mc particle
1536   //fSelectParticles=0: use charged primaries and pi0 and k0
1537   //fSelectParticles=1: use only charged primaries
1538   //fSelectParticles=2: use only pi0 and k0
1539
1540   if(fSelectParticles==0){
1541     if(charge==0){
1542       if(!(pdg==111||pdg==130||pdg==310))
1543         return false;
1544     }
1545     else{// charge !=0
1546       if(!prim)
1547         return false;
1548     }
1549   }
1550   else if (fSelectParticles==1){
1551     if (charge==0 || !prim){
1552       return false;
1553     }
1554   }
1555   else if(fSelectParticles==2){
1556     if(!(pdg==111||pdg==130||pdg==310))
1557       return false;
1558   }
1559
1560   return true;
1561
1562 }
1563
1564 //________________________________________________________________________
1565 Bool_t AliAnalysisTaskMinijet::CheckEvent(const Bool_t recVertex)
1566 {
1567   // This function tests the quality of an event (ESD/AOD) (rec and/or mc part)
1568   // recVertex=false:  check if Mc events and stack is there, Nmc>0
1569   // recVertex=false: " + check if there is a good, reconstructed SPD vertex
1570   // defined by |z|<fVertexCut(10cm), Contributer>0, no PileUpFromSPD(3,0,8)
1571
1572
1573   if(fMode==0){//esd
1574
1575     //mc
1576     if(fUseMC){
1577
1578       //mc event
1579       AliMCEvent *mcEvente = (AliMCEvent*) MCEvent();
1580       if (!mcEvente) {
1581         Error("CheckEvent", "Could not retrieve MC event");
1582         return false;
1583       }
1584
1585       //stack
1586       AliStack* stackg = MCEvent()->Stack();
1587       if(!stackg) return false;
1588       Int_t ntracksg = mcEvente->GetNumberOfTracks();
1589       if(ntracksg<0) return false;
1590
1591       //vertex
1592       //AliGenEventHeader*  headerg = MCEvent()->GenEventHeader();
1593       //TArrayF mcVg;
1594       //headerg->PrimaryVertex(mcVg);
1595       //  if(TMath::Abs(mcVg[0])<1e-8 && TMath::Abs(mcVg[0])<1e-8 &&
1596       // TMath::Abs(mcVg[0])<1e-8) return false;
1597       // Float_t vzMCg = mcVg[2];
1598       // if(TMath::Abs(vzMCg)>fVertexZCut) return false;
1599     }
1600
1601     //rec
1602     if(recVertex==true){
1603       if(fESDEvent->IsPileupFromSPD(3,0,8))return false;
1604
1605       //rec vertex
1606       // const AliESDVertex*    vertexESDg   = fESDEvent->GetPrimaryVertex(); // uses track or SPD vertexer
1607       // if(!vertexESDg) return false;
1608       // if(vertexESDg->GetNContributors()<=0)return false;
1609       // Float_t fVzg= vertexESDg->GetZ();
1610       // if(TMath::Abs(fVzg)>fVertexZCut) return false;
1611
1612       //rec spd vertex
1613       const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
1614       if(!vtxSPD) return false;
1615       if(vtxSPD->GetNContributors()<=0)return false;
1616       Float_t fVzSPD= vtxSPD->GetZ();
1617       if(TMath::Abs(fVzSPD)>fVertexZCut) return false;
1618     }
1619     return true;
1620
1621   }
1622
1623
1624   else if(fMode==1){ //aod
1625
1626     if(fUseMC){
1627       //mc
1628       //   AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1629       //         FindListObject(AliAODMCHeader::StdBranchName());
1630       //        Float_t vzMC = aodMCheader->GetVtxZ();
1631       //        if(TMath::Abs(vzMC)>fVertexZCut) return false;
1632
1633       //retreive MC particles from event
1634       TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1635         FindListObject(AliAODMCParticle::StdBranchName());
1636       if(!mcArray){
1637         Printf("No MC particle branch found");
1638         return false;
1639       }
1640     }
1641
1642     //rec
1643     if(recVertex==true){
1644       if(fAODEvent->IsPileupFromSPD(3,0.8))return false;
1645
1646       //      AliAODVertex*     vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertex();
1647       //      if(!vertex) return false;
1648       //      if(vertex->GetNContributors()<=0) return false;
1649       //      Double_t vzAOD=vertex->GetZ();
1650       //      if(TMath::Abs(vzAOD)<1e-9) return false;
1651       //      if(TMath::Abs(vzAOD)>fVertexZCut) return false;
1652
1653       AliAODVertex*     vertexSPD= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();
1654       if(!vertexSPD) return false;
1655       if(vertexSPD->GetNContributors()<=0) return false;
1656       Double_t vzSPD=vertexSPD->GetZ();
1657       //if(TMath::Abs(vzSPD)<1e-9) return false;
1658       if(TMath::Abs(vzSPD)>fVertexZCut) return false;
1659     }
1660     return true;
1661
1662   }
1663
1664   else {
1665     Printf("Wrong mode!");
1666     return false;
1667   }
1668
1669 }
1670
1671 //_____________________________________________________________________________
1672 const Double_t * AliAnalysisTaskMinijet::CreateLogAxis(const Int_t nbins,
1673                                                        const Double_t xmin,
1674                                                        const Double_t xmax)
1675 {
1676   // returns pointer to an array with can be used to build a logarithmic axis
1677   // it is user responsibility to delete the array
1678
1679   Double_t logxmin = TMath::Log10(xmin);
1680   Double_t logxmax = TMath::Log10(xmax);
1681   Double_t binwidth = (logxmax-logxmin)/nbins;
1682
1683   Double_t *xbins =  new Double_t[nbins+1];
1684
1685   xbins[0] = xmin;
1686   for (Int_t i=1;i<=nbins;i++) {
1687     xbins[i] = xmin + TMath::Power(10,logxmin+i*binwidth);
1688   }
1689
1690   return xbins;
1691 }
1692
1693 //_____________________________________________________________________________
1694 Bool_t AliAnalysisTaskMinijet::CheckLikeSign(const Short_t chargeEventAxis,
1695                                              const Short_t chargeOthers)
1696 {
1697   // compute if charge of two particles/tracks has same sign or different sign
1698
1699   if(fDebug>2) Printf("Charge1=%d, Charge2=%d",chargeEventAxis,chargeOthers);
1700
1701   if(chargeEventAxis<0){
1702     if(chargeOthers<0){
1703       return true;
1704     }
1705     else if(chargeOthers>0){
1706       return false;
1707     }
1708   }
1709
1710   else if(chargeEventAxis>0){
1711     if(chargeOthers>0){
1712       return true;
1713     }
1714     else if(chargeOthers<0){
1715       return false;
1716     }
1717   }
1718
1719   else{
1720     Printf("Error: Charge not lower nor higher as zero");
1721     return false;
1722   }
1723
1724   Printf("Error: Check values of Charge");
1725   return false;
1726 }
1727