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