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     //reset values
 476     fNMcPrimAccept=0;// number of accepted primaries
 477     fNRecAccept=0;   // number of accepted tracks
 478
 479     if(CheckEvent(false)){// all events, with and without reconstucted vertex
 480       ntracks  = LoopAOD  (pt, eta, phi, charge, nTracksTracklets, 2);//need to compute Nrec once more for all events
 481       ntracks  = LoopAODMC(pt, eta, phi, charge, nTracksTracklets, 1);//read tracks
 482       if(pt.size()){
 483         Analyse(pt, eta, phi, charge, fNRecAccept, nTracksTracklets[1],nTracksTracklets[2], 2); // step 2 = TrigAllMcNrec
 484
 485         Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 1);  // step 1 = TrigAllMcNmc
 486       }
 487
 488       //step 0 (include not triggered events) is not possible with AODs generated before October 2011
 489       //step 0 can be implemented for new AODs
 490
 491     }
 492   }//AOD loop
 493
 494
 495   //=================== ESD ===============
 496
 497
 498   if(fESDEvent){//ESD loop
 499
 500     //reset values
 501     fNMcPrimAccept=0;// number of accepted primaries
 502     fNRecAccept=0;   // number of accepted tracks
 503
 504     // instead of task->SelectCollisionCandidate(mask) in AddTask macro
 505     Bool_t isSelected = (((AliInputEventHandler*)
 506                           (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 507                          ->IsEventSelected() &    AliVEvent::kMB);
 508
 509
 510     if(fDebug){
 511       Printf("IsSelected = %d", isSelected);
 512       Printf("CheckEvent(true)= %d", CheckEvent(true));
 513       Printf("CheckEvent(false)= %d", CheckEvent(false));
 514     }
 515
 516     //check trigger
 517     if(isSelected){ // has offline trigger
 518
 519       if(CheckEvent(true)){//step 5 = TrigVtxRecNrec, step 4 = TrigVtxRecMcPropNrec ,step 3 = TrigVtxMcNrec
 520
 521         if(!fMcOnly){
 522           //step 5 = TrigVtxRecNrec
 523           ntracks = LoopESD(pt, eta, phi, charge,nTracksTracklets, 5);//read tracks
 524           if(pt.size()){ //(internally ntracks=fNRecAccept)
 525             Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 5);//analyse
 526           }
 527
 528           if(fUseMC){
 529             // step 4 = TrigVtxRecMcPropNrec
 530             ntracks = LoopESDRecMcProp(pt, eta, phi, charge, nTracksTracklets, 4);//read tracks
 531             if(pt.size()){//(internally ntracks=fNRecAccept)
 532               Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1], nTracksTracklets[2], 4);//analyse
 533             }
 534           }
 535         }
 536
 537         if(fUseMC){
 538           // step 3 = TrigVtxMcNrec
 539           ntracks = LoopESDMC(pt, eta, phi, charge, nTracksTracklets, 3);//read tracks
 540           if(pt.size()){//(internally ntracks=fNRecAccept)
 541             Analyse(pt, eta, phi, charge, ntracks, nTracksTracklets[1],nTracksTracklets[2], 3);//analyse
 542           }
 543         }
 544
 545
 546       }//check event (true)
 547
 548       //reset values
 549       fNMcPrimAccept=0;// number of accepted primaries
 550       fNRecAccept=0;   // number of accepted tracks
 551
 552       if(CheckEvent(false)){// all events, with and without reconstucted vertex
 553         ntracks  = LoopESD  (pt, eta, phi, charge, nTracksTracklets, 2);//need to compute Nrec once more for all events
 554         ntracks  = LoopESDMC(pt, eta, phi, charge, nTracksTracklets, 1);//read tracks
 555         if(pt.size()){
 556           Analyse(pt, eta, phi, charge, fNRecAccept, nTracksTracklets[1],nTracksTracklets[2], 2); // step 2 = TrigAllMcNrec
 557
 558           Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 1);  // step 1 = TrigAllMcNmc
 559
 560           Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);  //first part of step 0 // step 0 = AllAllMcNmc
 561         }
 562
 563
 564       }
 565     }// triggered event
 566
 567     else { // not selected by physics selection task = not triggered
 568       if(CheckEvent(false)){
 569         ntracks  = LoopESDMC(pt, eta, phi, charge, nTracksTracklets, 0);//read tracks
 570         if(pt.size())Analyse(pt, eta, phi, charge, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);  //second part of step 0 // step 0 = AllAllMcNmc
 571       }
 572     }
 573
 574
 575   }//ESD loop
 576
 577
 578
 579
 580
 581 }
 582
 583
 584 //________________________________________________________________________
 585 Int_t AliAnalysisTaskMinijet::LoopESD( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 586                                        vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
 587                                        vector<Int_t> &nTracksTracklets, const Int_t step)
 588 {
 589   // gives back the number of esd tracks and pointer to arrays with track
 590   // properties (pt, eta, phi)
 591   // Uses TPC tracks with SPD vertex from now on
 592
 593   ptArray.clear();
 594   etaArray.clear();
 595   phiArray.clear();
 596   chargeArray.clear();
 597   nTracksTracklets.clear();
 598
 599   const AliESDVertex*   vtxSPD   = fESDEvent->GetPrimaryVertexSPD(); // uses track or SPD vertexer
 600   fVertexZ[step]->Fill(vtxSPD->GetZ());
 601
 602   // Retreive the number of all tracks for this event.
 603   Int_t ntracks = fESDEvent->GetNumberOfTracks();
 604   if(fDebug>1)  Printf("all ESD tracks: %d", ntracks);
 605
 606   //first loop to check how many tracks are accepted
 607   //------------------
 608   Int_t nAcceptedTracks=0;
 609   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 610
 611     AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
 612     if (!esdTrack) {
 613       Error("LoopESD", "Could not receive track %d", iTracks);
 614       continue;
 615     }
 616
 617     // use TPC only tracks with non default SPD vertex
 618     AliESDtrack *track = AliESDtrackCuts::
 619       GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
 620     if(!track) continue;
 621     if(!fCuts->AcceptTrack(track)) {
 622       delete track;
 623       continue;// apply TPC track cuts before constrain to SPD vertex
 624     }
 625     if(track->Pt()>0.){
 626       // only constrain tracks above threshold
 627       AliExternalTrackParam exParam;
 628       // take the B-field from the ESD, no 3D fieldMap available at this point
 629       Bool_t relate = false;
 630       relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
 631                                         kVeryBig,&exParam);
 632       if(!relate){
 633         delete track;
 634         continue;
 635       }
 636       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
 637                  exParam.GetCovariance());
 638     }
 639     else{
 640       delete track;
 641       continue;// only if tracks have pt<=0
 642     }
 643
 644     if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax) {
 645       ptArray.push_back(track->Pt());
 646       etaArray.push_back(track->Eta());
 647       phiArray.push_back(track->Phi());
 648       chargeArray.push_back(track->Charge());
 649       ++nAcceptedTracks;
 650       fHistPt->Fill(track->Pt());
 651     }
 652
 653     // TPC only track memory needs to be cleaned up
 654     if(track)
 655       delete track;
 656
 657   }
 658
 659   //need to be checked
 660   if(nAcceptedTracks==0) return 0;
 661
 662   //tracklet loop
 663   Int_t ntrackletsAccept=0;
 664   AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
 665   Int_t ntracklets = mult->GetNumberOfTracklets();
 666   for(Int_t i=0;i< ntracklets;i++){
 667     if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
 668       ntrackletsAccept++;
 669     }
 670   }
 671   nTracksTracklets.push_back(nAcceptedTracks);
 672   nTracksTracklets.push_back(ntrackletsAccept);
 673   nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
 674                                               //where here also neutral particles are counted.
 675
 676
 677   fVzEvent=vtxSPD->GetZ(); // needed for correction map
 678   if(step==5 || step ==2) fNRecAccept=nAcceptedTracks;
 679   return fNRecAccept;
 680
 681
 682 }
 683
 684 //________________________________________________________________________
 685 Int_t AliAnalysisTaskMinijet::LoopESDRecMcProp( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 686                                                 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
 687                                                 vector<Int_t> &nTracksTracklets, const Int_t step)
 688 {
 689   // gives back the number of esd tracks and pointer to arrays with track
 690   // properties (pt, eta, phi) of mc particles if available
 691   // Uses TPC tracks with SPD vertex from now on
 692
 693   ptArray.clear();
 694   etaArray.clear();
 695   phiArray.clear();
 696   chargeArray.clear();
 697   nTracksTracklets.clear();
 698
 699
 700   AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
 701   if (!mcEvent) {
 702     Error("LoopESDRecMcProp", "Could not retrieve MC event");
 703     return 0;
 704   }
 705   AliStack* stack = MCEvent()->Stack();
 706   if(!stack) return 0;
 707
 708
 709   // Retreive the number of all tracks for this event.
 710   Int_t ntracks = fESDEvent->GetNumberOfTracks();
 711   if(fDebug>1)Printf("all ESD tracks: %d", ntracks);
 712
 713   const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
 714   fVertexZ[step]->Fill(vtxSPD->GetZ());
 715
 716   //track loop
 717   Int_t nAcceptedTracks=0;
 718   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 719
 720     AliVParticle *vtrack = fESDEvent->GetTrack(iTracks);
 721     AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
 722     if (!esdTrack) {
 723       Error("LoopESDRecMcProp", "Could not receive track %d", iTracks);
 724       continue;
 725     }
 726
 727     // use TPC only tracks with non default SPD vertex
 728     AliESDtrack *track = AliESDtrackCuts::
 729       GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
 730     if(!track) continue;
 731     if(!fCuts->AcceptTrack(track)) {
 732       delete track;
 733       continue;// apply TPC track cuts before constrain to SPD vertex
 734     }
 735     if(track->Pt()>0.){
 736       // only constrain tracks above threshold
 737       AliExternalTrackParam exParam;
 738       // take the B-field from the ESD, no 3D fieldMap available at this point
 739       Bool_t relate = false;
 740       relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
 741                                         kVeryBig,&exParam);
 742       if(!relate){
 743         delete track;
 744         continue;
 745       }
 746       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
 747                  exParam.GetCovariance());
 748     }
 749     else{
 750       delete track;
 751       continue;// only if tracks have pt<=0
 752     }
 753
 754     //count tracks, if available, use mc particle properties
 755     if(vtrack->GetLabel()<0){
 756       if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax){
 757         ptArray.push_back(track->Pt());
 758         etaArray.push_back(track->Eta());
 759         phiArray.push_back(track->Phi());
 760         chargeArray.push_back(track->Charge());
 761         ++nAcceptedTracks;
 762       }
 763     }
 764     else{
 765       TParticle *partOfTrack = stack->Particle(vtrack->GetLabel());
 766       if (TMath::Abs(partOfTrack->Eta())<fEtaCut && partOfTrack->Pt()>fPtMin && partOfTrack->Pt()<fPtMax) {
 767         ptArray.push_back(partOfTrack->Pt());
 768         etaArray.push_back(partOfTrack->Eta());
 769         phiArray.push_back(partOfTrack->Phi());
 770         chargeArray.push_back(vtrack->Charge());
 771         ++nAcceptedTracks;
 772       }
 773     }
 774
 775     // TPC only track memory needs to be cleaned up
 776     if(track)
 777       delete track;
 778
 779   }
 780
 781   if(nAcceptedTracks==0) return 0;
 782
 783   //tracklet loop
 784   Int_t ntrackletsAccept=0;
 785   AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
 786   Int_t ntracklets = mult->GetNumberOfTracklets();
 787   for(Int_t i=0;i< ntracklets;i++){
 788     if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
 789       ntrackletsAccept++;
 790     }
 791   }
 792
 793   nTracksTracklets.push_back(nAcceptedTracks);
 794   nTracksTracklets.push_back(ntrackletsAccept);
 795   nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
 796                                               //where here also neutral particles are counted.
 797
 798
 799   //get mc vertex for correction maps
 800   AliGenEventHeader*  header = MCEvent()->GenEventHeader();
 801   TArrayF mcV;
 802   header->PrimaryVertex(mcV);
 803   fVzEvent= mcV[2];
 804
 805   return fNRecAccept; // give back reconstructed value
 806
 807
 808 }
 809
 810
 811
 812
 813 //________________________________________________________________________
 814 Int_t AliAnalysisTaskMinijet::LoopESDMC(vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 815                                         vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
 816                                         vector<Int_t> &nTracksTracklets, const Int_t step)
 817 {
 818   // gives back the number of charged prim MC particle and pointer to arrays
 819   // with particle properties (pt, eta, phi)
 820
 821   ptArray.clear();
 822   etaArray.clear();
 823   phiArray.clear();
 824   chargeArray.clear();
 825   nTracksTracklets.clear();
 826
 827   fNMcPrimAccept=0;
 828
 829   AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
 830   if (!mcEvent) {
 831     Error("LoopESDMC", "Could not retrieve MC event");
 832     return 0;
 833   }
 834
 835   AliStack* stack = MCEvent()->Stack();
 836   if(!stack) return 0;
 837
 838   Int_t ntracks = mcEvent->GetNumberOfTracks();
 839   if(fDebug>1)Printf("MC particles: %d", ntracks);
 840
 841   //vertex
 842   AliGenEventHeader*  header = MCEvent()->GenEventHeader();
 843   TArrayF mcV;
 844   header->PrimaryVertex(mcV);
 845   Float_t vzMC = mcV[2];
 846   fVertexZ[step]->Fill(vzMC);
 847
 848   //----------------------------------
 849   //first track loop to check how many chared primary tracks are available
 850   Int_t nChargedPrimaries=0;
 851   Int_t nAllPrimaries=0;
 852
 853   Int_t nPseudoTracklets=0;
 854   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 855     AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
 856     if (!track) {
 857       Error("LoopESDMC", "Could not receive track %d", iTracks);
 858       continue;
 859     }
 860
 861
 862     if(//count also charged particles in case of fSelectParticles==2 (only neutral)
 863        !SelectParticlePlusCharged(
 864                                   track->Charge(),
 865                                   track->PdgCode(),
 866                                   stack->IsPhysicalPrimary(track->Label())
 867                                   )
 868        )
 869       continue;
 870
 871     //count number of pseudo tracklets
 872     if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.0) ++nPseudoTracklets; //0.035
 873     //same cuts as on ESDtracks
 874     if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
 875
 876     //count all primaries
 877     ++nAllPrimaries;
 878     //count charged primaries
 879     if (track->Charge() != 0) ++nChargedPrimaries;
 880
 881     if(fDebug>2) Printf("PDG=%d, IsPrim=%d",  track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
 882
 883
 884   }
 885   //----------------------------------
 886
 887   if(fDebug>2){
 888     Printf("All in acceptance=%d",nAllPrimaries);
 889     Printf("Charged in acceptance =%d",nChargedPrimaries);
 890   }
 891
 892   fChargedPi0->Fill(nAllPrimaries,nChargedPrimaries);
 893
 894   if(nAllPrimaries==0) return 0;
 895   if(nChargedPrimaries==0) return 0;
 896
 897   //track loop
 898   //Int_t iChargedPiK=0;
 899   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 900     AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
 901     if (!track) {
 902       Error("LoopESDMC", "Could not receive track %d", iTracks);
 903       continue;
 904     }
 905
 906     if(!SelectParticle(
 907                        track->Charge(),
 908                        track->PdgCode(),
 909                        stack->IsPhysicalPrimary(track->Label())
 910                        )
 911        )
 912       continue;
 913
 914
 915     //same cuts as on ESDtracks
 916     if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin  || track->Pt()>fPtMax) continue;
 917
 918     if(fDebug>2) Printf("After: PDG=%d, IsPrim=%d",  track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
 919
 920
 921     fHistPtMC->Fill(track->Pt());
 922     //fills arrays with track properties
 923     ptArray.push_back(track->Pt());
 924     etaArray.push_back(track->Eta());
 925     phiArray.push_back(track->Phi());
 926     chargeArray.push_back(track->Charge());
 927
 928
 929   } //track loop
 930
 931   nTracksTracklets.push_back(nChargedPrimaries);
 932   nTracksTracklets.push_back(nPseudoTracklets);
 933   if(fSelectParticles!=2){
 934     nTracksTracklets.push_back(nAllPrimaries);
 935   }
 936   else{
 937     nTracksTracklets.push_back(nAllPrimaries-nChargedPrimaries); // only neutral
 938   }
 939
 940
 941   fNMcPrimAccept=nChargedPrimaries;
 942
 943   if(step==1){
 944     fNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
 945     fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
 946   }
 947
 948   fVzEvent= mcV[2];
 949   return fNRecAccept;
 950
 951 }
 952
 953 //________________________________________________________________________
 954 Int_t AliAnalysisTaskMinijet::LoopAOD( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 955                                        vector<Float_t> &phiArray,  vector<Short_t> &chargeArray,
 956                                        vector<Int_t> &nTracksTracklets, const Int_t step)
 957 {
 958   // gives back the number of AOD tracks and pointer to arrays with track
 959   // properties (pt, eta, phi)
 960
 961   ptArray.clear();
 962   etaArray.clear();
 963   phiArray.clear();
 964   chargeArray.clear();
 965   nTracksTracklets.clear();
 966
 967   TClonesArray *mcArray=0x0;
 968   if(fAnalysePrimOnly){
 969     mcArray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
 970   }
 971
 972
 973   AliAODVertex* vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
 974   Double_t vzAOD=vertex->GetZ();
 975   fVertexZ[step]->Fill(vzAOD);
 976
 977   // Retreive the number of tracks for this event.
 978   Int_t ntracks = fAODEvent->GetNumberOfTracks();
 979   if(fDebug>1) Printf("AOD tracks: %d", ntracks);
 980
 981
 982   Int_t nAcceptedTracks=0;
 983   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 984     AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
 985     if (!track) {
 986       Error("LoopAOD", "Could not receive track %d", iTracks);
 987       continue;
 988     }
 989
 990     AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
 991
 992     //use only tracks from primaries
 993     if(fAnalysePrimOnly){
 994       if(vtrack->GetLabel()<0)continue;
 995       if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
 996     }
 997
 998     if(track->TestFilterBit(128) && TMath::Abs(track->Eta())<fEtaCut
 999        && track->Pt()>fPtMin && track->Pt()<fPtMax){
1000
1001       nAcceptedTracks++;
1002
1003       // Printf("dca= %f", track->DCA());
1004       //save track properties in vector
1005       ptArray.push_back(track->Pt());
1006       etaArray.push_back(track->Eta());
1007       phiArray.push_back(track->Phi());
1008       chargeArray.push_back(track->Charge());
1009       fHistPt->Fill(track->Pt());
1010
1011     }
1012   }
1013   //need to check this option for MC
1014   if(nAcceptedTracks==0) return 0;
1015
1016
1017   //tracklet loop
1018   Int_t ntrackletsAccept=0;
1019   AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1020   for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1021     if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut){
1022       ++ntrackletsAccept;
1023     }
1024   }
1025
1026
1027   nTracksTracklets.push_back(nAcceptedTracks);
1028   nTracksTracklets.push_back(ntrackletsAccept);
1029   nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1030                                               //where here also neutral particles are counted.
1031
1032
1033   fVzEvent= vzAOD;
1034   if(step==5 || step==2)fNRecAccept = nAcceptedTracks; // needed for MC case //step5 = TrigVtxRecNrec
1035   return fNRecAccept; // at the moment, always return reconstructed multiplicity
1036
1037 }
1038
1039 //________________________________________________________________________
1040 Int_t AliAnalysisTaskMinijet::LoopAODRecMcProp( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1041                                                 vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1042                                                 vector<Int_t> &nTracksTracklets, const Int_t step)
1043 {
1044   // gives back the number of AOD tracks and pointer to arrays with track
1045   // properties (pt, eta, phi)
1046
1047
1048   ptArray.clear();
1049   etaArray.clear();
1050   phiArray.clear();
1051   chargeArray.clear();
1052   nTracksTracklets.clear();
1053
1054
1055   // Retreive the number of tracks for this event.
1056   Int_t ntracks = fAODEvent->GetNumberOfTracks();
1057   if(fDebug>1) Printf("AOD tracks: %d", ntracks);
1058
1059
1060   //get array of mc particles
1061   TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1062     FindListObject(AliAODMCParticle::StdBranchName());
1063   if(!mcArray){
1064     Printf("No MC particle branch found");
1065     return kFALSE;
1066   }
1067
1068   AliAODVertex* vtx= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1069   Double_t vzAOD=vtx->GetZ();
1070   fVertexZ[step]->Fill(vzAOD);
1071
1072   Int_t nAcceptedTracks=0;
1073   for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1074     AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
1075
1076     AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
1077
1078     if (!track) {
1079       Error("LoopAODRecMcProp", "Could not receive track %d", iTracks);
1080       continue;
1081     }
1082
1083     //use only tracks from primaries
1084     if(fAnalysePrimOnly){
1085       if(vtrack->GetLabel()<0)continue;
1086       if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
1087     }
1088
1089     if(track->TestFilterBit(128) &&  TMath::Abs(track->Eta())<fEtaCut &&
1090        track->Pt()>fPtMin && track->Pt()<fPtMax){
1091
1092       nAcceptedTracks++;
1093
1094       //save track properties in vector
1095       if(vtrack->GetLabel()<0){ //fake tracks
1096         //      Printf("Fake track");
1097         //      continue;
1098         ptArray.push_back(track->Pt());
1099         etaArray.push_back(track->Eta());
1100         phiArray.push_back(track->Phi());
1101         chargeArray.push_back(track->Charge());
1102       }
1103       else{//mc properties
1104         AliAODMCParticle *partOfTrack = (AliAODMCParticle*)mcArray->At(vtrack->GetLabel());
1105
1106         ptArray.push_back(partOfTrack->Pt());
1107         etaArray.push_back(partOfTrack->Eta());
1108         phiArray.push_back(partOfTrack->Phi());
1109         chargeArray.push_back(vtrack->Charge());//partOfTrack?
1110       }
1111
1112     }
1113   }
1114   //need to check this option for MC
1115   if(nAcceptedTracks==0) return 0;
1116
1117   //tracklet loop
1118   Int_t ntrackletsAccept=0;
1119   AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1120   for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1121     if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut ){
1122       ++ntrackletsAccept;
1123     }
1124   }
1125
1126
1127   nTracksTracklets.push_back(nAcceptedTracks);
1128   nTracksTracklets.push_back(ntrackletsAccept);
1129   nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1130                                               //where here also neutral particles are counted.
1131
1132
1133   //check vertex (mc)
1134   AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1135     FindListObject(AliAODMCHeader::StdBranchName());
1136   Float_t vzMC = aodMCheader->GetVtxZ();
1137
1138   fVzEvent= vzMC;
1139   return fNRecAccept;//this is the rec value from step 5
1140
1141 }
1142
1143
1144
1145 //________________________________________________________________________
1146 Int_t AliAnalysisTaskMinijet::LoopAODMC( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1147                                          vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1148                                          vector<Int_t> &nTracksTracklets, const Int_t step)
1149 {
1150   // gives back the number of AOD MC particles and pointer to arrays with particle
1151   // properties (pt, eta, phi)
1152
1153   ptArray.clear();
1154   etaArray.clear();
1155   phiArray.clear();
1156   chargeArray.clear();
1157   nTracksTracklets.clear();
1158
1159   //check vertex
1160   AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1161     FindListObject(AliAODMCHeader::StdBranchName());
1162   Float_t vzMC = aodMCheader->GetVtxZ();
1163   fVertexZ[step]->Fill(vzMC);
1164
1165
1166   //retreive MC particles from event
1167   TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1168     FindListObject(AliAODMCParticle::StdBranchName());
1169   if(!mcArray){
1170     Printf("No MC particle branch found");
1171     return kFALSE;
1172   }
1173
1174   Int_t ntracks = mcArray->GetEntriesFast();
1175   if(fDebug>1)Printf("MC particles: %d", ntracks);
1176
1177
1178   // Track loop: chek how many particles will be accepted
1179   //Float_t vzMC=0.;
1180   Int_t nChargedPrim=0;
1181   Int_t nAllPrim=0;
1182   Int_t nPseudoTracklets=0;
1183   for (Int_t it = 0; it < ntracks; it++) {
1184     AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1185     if (!track) {
1186       Error("LoopAODMC", "Could not receive particle %d", it);
1187       continue;
1188     }
1189
1190     if(!SelectParticlePlusCharged(
1191                                   track->Charge(),
1192                                   track->PdgCode(),
1193                                   track->IsPhysicalPrimary()
1194                                   )
1195        )
1196       continue;
1197
1198     if(TMath::Abs(track->Eta())<fEtaCut && track->Pt()>0.0)++nPseudoTracklets; //0.035
1199     if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1200
1201     nAllPrim++;
1202     if(track->Charge()!=0) nChargedPrim++;
1203
1204   }
1205
1206
1207   if(nAllPrim==0) return 0;
1208   if(nChargedPrim==0) return 0;
1209
1210   //generate array with size of number of accepted tracks
1211   fChargedPi0->Fill(nAllPrim,nChargedPrim);
1212
1213
1214   // Track loop: fill arrays for accepted tracks
1215   // Int_t iChargedPrim=0;
1216   for (Int_t it = 0; it < ntracks; it++) {
1217     AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1218     if (!track) {
1219       Error("LoopAODMC", "Could not receive particle %d", it);
1220       continue;
1221     }
1222
1223     if(!SelectParticle(
1224                        track->Charge(),
1225                        track->PdgCode(),
1226                        track->IsPhysicalPrimary()
1227                        )
1228        )
1229       continue;
1230
1231     if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1232
1233     fHistPtMC->Fill(track->Pt());
1234     ptArray.push_back(track->Pt());
1235     etaArray.push_back(track->Eta());
1236     phiArray.push_back(track->Phi());
1237     chargeArray.push_back(track->Charge());
1238   }
1239
1240   nTracksTracklets.push_back(nChargedPrim);
1241   nTracksTracklets.push_back(nPseudoTracklets);
1242   if(fSelectParticles!=2){
1243     nTracksTracklets.push_back(nAllPrim);
1244   }
1245   else{
1246     nTracksTracklets.push_back(nAllPrim-nChargedPrim); // only neutral
1247   }
1248
1249
1250
1251   fVzEvent= vzMC;
1252   fNMcPrimAccept = nChargedPrim;
1253   if(step==1){ // step 1 = Trig All Mc Nmc
1254     fNmcNch->Fill( fNMcPrimAccept,fNRecAccept);
1255     fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1256   }
1257   return fNRecAccept; // rec value from step 5 or step 2
1258
1259
1260 }
1261
1262 //________________________________________________________________________
1263 void AliAnalysisTaskMinijet::Analyse(const vector<Float_t> &pt,
1264                                      const vector<Float_t> &eta,
1265                                      const vector<Float_t> &phi,
1266                                      const vector<Short_t> &charge,
1267                                      const Int_t ntracksCharged,
1268                                      const Int_t ntracklets,
1269                                      const Int_t nAll,
1270                                      const Int_t step)
1271 {
1272
1273   // analyse track properties (comming from either ESDs or AODs) in order to compute
1274   // mini jet activity (chared tracks) as function of charged multiplicity
1275
1276   fStep->Fill(step);
1277
1278   if(fDebug){
1279     Printf("Analysis Step=%d", step);
1280     if(fDebug>2){
1281       Printf("nAll=%d",nAll);
1282       Printf("nCharged=%d",ntracksCharged);
1283       for (Int_t i = 0; i < nAll; i++) {
1284         Printf("pt[%d]=%f",i,pt[i]);
1285       }
1286     }
1287   }
1288
1289   //calculation of mean pt for all tracks in case of step==0
1290   if(step==5 || step ==2){ // rec step
1291     Double_t meanPt=0.;
1292     Double_t leadingPt=0.;
1293     for (Int_t i = 0; i < nAll; i++) {
1294       if(pt[i]<0.01)continue;
1295       meanPt+=pt[i];
1296       if(leadingPt<pt[i])leadingPt=pt[i];
1297     }
1298     meanPt=meanPt/nAll;
1299     fMeanPtRec=meanPt;
1300     fLeadingPtRec=leadingPt;
1301   }
1302
1303   Double_t propEvent[4] = {ntracksCharged,fVzEvent,fMeanPtRec, fLeadingPtRec}; //vz: {rec, mc, mc}, meanPt and Nrec is always rec value
1304   fMapEvent[step]->Fill(propEvent);
1305
1306   fNcharge[step]->Fill(ntracksCharged);
1307
1308   Float_t ptEventAxis=0;  // pt event axis
1309   Float_t etaEventAxis=0; // eta event axis
1310   Float_t phiEventAxis=0; // phi event axis
1311   Short_t chargeEventAxis=0; // charge event axis
1312
1313   Float_t ptOthers  = 0; // pt others // for all other tracks around event axis -> see loop
1314   Float_t etaOthers = 0; // eta others
1315   Float_t phiOthers = 0; // phi others
1316   Short_t chargeOthers = 0; // charge others
1317
1318   Int_t   *pindexInnerEta  = new Int_t  [nAll+1];
1319   Float_t *ptInnerEta      = new Float_t[nAll+1];
1320
1321
1322
1323   for (Int_t i = 0; i < nAll; i++) {
1324
1325     if(pt[i]<0.01)continue;
1326
1327     //fill single particle correction for first step of pair correction
1328     Double_t propAll[3] = {eta[i],pt[i],ntracksCharged };
1329     fMapAll[step]->Fill(propAll);
1330
1331
1332     //filling of simple check plots
1333     if(pt[i]<0.7) continue;
1334     fPt[step]    -> Fill( pt[i]);
1335     fEta[step]   -> Fill(eta[i]);
1336     fPhi[step]   -> Fill(phi[i]);
1337     fPhiEta[step]-> Fill(phi[i], eta[i]);
1338
1339     pindexInnerEta[i]=0; //set all values to zero
1340     //fill new array for eta check
1341     ptInnerEta[i]= pt[i];
1342
1343   }
1344
1345
1346
1347   // define event axis: leading or random track (with pt>fTriggerPtCut)
1348   // ---------------------------------------
1349   Int_t highPtTracks=0;
1350   Int_t highPtTracksInnerEta=0;
1351   // Double_t highPtTracksInnerEtaCorr=0;
1352   Int_t mult09=0;
1353
1354   //count high pt tracks and high pt tracks in acceptance for seeds
1355   for(Int_t i=0;i<nAll;i++){
1356
1357     if(pt[i]<0.01)continue;
1358
1359     if(TMath::Abs(eta[i])<0.9){
1360       mult09++;
1361     }
1362
1363     if(pt[i]>fTriggerPtCut) {
1364       highPtTracks++;
1365     }
1366
1367     // seed should be place in middle of acceptance, that complete cone is in acceptance
1368     if(pt[i]>fTriggerPtCut && TMath::Abs(eta[i])<fEtaCutSeed && charge[i]!=0){
1369
1370       // Printf("eta=%f", eta[i]);
1371       highPtTracksInnerEta++;
1372
1373     }
1374     else{
1375       ptInnerEta[i]=0;
1376     }
1377   }
1378
1379
1380   //sort array in order to get highest pt tracks first
1381   //index can be computed with pindexInnerEta[number]
1382   if(nAll) TMath::Sort(nAll, ptInnerEta, pindexInnerEta, kTRUE);
1383
1384   //     plot of multiplicity distributions
1385   fNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1386   fPNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1387
1388   if(ntracklets){
1389     fNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1390     fNchTracklet[step]->Fill(ntracklets, ntracksCharged);
1391     fPNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1392   }
1393
1394   //analysis can only be performed with event axis, defined by high pt track
1395
1396
1397   if(highPtTracks>0 && highPtTracksInnerEta>0){
1398
1399     // build pairs in two track loops
1400     // loop over all possible trigger particles (defined by pt_trig and eta_acceptance)
1401     for(Int_t axis=0;(axis<nAll) && (pt[pindexInnerEta[axis]]>=fTriggerPtCut); axis++){
1402
1403       //EventAxisRandom track properties
1404       ptEventAxis  = pt [pindexInnerEta[axis]];
1405       etaEventAxis = eta[pindexInnerEta[axis]];
1406       phiEventAxis = phi[pindexInnerEta[axis]];
1407       chargeEventAxis = charge[pindexInnerEta[axis]];
1408       fPtSeed[step]    -> Fill( ptEventAxis);
1409       fEtaSeed[step]   -> Fill(etaEventAxis);
1410       fPhiSeed[step]   -> Fill(phiEventAxis);
1411
1412
1413       Double_t prop[3] = {etaEventAxis,ptEventAxis,ntracksCharged };
1414       fMapSingleTrig[step]->Fill(prop);
1415
1416       //associated tracks
1417       for (Int_t iTrack = axis+1; iTrack < nAll; iTrack++) {
1418
1419         if(pt[pindexInnerEta[iTrack]]<fAssociatePtCut) continue;
1420
1421         if(fSelectParticlesAssoc==1){
1422           if(charge[pindexInnerEta[iTrack]]==0)continue;
1423         }
1424         if(fSelectParticlesAssoc==2){
1425           if(charge[pindexInnerEta[iTrack]]!=0)continue;
1426         }
1427
1428
1429         ptOthers   = pt [pindexInnerEta[iTrack]];
1430         etaOthers  = eta[pindexInnerEta[iTrack]];
1431         phiOthers  = phi[pindexInnerEta[iTrack]];
1432         chargeOthers = charge[pindexInnerEta[iTrack]];
1433
1434
1435         //plot only properties of tracks with pt>ptassoc
1436         fPtOthers[step]    -> Fill( ptOthers);
1437         fEtaOthers[step]   -> Fill(etaOthers);
1438         fPhiOthers[step]   -> Fill(phiOthers);
1439         fPtEtaOthers[step]   -> Fill(ptOthers, etaOthers);
1440
1441         //      if(fDebug>2)Printf("%f, %f", pt[pindexInnerEta[axis]], pt[pindexInnerEta[iTrack]]);
1442
1443         Float_t dPhi = phiOthers-phiEventAxis;
1444         if(dPhi>1.5*TMath::Pi()) dPhi = dPhi-2*TMath::Pi();
1445         else if(dPhi<-0.5*TMath::Pi())dPhi=dPhi+2*TMath::Pi();
1446         Float_t dEta=etaOthers-etaEventAxis;
1447
1448
1449         fDPhiDEtaEventAxis[step]->Fill(dPhi, dEta);
1450         fDPhiEventAxis[step]->Fill(dPhi);
1451
1452         //check outliers
1453         if(ptEventAxis< 0.4 || ptEventAxis > 100) Printf("particles out of range pt");
1454         if(ntracksCharged<0 || ntracksCharged>150) Printf("particles out of range ncharge");
1455         if(TMath::Abs(dEta)>2*fEtaCut) {
1456           Printf("particles out of range dEta");
1457           Printf("eta1=%f, eta2=%f", etaOthers, etaEventAxis);
1458           Printf("step=%d",step);
1459         }
1460         if(dPhi<-0.5*TMath::Pi() || dPhi>1.5*TMath::Pi()){
1461           Printf("particles out of range dPhi");
1462           Printf("phi1=%f, phi2=%f", phiOthers, phiEventAxis);
1463         }
1464
1465         Bool_t isLikeSign = CheckLikeSign(chargeEventAxis, chargeOthers);
1466
1467         Double_t prop6[6] = {ptEventAxis,ptOthers,dEta,dPhi,ntracksCharged, isLikeSign };
1468         fMapPair[step]->Fill(prop6);
1469
1470       }// end of inner track loop
1471
1472     } //end of outer track loop
1473
1474     fTriggerNch[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1475     fTriggerNchSeeds[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1476     fTriggerTracklet[step]->Fill(ntracklets);
1477
1478
1479   }//if there is at least one high pt track
1480
1481
1482   if(pindexInnerEta){// clean up array memory used for TMath::Sort
1483     delete[] pindexInnerEta;
1484     pindexInnerEta=0;
1485   }
1486
1487   if(ptInnerEta){// clean up array memory used for TMath::Sort
1488     delete[] ptInnerEta;
1489     ptInnerEta=0;
1490   }
1491
1492 }
1493
1494
1495
1496 //________________________________________________________________________
1497 void AliAnalysisTaskMinijet::Terminate(Option_t*)
1498 {
1499   //terminate function is called at the end
1500   //can be used to draw histograms etc.
1501
1502
1503 }
1504
1505 //________________________________________________________________________
1506 Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(const Short_t charge, const Int_t pdg, Bool_t prim)
1507 {
1508   //selection of mc particle
1509   //fSelectParticles=0: use charged primaries and pi0 and k0
1510   //fSelectParticles=1: use only charged primaries
1511   //fSelectParticles=2: use only pi0 and k0
1512
1513   if(fSelectParticles==0 || fSelectParticles==2){ // in case of 2: need to count also charged particles
1514     if(charge==0){
1515       if(!(pdg==111||pdg==130||pdg==310))
1516         return false;
1517     }
1518     else{// charge !=0
1519       if(!prim)
1520         return false;
1521     }
1522   }
1523
1524   else if(fSelectParticles==1){
1525     if (charge==0 || !prim){
1526       return false;
1527     }
1528   }
1529
1530   else{
1531     Printf("Error: wrong selection of charged/pi0/k0");
1532     return 0;
1533   }
1534
1535   return true;
1536
1537 }
1538
1539 //________________________________________________________________________
1540 Bool_t AliAnalysisTaskMinijet::SelectParticle(const Short_t charge, const Int_t pdg, const Bool_t prim)
1541 {
1542   //selection of mc particle
1543   //fSelectParticles=0: use charged primaries and pi0 and k0
1544   //fSelectParticles=1: use only charged primaries
1545   //fSelectParticles=2: use only pi0 and k0
1546
1547   if(fSelectParticles==0){
1548     if(charge==0){
1549       if(!(pdg==111||pdg==130||pdg==310))
1550         return false;
1551     }
1552     else{// charge !=0
1553       if(!prim)
1554         return false;
1555     }
1556   }
1557   else if (fSelectParticles==1){
1558     if (charge==0 || !prim){
1559       return false;
1560     }
1561   }
1562   else if(fSelectParticles==2){
1563     if(!(pdg==111||pdg==130||pdg==310))
1564       return false;
1565   }
1566
1567   return true;
1568
1569 }
1570
1571 //________________________________________________________________________
1572 Bool_t AliAnalysisTaskMinijet::CheckEvent(const Bool_t recVertex)
1573 {
1574   // This function tests the quality of an event (ESD/AOD) (rec and/or mc part)
1575   // recVertex=false:  check if Mc events and stack is there, Nmc>0
1576   // recVertex=false: " + check if there is a good, reconstructed SPD vertex
1577   // defined by |z|<fVertexCut(10cm), Contributer>0, no PileUpFromSPD(3,0,8)
1578
1579
1580   if(fMode==0){//esd
1581
1582     //mc
1583     if(fUseMC){
1584
1585       //mc event
1586       AliMCEvent *mcEvente = (AliMCEvent*) MCEvent();
1587       if (!mcEvente) {
1588         Error("CheckEvent", "Could not retrieve MC event");
1589         return false;
1590       }
1591
1592       //stack
1593       AliStack* stackg = MCEvent()->Stack();
1594       if(!stackg) return false;
1595       Int_t ntracksg = mcEvente->GetNumberOfTracks();
1596       if(ntracksg<0) return false;
1597
1598       //vertex
1599       //AliGenEventHeader*  headerg = MCEvent()->GenEventHeader();
1600       //TArrayF mcVg;
1601       //headerg->PrimaryVertex(mcVg);
1602       //  if(TMath::Abs(mcVg[0])<1e-8 && TMath::Abs(mcVg[0])<1e-8 &&
1603       // TMath::Abs(mcVg[0])<1e-8) return false;
1604       // Float_t vzMCg = mcVg[2];
1605       // if(TMath::Abs(vzMCg)>fVertexZCut) return false;
1606     }
1607
1608     //rec
1609     if(recVertex==true){
1610       if(fESDEvent->IsPileupFromSPD(3,0,8))return false;
1611
1612       //rec vertex
1613       // const AliESDVertex*    vertexESDg   = fESDEvent->GetPrimaryVertex(); // uses track or SPD vertexer
1614       // if(!vertexESDg) return false;
1615       // if(vertexESDg->GetNContributors()<=0)return false;
1616       // Float_t fVzg= vertexESDg->GetZ();
1617       // if(TMath::Abs(fVzg)>fVertexZCut) return false;
1618
1619       //rec spd vertex
1620       const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
1621       if(!vtxSPD) return false;
1622       if(vtxSPD->GetNContributors()<=0)return false;
1623       Float_t fVzSPD= vtxSPD->GetZ();
1624       if(TMath::Abs(fVzSPD)>fVertexZCut) return false;
1625     }
1626     return true;
1627
1628   }
1629
1630
1631   else if(fMode==1){ //aod
1632
1633     if(fUseMC){
1634       //mc
1635       //   AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1636       //         FindListObject(AliAODMCHeader::StdBranchName());
1637       //        Float_t vzMC = aodMCheader->GetVtxZ();
1638       //        if(TMath::Abs(vzMC)>fVertexZCut) return false;
1639
1640       //retreive MC particles from event
1641       TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1642         FindListObject(AliAODMCParticle::StdBranchName());
1643       if(!mcArray){
1644         Printf("No MC particle branch found");
1645         return false;
1646       }
1647     }
1648
1649     //rec
1650     if(recVertex==true){
1651       if(fAODEvent->IsPileupFromSPD(3,0.8))return false;
1652
1653       //      AliAODVertex*     vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertex();
1654       //      if(!vertex) return false;
1655       //      if(vertex->GetNContributors()<=0) return false;
1656       //      Double_t vzAOD=vertex->GetZ();
1657       //      if(TMath::Abs(vzAOD)<1e-9) return false;
1658       //      if(TMath::Abs(vzAOD)>fVertexZCut) return false;
1659
1660       AliAODVertex*     vertexSPD= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();
1661       if(!vertexSPD) return false;
1662       if(vertexSPD->GetNContributors()<=0) return false;
1663       Double_t vzSPD=vertexSPD->GetZ();
1664       //if(TMath::Abs(vzSPD)<1e-9) return false;
1665       if(TMath::Abs(vzSPD)>fVertexZCut) return false;
1666     }
1667     return true;
1668
1669   }
1670
1671   else {
1672     Printf("Wrong mode!");
1673     return false;
1674   }
1675
1676 }
1677
1678 //_____________________________________________________________________________
1679 const Double_t * AliAnalysisTaskMinijet::CreateLogAxis(const Int_t nbins,
1680                                                        const Double_t xmin,
1681                                                        const Double_t xmax)
1682 {
1683   // returns pointer to an array with can be used to build a logarithmic axis
1684   // it is user responsibility to delete the array
1685
1686   Double_t logxmin = TMath::Log10(xmin);
1687   Double_t logxmax = TMath::Log10(xmax);
1688   Double_t binwidth = (logxmax-logxmin)/nbins;
1689
1690   Double_t *xbins =  new Double_t[nbins+1];
1691
1692   xbins[0] = xmin;
1693   for (Int_t i=1;i<=nbins;i++) {
1694     xbins[i] = xmin + TMath::Power(10,logxmin+i*binwidth);
1695   }
1696
1697   return xbins;
1698 }
1699
1700 //_____________________________________________________________________________
1701 Bool_t AliAnalysisTaskMinijet::CheckLikeSign(const Short_t chargeEventAxis,
1702                                              const Short_t chargeOthers)
1703 {
1704   // compute if charge of two particles/tracks has same sign or different sign
1705
1706   if(fDebug>2) Printf("Charge1=%d, Charge2=%d",chargeEventAxis,chargeOthers);
1707
1708   if(chargeEventAxis<0){
1709     if(chargeOthers<0){
1710       return true;
1711     }
1712     else if(chargeOthers>0){
1713       return false;
1714     }
1715   }
1716
1717   else if(chargeEventAxis>0){
1718     if(chargeOthers>0){
1719       return true;
1720     }
1721     else if(chargeOthers<0){
1722       return false;
1723     }
1724   }
1725
1726   else{
1727     Printf("Error: Charge not lower nor higher as zero");
1728     return false;
1729   }
1730
1731   Printf("Error: Check values of Charge");
1732   return false;
1733 }
1734