PWGCF/Correlations/DPhi/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 #include "AliCentrality.h"
  30
  31 #include "AliAnalysisManager.h"
  32 #include "AliInputEventHandler.h"
  33
  34 #include <vector>
  35 #include <algorithm>
  36 #include <functional>
  37 using namespace std;
  38
  39 #include "AliAnalysisTaskMinijet.h"
  40
  41 // Analysis task for two-particle correlations using all particles over pt threshold
  42 // pt_trig threshold for trigger particle (event axis) and pt_assoc for possible associated particles.
  43 // Extract mini-jet yield and fragmentation properties via Delta-Phi histograms of these correlations
  44 // post processing of analysis output via macro plot3and2Gaus.C
  45 // Can use ESD or AOD, reconstructed and Monte Carlo data as input
  46 // Author: Eva Sicking, modifications by Emilia Leogrande
  47
  48
  49 ClassImp(AliAnalysisTaskMinijet)
  50
  51 //________________________________________________________________________
  52 AliAnalysisTaskMinijet::AliAnalysisTaskMinijet(const char *name)
  53     : AliAnalysisTaskSE(name),
  54     fUseMC(kFALSE),
  55     fMcOnly(kFALSE),
  56     fBSign(0),
  57     fAnalysePrimOnly(kFALSE),// not used
  58     fPtMin(0.2),
  59     fPtMax(50.0),
  60     fCuts(0),
  61     fTriggerPtCut(0.7),
  62     fAssociatePtCut(0.4),
  63     fMode(0),
  64     fTriggerType(1),
  65     fFilterBit(128),
  66     fVertexZCut(10.),
  67     fEtaCut(0.9),
  68     fEtaCutSeed(0.9),
  69     fSelectParticles(1),
  70     fSelectParticlesAssoc(1),
  71     fCheckSDD(true),
  72     fSelOption(1),
  73     fCorrStrangeness(true),
  74     fThreeParticleCorr(false),
  75     fRejectChunks(false),
  76     fNTPC(10),
  77     fESDEvent(0),
  78     fAODEvent(0),
  79     fNMcPrimAccept(0),
  80     fNRecAccept(0),
  81     fNRecAcceptStrangeCorr(0),
  82     fNMcPrimAcceptTracklet(0),
  83     fNRecAcceptTracklet(0),
  84     fVzEvent(0),
  85     fMeanPtRec(0),
  86     fLeadingPtRec(0),
  87     fHists(0),
  88     fStep(0),
  89     fEventStat(0),
  90     fHistPt(0),
  91     fHistPtMC(0),
  92     fNContrNtracklets(0),
  93     fNContrNtracks(0),
  94     fCorruptedChunks(0),
  95     fCorruptedChunksAfter(0),
  96     fNmcNch(0),
  97     fPNmcNch(0),
  98     fNmcNchVtx(0),
  99     fNmcNchVtxStrangeCorr(0),
 100     fPNmcNchVtx(0),
 101     fNmcNchTracklet(0),
 102     fPNmcNchTracklet(0),
 103     fNmcNchVtxTracklet(0),
 104     fPNmcNchVtxTracklet(0),
 105     fChargedPi0(0),
 106     fVertexCheck(0),
 107     fPropagateDca(0),
 108     fCentralityMethod("")
 109 {
 110
 111     //Constructor
 112
 113     for(Int_t i = 0;i< 8;i++){
 114         fMapSingleTrig[i]         =  0;
 115         fMapPair[i]               =  0;
 116         fMapEvent[i]              =  0;
 117         fMapAll[i]                =  0;
 118         fMapThree[i]              =  0;
 119
 120         fVertexZ[i]               =  0;
 121
 122         fNcharge[i]               =  0;
 123         fPt[i]                    =  0;
 124         fEta[i]                   =  0;
 125         fPhi[i]                   =  0;
 126         fDcaXY[i]                 =  0;
 127         fDcaZ[i]                  =  0;
 128
 129         fPtSeed[i]                =  0;
 130         fEtaSeed[i]               =  0;
 131         fPhiSeed[i]               =  0;
 132
 133         fPtOthers[i]              =  0;
 134         fEtaOthers[i]             =  0;
 135         fPhiOthers[i]             =  0;
 136         fPtEtaOthers[i]           =  0;
 137
 138         fPhiEta[i]                =  0;
 139
 140         fDPhiDEtaEventAxis[i]     =  0;
 141         fDPhiDEtaEventAxisSeeds[i]=  0;
 142         fTriggerNch[i]            =  0;
 143
 144         fTriggerNchSeeds[i]       =  0;
 145         fTriggerTracklet[i]       =  0;
 146
 147         fNch07Nch[i]              =  0;
 148         fPNch07Nch[i]             =  0;
 149
 150         fNch07Tracklet[i]         =  0;
 151         fNchTracklet[i]           =  0;
 152         fPNch07Tracklet[i]        =  0;
 153         fDPhiEventAxis[i]         =  0;
 154         fDPhi1DPhi2[i]            =  0;
 155     }
 156     DefineOutput(1, TList::Class());
 157 }
 158
 159 //________________________________________________________________________
 160 AliAnalysisTaskMinijet::~AliAnalysisTaskMinijet()
 161 {
 162     // Destructor
 163
 164     if (fHists && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHists;
 165 }
 166
 167 //________________________________________________________________________
 168 void AliAnalysisTaskMinijet::UserCreateOutputObjects()
 169 {
 170     // Create histograms
 171     // Called once
 172     if(fDebug) Printf("In User Create Output Objects.");
 173
 174     Int_t nbinsCentr = 0;
 175     Float_t minbinCentr=0, maxbinCentr=0;
 176
 177     if (fCentralityMethod.Length() > 0)
 178     {
 179         nbinsCentr = 100;
 180         minbinCentr=0;
 181         maxbinCentr=100;
 182     }
 183     else
 184     {
 185         nbinsCentr = 101;
 186         minbinCentr=-0.5;
 187         maxbinCentr=100.5;
 188     }
 189
 190     fStep = new TH1F("fStep", "fStep", 10, -0.5, 9.5);
 191     fEventStat = new TH1F("fEventStat", "fEventStat", 10, -0.5, 9.5);
 192     fHistPt = new TH1F("fHistPt", "P_{T} distribution REC", 150, 0.1, 3.1);
 193     fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 194     fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 195     fHistPt->SetMarkerStyle(kFullCircle);
 196     fNContrNtracklets = new TH2F ("fNContrNtracklets", ";N_{tracklets};N_{vtx contrib}", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
 197     fNContrNtracks    = new TH2F ("fNContrNtracks", ";N_{tracks};N_{vtx contrib}", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
 198     fCorruptedChunks  = new TH2F ("fCorruptedChunks",
 199                                   ";N_{tracks,TPC};N_{tracks,ITS-TPC}", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
 200     fCorruptedChunksAfter  = new TH2F ("fCorruptedChunksAfter",
 201                                        ";N_{tracks,TPC};N_{tracks,ITS-TPC}", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
 202
 203
 204
 205
 206     if (fUseMC) {
 207         fHistPtMC = new TH1F("fHistPtMC", "P_{T} distribution MC", 150, 0.1, 3.1);
 208         fHistPtMC->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 209         fHistPtMC->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 210         fHistPtMC->SetMarkerStyle(kFullCircle);
 211
 212         fNmcNch = new TH2F("fNmcNch", "fNmcNch", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
 213         fPNmcNch = new TProfile("pNmcNch", "pNmcNch", nbinsCentr,minbinCentr,maxbinCentr);
 214         fNmcNchVtx = new TH2F("fNmcNchVtx", "fNmcNchVtx", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
 215         fNmcNchVtxStrangeCorr = new TH2F("fNmcNchVtxStrangeCorr", "fNmcNchVtxStrangeCorr", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
 216         fPNmcNchVtx = new TProfile("pNmcNchVtx", "pNmcNchVtx", nbinsCentr,minbinCentr,maxbinCentr);
 217
 218         fNmcNchTracklet = new TH2F("fNmcNchTracklet", "fNmcNchTracklet", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
 219         fPNmcNchTracklet = new TProfile("pNmcNchTracklet", "pNmcNchTracklet", nbinsCentr,minbinCentr,maxbinCentr);
 220         fNmcNchVtxTracklet = new TH2F("fNmcNchVtxTracklet", "fNmcNchVtxTracklet", nbinsCentr,minbinCentr,maxbinCentr,nbinsCentr,minbinCentr,maxbinCentr);
 221         fPNmcNchVtxTracklet = new TProfile("pNmcNchVtxTracklet", "pNmcNchVtxTracklet", nbinsCentr,minbinCentr,maxbinCentr);
 222
 223
 224
 225     }
 226
 227     fChargedPi0  = new TH2F("fChargedPi0", "fChargedPi0", 200, -0.5, 199.5, 200, -0.5, 199.5);
 228     fVertexCheck = new TH1F("fVertexCheck", "fVertexCheck", 100, -0.5, 99.5);
 229     fPropagateDca = new TH1F("fPropagateDca", "fPropagateDca", 2, -0.5, 1.5);
 230
 231     //----------------------
 232     //bins for pt in THnSpare
 233     Double_t ptMin = 0.0, ptMax = 100.;
 234     Int_t nPtBins = 37;
 235     Double_t binsPt[]  = {0.0, 0.1, 0.2, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8,
 236         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,
 237         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};
 238
 239     //Int_t nPtBinsAll = 100;
 240     //Double_t ptMinAll = 1.e-2, ptMaxAll = 100.;
 241     //Double_t *binsPtAll = 0;
 242     //binsPtAll = (Double_t*)CreateLogAxis(nPtBinsAll,ptMinAll,ptMaxAll);
 243
 244
 245     //   Double_t ptMin2 = 0.0, ptMax2 = 100.;
 246     //   Int_t nPtBins2 = 9;
 247     //   Double_t binsPt2[]  = {0.1, 0.4, 0.7, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 100.0};
 248
 249     //  Int_t nPtBins2 = 10;
 250     //   Double_t ptMin2 = 0.4, ptMax2 = 100.;
 251     //   Double_t *binsPt2 = 0;
 252     //   binsPt2 = CreateLogAxis(nPtBins2,ptMin2,ptMax2);
 253
 254     //3 dim matrix
 255     Int_t binsEffHisto[3]   = {Int_t(fEtaCut*20),  nPtBins,       nbinsCentr };
 256     Double_t minEffHisto[3] = {-fEtaCut,           ptMin,        minbinCentr };
 257     Double_t maxEffHisto[3] = {fEtaCut,            ptMax,        maxbinCentr};
 258
 259     //5 dim matrix
 260     Int_t binsEffHisto5[6]   = {  nPtBins,   nPtBins,    1,                              90,         nbinsCentr ,      2 };
 261     Double_t minEffHisto5[6] = {  ptMin,     ptMin,     -2*fEtaCut,          -0.5*TMath::Pi(),       minbinCentr ,   -0.5 };
 262     Double_t maxEffHisto5[6] = {  ptMax,     ptMax,      2*fEtaCut,           1.5*TMath::Pi(),      maxbinCentr ,    1.5 };
 263
 264
 265     //4 dim matrix
 266     Int_t binsEvent[4]   = {   nbinsCentr,        20,   50,  nPtBins };
 267     Double_t minEvent[4] = {   minbinCentr,       -10,    0,    ptMin };
 268     Double_t maxEvent[4] = {   maxbinCentr,        10,   10,    ptMax };
 269
 270     //3 dim matrix
 271     Int_t binsAll[3]   = {Int_t(fEtaCut*20),  nPtBins,       nbinsCentr };
 272     Double_t minAll[3] = {-fEtaCut,           ptMin,         minbinCentr };
 273     Double_t maxAll[3] = {fEtaCut,            ptMax,        maxbinCentr };
 274
 275     //3 dim matrix
 276     Int_t binsThree[3]   = {              90,               90,    nbinsCentr};
 277     Double_t minThree[3] = {-0.5*TMath::Pi(), -0.5*TMath::Pi(),    minbinCentr};
 278     Double_t maxThree[3] = { 1.5*TMath::Pi(),  1.5*TMath::Pi(),    maxbinCentr};
 279
 280
 281
 282
 283
 284     //--------------------
 285     TString dataType[2] ={"ESD", "AOD"};
 286     TString labels[8]={Form("%sAllAllMcNmc",dataType[fMode].Data()),
 287         Form("%sTrigAllMcNmc",dataType[fMode].Data()),
 288         Form("%sTrigVtxMcNmc",dataType[fMode].Data()),
 289         Form("%sTrigVtxMcNrec",dataType[fMode].Data()),
 290         Form("%sTrigVtxRecMcPropNrec",dataType[fMode].Data()),
 291         Form("%sTrigVtxRecNrec",dataType[fMode].Data()),
 292         Form("%sTrigVtxRecMcPropNrecStrangeCorr",dataType[fMode].Data()),
 293         Form("%sTrigVtxRecNrecStrangeCorr",dataType[fMode].Data())};
 294
 295
 296     for(Int_t i=0;i<8;i++){
 297
 298         fMapSingleTrig[i] = new THnSparseD(Form("fMapSingleTrig%s", labels[i].Data()),"eta:pt:Nrec",
 299                                            3,binsEffHisto,minEffHisto,maxEffHisto);
 300         fMapSingleTrig[i]->SetBinEdges(1,binsPt);
 301         fMapSingleTrig[i]->GetAxis(0)->SetTitle("#eta");
 302         fMapSingleTrig[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
 303         fMapSingleTrig[i]->GetAxis(2)->SetTitle("N_{rec}");
 304         fMapSingleTrig[i]->Sumw2();
 305
 306         fMapPair[i] = new THnSparseD(Form("fMapPair%s", labels[i].Data()),"pt_trig:pt_assoc:DeltaEta:DeltaPhi:Nrec:likesign",
 307                                      6,binsEffHisto5,minEffHisto5,maxEffHisto5);
 308         fMapPair[i]->SetBinEdges(0,binsPt);
 309         fMapPair[i]->SetBinEdges(1,binsPt);
 310         fMapPair[i]->GetAxis(0)->SetTitle("p_{T, trig} (GeV/c)");
 311         fMapPair[i]->GetAxis(1)->SetTitle("p_{T, assoc} (GeV/c)");
 312         fMapPair[i]->GetAxis(2)->SetTitle("#Delta #eta");
 313         fMapPair[i]->GetAxis(3)->SetTitle("#Delta #phi");
 314         fMapPair[i]->GetAxis(4)->SetTitle("N_{rec}");
 315         fMapPair[i]->GetAxis(5)->SetTitle("Like-sign or Unlike-sign");
 316         fMapPair[i]->Sumw2();
 317
 318
 319         fMapEvent[i] = new THnSparseD(Form("fMapEvent%s", labels[i].Data()),"Nrec:vertexZ:meanPt:leadingPt",
 320                                       4,binsEvent,minEvent,maxEvent);
 321         fMapEvent[i]->GetAxis(0)->SetTitle("N_{rec}");
 322         fMapEvent[i]->GetAxis(1)->SetTitle("z_{vertex} (cm)");
 323         fMapEvent[i]->GetAxis(2)->SetTitle("meanPt");
 324         fMapEvent[i]->SetBinEdges(3,binsPt);
 325         fMapEvent[i]->GetAxis(3)->SetTitle("leadingPt");
 326         fMapEvent[i]->Sumw2();
 327
 328         fMapAll[i] = new THnSparseD(Form("fMapAll%s", labels[i].Data()),"eta:pt:Nrec",
 329                                     3,binsAll,minAll,maxAll);
 330         fMapAll[i]->SetBinEdges(1,binsPt);
 331         fMapAll[i]->GetAxis(0)->SetTitle("#eta");
 332         fMapAll[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
 333         fMapAll[i]->GetAxis(2)->SetTitle("N_{rec}");
 334         fMapAll[i]->Sumw2();
 335
 336
 337         fMapThree[i] = new THnSparseD(Form("fMapThree%s", labels[i].Data()),"dphi1:dphi2:Nrec",
 338                                       3,binsThree,minThree,maxThree);
 339         fMapThree[i]->GetAxis(0)->SetTitle("#Delta#varphi_{1}");
 340         fMapThree[i]->GetAxis(1)->SetTitle("#Delta#varphi_{2}");
 341         fMapThree[i]->GetAxis(2)->SetTitle("N_{rec}");
 342         fMapThree[i]->Sumw2();
 343
 344
 345         fVertexZ[i]                  = new TH1F(Form("fVertexZ%s",labels[i].Data()),
 346                                                 Form("fVertexZ%s",labels[i].Data()) ,
 347                                                 220, -11., 11.);
 348         fPt[i]                       = new TH1F(Form("fPt%s",labels[i].Data()),
 349                                                 Form("fPt%s",labels[i].Data()) ,
 350                                                 200, 0., 50);
 351         fNcharge[i]                  = new TH1F(Form("fNcharge%s",labels[i].Data()),
 352                                                 Form("fNcharge%s",labels[i].Data()) ,
 353                                                 250, -0.5, 249.5);
 354         fEta[i]                      = new TH1F (Form("fEta%s",labels[i].Data()),
 355                                                  Form("fEta%s",labels[i].Data()) ,
 356                                                  100, -2., 2);
 357         fPhi[i]                      = new TH1F(Form("fPhi%s",labels[i].Data()),
 358                                                 Form("fPhi%s",labels[i].Data()) ,
 359                                                 360, 0.,2*TMath::Pi());
 360         fDcaXY[i]                    = new TH1F(Form("fDcaXY%s",labels[i].Data()),
 361                                                 Form("fDcaXY%s",labels[i].Data()) ,
 362                                                 200, -0.3,0.3);
 363         fDcaZ[i]                     = new TH1F(Form("fDcaZ%s",labels[i].Data()),
 364                                                 Form("fDcaZ%s",labels[i].Data()) ,
 365                                                 200, -2.2,2.2);
 366         fPtSeed[i]                       = new TH1F(Form("fPSeedt%s",labels[i].Data()),
 367                                                     Form("fPtSeed%s",labels[i].Data()) ,
 368                                                     500, 0., 50);
 369         fEtaSeed[i]                      = new TH1F (Form("fEtaSeed%s",labels[i].Data()),
 370                                                      Form("fEtaSeed%s",labels[i].Data()) ,
 371                                                      100, -2., 2);
 372         fPhiSeed[i]                      = new TH1F(Form("fPhiSeed%s",labels[i].Data()),
 373                                                     Form("fPhiSeed%s",labels[i].Data()) ,
 374                                                     360, 0.,2*TMath::Pi());
 375
 376         fPtOthers[i]                       = new TH1F(Form("fPOtherst%s",labels[i].Data()),
 377                                                       Form("fPtOthers%s",labels[i].Data()) ,
 378                                                       500, 0., 50);
 379         fEtaOthers[i]                      = new TH1F (Form("fEtaOthers%s",labels[i].Data()),
 380                                                        Form("fEtaOthers%s",labels[i].Data()) ,
 381                                                        100, -2., 2);
 382         fPhiOthers[i]                      = new TH1F(Form("fPhiOthers%s",labels[i].Data()),
 383                                                       Form("fPhiOthers%s",labels[i].Data()) ,
 384                                                       360, 0.,2*TMath::Pi());
 385         fPtEtaOthers[i]                      = new TH2F(Form("fPtEtaOthers%s",labels[i].Data()),
 386                                                         Form("fPtEtaOthers%s",labels[i].Data()) ,
 387                                                         500, 0., 50, 100, -1., 1);
 388         fPhiEta[i]                   = new TH2F(Form("fPhiEta%s",labels[i].Data()),
 389                                                 Form("fPhiEta%s",labels[i].Data()) ,
 390                                                 180, 0., 2*TMath::Pi(), 100, -1.,1.);
 391         fDPhiDEtaEventAxis[i]        = new TH2F(Form("fDPhiDEtaEventAxis%s",labels[i].Data()),
 392                                                 Form("fDPhiDEtaEventAxis%s",labels[i].Data()) ,
 393                                                 180, -0.5* TMath::Pi(), 1.5*TMath::Pi(), 9, -1.8,1.8);
 394         fTriggerNch[i]               = new TH1F(Form("fTriggerNch%s",labels[i].Data()),
 395                                                 Form("fTriggerNch%s",labels[i].Data()) ,
 396                                                 250, -0.5, 249.5);
 397         fTriggerNchSeeds[i]          = new TH2F(Form("fTriggerNchSeeds%s",labels[i].Data()),
 398                                                 Form("fTriggerNchSeeds%s",labels[i].Data()) ,
 399                                                 250, -0.5, 249.5, 100, -0.5, 99.5);
 400         fTriggerTracklet[i]          = new TH1F(Form("fTriggerTracklet%s",labels[i].Data()),
 401                                                 Form("fTriggerTracklet%s",labels[i].Data()) ,
 402                                                 250, -0.5, 249.5);
 403         fNch07Nch[i]                 = new TH2F(Form("fNch07Nch%s",labels[i].Data()),
 404                                                 Form("fNch07Nch%s",labels[i].Data()) ,
 405                                                 250, -2.5, 247.5,250, -2.5, 247.5);
 406         fPNch07Nch[i]                 = new TProfile(Form("pNch07Nch%s",labels[i].Data()),
 407                                                      Form("pNch07Nch%s",labels[i].Data()) ,
 408                                                      250, -2.5, 247.5);
 409         fNch07Tracklet[i]            = new TH2F(Form("fNch07Tracklet%s",labels[i].Data()),
 410                                                 Form("fNch07Tracklet%s",labels[i].Data()) ,
 411                                                 250, -2.5, 247.5,250, -2.5, 247.5);
 412         fNchTracklet[i]              = new TH2F(Form("fNchTracklet%s",labels[i].Data()),
 413                                                 Form("fNchTracklet%s",labels[i].Data()) ,
 414                                                 250, -2.5, 247.5,250, -2.5, 247.5);
 415         fPNch07Tracklet[i]            = new TProfile(Form("pNch07Tracklet%s",labels[i].Data()),
 416                                                      Form("pNch07Tracklet%s",labels[i].Data()) ,
 417                                                      250, -2.5, 247.5);
 418         fDPhiEventAxis[i]          = new TH1F(Form("fDPhiEventAxis%s",
 419                                                    labels[i].Data()),
 420                                               Form("fDPhiEventAxis%s",
 421                                                    labels[i].Data()) ,
 422                                               180, -0.5*TMath::Pi(), 1.5*TMath::Pi());
 423         fDPhi1DPhi2[i] = new TH2F(Form("fDPhi1DPhi2%s",labels[i].Data()),
 424                                   Form("fDPhi1DPhi2%s",labels[i].Data()),
 425                                   180, -0.5* TMath::Pi(), 1.5*TMath::Pi(),
 426                                   180, -0.5* TMath::Pi(), 1.5*TMath::Pi());
 427     }
 428
 429     fHists = new TList();
 430     fHists->SetOwner();
 431
 432     fHists->Add(fStep);
 433     fHists->Add(fEventStat);
 434     fHists->Add(fHistPt);
 435     fHists->Add(fNContrNtracklets);
 436     fHists->Add(fNContrNtracks);
 437     fHists->Add(fCorruptedChunks);
 438     fHists->Add(fCorruptedChunksAfter);
 439
 440     if(fUseMC){
 441         fHists->Add(fHistPtMC);
 442
 443         fHists->Add(fNmcNch);
 444         fHists->Add(fPNmcNch);
 445         fHists->Add(fNmcNchVtx);
 446         fHists->Add(fNmcNchVtxStrangeCorr);
 447         fHists->Add(fPNmcNchVtx);
 448
 449         fHists->Add(fNmcNchTracklet);
 450         fHists->Add(fPNmcNchTracklet);
 451         fHists->Add(fNmcNchVtxTracklet);
 452         fHists->Add(fPNmcNchVtxTracklet);
 453     }
 454     fHists->Add(fChargedPi0);
 455     fHists->Add(fVertexCheck);
 456     fHists->Add(fPropagateDca);
 457
 458     for(Int_t i=0;i<8;i++){
 459         fHists->Add(fMapSingleTrig[i]);
 460         fHists->Add(fMapPair[i]);
 461         fHists->Add(fMapEvent[i]);
 462         fHists->Add(fMapAll[i]);
 463         fHists->Add(fMapThree[i]);
 464         fHists->Add(fVertexZ[i]);
 465         fHists->Add(fPt[i]);
 466         fHists->Add(fNcharge[i]);
 467         fHists->Add(fEta[i]);
 468         fHists->Add(fPhi[i]);
 469         fHists->Add(fDcaXY[i]);
 470         fHists->Add(fDcaZ[i]);
 471         fHists->Add(fPtSeed[i]);
 472         fHists->Add(fEtaSeed[i]);
 473         fHists->Add(fPhiSeed[i]);
 474         fHists->Add(fPtOthers[i]);
 475         fHists->Add(fEtaOthers[i]);
 476         fHists->Add(fPhiOthers[i]);
 477         fHists->Add(fPtEtaOthers[i]);
 478         fHists->Add(fPhiEta[i]);
 479         fHists->Add(fDPhiDEtaEventAxis[i]);
 480         fHists->Add(fTriggerNch[i]);
 481         fHists->Add(fTriggerNchSeeds[i]);
 482         fHists->Add(fTriggerTracklet[i]);
 483         fHists->Add(fNch07Nch[i]);
 484         fHists->Add(fPNch07Nch[i]);
 485         fHists->Add(fNch07Tracklet[i]);
 486         fHists->Add(fNchTracklet[i]);
 487         fHists->Add(fPNch07Tracklet[i]);
 488         fHists->Add(fDPhiEventAxis[i]);
 489         fHists->Add(fDPhi1DPhi2[i]);
 490     }
 491
 492     PostData(1, fHists);
 493
 494 }
 495
 496 //________________________________________________________________________
 497 void AliAnalysisTaskMinijet::UserExec(Option_t *)
 498 {
 499     // Main function, called for each event
 500     // Kinematics-only, ESD and AOD can be processed.
 501     // Data is read (ReadEventESD, ReadEventAOD...) and then analysed (Analyse).
 502     //  - in case of MC with full detector simulation, all correction steps(0-5) can be processed
 503     //  - for Data, only step 5 is performed
 504     //  - for kinematics-only, only step 0 is processed
 505
 506     //             - trigger -               - vertex -                       - tracks -                                       - multiplicity -
 507     // step 7 =  Triggered events, reconstructed accepted vertex,  reconstructed tracks + strangness corr,               reconstructed multiplicity+strangeCorr
 508     // step 6 =  Triggered events, reconstructed accepted vertex,  reconstructed tracks with MC prop + stangess corr,    reconstructed multiplicity+strangeCorr
 509     // step 5 =  Triggered events, reconstructed accepted vertex,  reconstructed tracks,                                 reconstructed multiplicity
 510     // step 4 =  Triggered events, reconstructed accepted vertex,  reconstructed tracks with MC properties,              reconstructed multiplicity
 511     // step 3 =  Triggered events, reconstructed accepted vertex,  mc primary particles,                                 reconstructed multiplicity
 512     // step 2 =  Triggered events, reconstructed accepted vertex,  mc primary particles,                                 true multiplicity
 513     // step 1 =  Triggered events, all                             mc primary particles,                                 true multiplicity
 514     // step 0 =  All events,       all                             mc primary particles,                                 true multiplicity
 515
 516
 517     if(fDebug) Printf("UserExec: Event starts");
 518
 519     AliVEvent *event = InputEvent();
 520     if (!event) {
 521         Error("UserExec", "Could not retrieve event");
 522         return;
 523     }
 524     fBSign= event->GetMagneticField();
 525
 526     //get events, either ESD or AOD
 527     fESDEvent = dynamic_cast<AliESDEvent*> (InputEvent());
 528     fAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
 529
 530     vector<Float_t> pt;
 531     vector<Float_t> eta;
 532     vector<Float_t> phi;
 533     vector<Short_t> charge;
 534     vector<Float_t> strangenessWeight;
 535     vector<Float_t> px;
 536     vector<Float_t> py;
 537     vector<Float_t> pz;
 538     vector<Float_t> theta;
 539
 540
 541     //number of accepted tracks and tracklets
 542     Double_t ntracks = 0;  //return value for reading functions for ESD and AOD
 543     //Int_t ntracksRemove = 0;  //return value for reading functions for ESD and AOD
 544     vector<Double_t> nTracksTracklets; // [0]=nAccepted,1=ntracklets,2=nall(also neutral in case of mc,
 545     //for real nall=ncharged)
 546
 547     if(!fAODEvent && !fESDEvent)return;
 548
 549     //Centrality
 550     Double_t centrality = 0;
 551     if (fCentralityMethod.Length() > 0)
 552     {
 553         AliCentrality *centralityObj = 0;
 554         if (fAODEvent)
 555             centralityObj = fAODEvent->GetHeader()->GetCentralityP();
 556         else if (fESDEvent)
 557             centralityObj = fESDEvent->GetCentrality();
 558         if (centralityObj)
 559         {
 560             centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
 561             AliInfo(Form("Centrality is %f", centrality));
 562         }
 563         else
 564         {
 565             Printf("WARNING: Centrality object is 0");
 566             centrality = -1;
 567         }
 568     }
 569
 570     // SDD check for LHC11a
 571     if (fCheckSDD) {
 572
 573         //ESD
 574         if(fESDEvent){
 575             if(fSelOption==0){
 576                 const AliMultiplicity *mul = fESDEvent->GetMultiplicity();
 577                 Int_t nClu3 = mul->GetNumberOfITSClusters(2);
 578                 Int_t nClu4 = mul->GetNumberOfITSClusters(3);
 579                 if(nClu3==0 &&  nClu4==0) return;
 580             }
 581             else if (fSelOption==1){
 582                 TString trcl = fESDEvent->GetFiredTriggerClasses().Data();
 583                 if (!(trcl.Contains("CINT1-B-NOPF-ALLNOTRD"))) return;
 584             }
 585         }
 586
 587         //AOD
 588         if(fAODEvent){
 589             if(fSelOption==0){
 590                 Bool_t useEvent = false;
 591                 Int_t nTracks = fAODEvent->GetNTracks();
 592                 for(Int_t itrack=0; itrack<nTracks; itrack++) {
 593                     AliAODTrack * track = fAODEvent->GetTrack(itrack);
 594                     if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
 595                         useEvent=true;
 596                         break;
 597                     }
 598                 }
 599                 if (!useEvent) return;
 600             }
 601             else if(fSelOption==1){
 602                 TString trcl = fAODEvent->GetFiredTriggerClasses().Data();
 603                 if (!(trcl.Contains("CINT1-B-NOPF-ALLNOTRD"))) return;
 604             }
 605         }
 606     }
 607
 608     //reset values
 609     fNMcPrimAccept=0;// number of accepted primaries
 610     fNRecAccept=0;   // number of accepted tracks
 611     fNRecAcceptStrangeCorr=0;   // number of accepted tracks + strangeness correction
 612     fNMcPrimAcceptTracklet=0;// number of accepted primaries (no pt cut)
 613     fNRecAcceptTracklet=0;   // number of accepted tracklets
 614
 615     // instead of task->SelectCollisionCandidate(mask) in AddTask macro
 616     Bool_t isSelected = (((AliInputEventHandler*)
 617                           (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 618                          ->IsEventSelected() &  fTriggerType);
 619
 620
 621     if(fDebug){
 622         Printf("IsSelected = %d", isSelected);
 623         Printf("CheckEvent(true)= %d", CheckEvent(true));
 624         Printf("CheckEvent(false)= %d", CheckEvent(false));
 625     }
 626
 627     fEventStat->Fill(0);//all events
 628
 629     //check trigger
 630     if(isSelected){ // has offline trigger
 631
 632         fEventStat->Fill(1);//triggered event
 633
 634         if(CheckEvent(true)){//step 5 = TrigVtxRecNrec, step 4 = TrigVtxRecMcPropNrec ,step 3 = TrigVtxMcNrec
 635
 636             fEventStat->Fill(2);//triggered event with vertex
 637
 638             if(!fMcOnly){
 639                 //step 5 = TrigVtxRecNrec
 640
 641                 // read tracks
 642                 if(fESDEvent)     ntracks = ReadEventESD(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 5);
 643                 else if(fAODEvent)ntracks = ReadEventAOD(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 5);
 644                 else AliInfo("Fatal Error");
 645
 646                 if (fCentralityMethod.Length() > 0)
 647                     ntracks = centrality;
 648
 649                 // analyse
 650                 if(pt.size()){ //(internally ntracks=fNRecAccept)
 651                     fEventStat->Fill(3);//triggered event with vertex and one reconstructed track in acceptance
 652                         Analyse(pt, eta, phi, charge, strangenessWeight, ntracks, nTracksTracklets[1], nTracksTracklets[2], 5);//step 5 = TrigVtxRecNrec
 653                     }
 654
 655                 if(fCorrStrangeness){
 656                     //step 7 = TrigVtxRecNrecStrangeCorr
 657
 658                     // read tracks
 659                     if(fESDEvent)     ntracks = ReadEventESD(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 7);//stagness version not yet implemented
 660                     else if(fAODEvent)ntracks = ReadEventAOD(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 7);
 661                     else AliInfo("Fatal Error");
 662
 663                     if (fCentralityMethod.Length() > 0)
 664                         ntracks = centrality;
 665
 666                     // analyse
 667                     if(pt.size()){ //(internally ntracks=fNRecAccept)
 668                         Analyse(pt, eta, phi, charge, strangenessWeight, fNRecAcceptStrangeCorr, nTracksTracklets[1], nTracksTracklets[2], 7);//step 7 = TrigVtxRecNrecStrangeCorr
 669                     }
 670                 }
 671
 672                 if(fUseMC){
 673                     // step 4 = TrigVtxRecMcPropNrec
 674
 675                     // read tracks
 676                     if(fESDEvent)       ntracks = ReadEventESDRecMcProp(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 4);
 677                     else if(fAODEvent)  ntracks = ReadEventAODRecMcProp(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 4);
 678                     else AliInfo("Fatal Error");
 679
 680                     if (fCentralityMethod.Length() > 0)
 681                         ntracks = centrality;
 682
 683                     //analyse
 684                     if(pt.size()){//(internally ntracks=fNRecAccept)
 685                         Analyse(pt, eta, phi, charge,strangenessWeight, ntracks, nTracksTracklets[1], nTracksTracklets[2], 4); //step 4 = TrigVtxRecMcPropNrec
 686                     }
 687
 688
 689                     if(fCorrStrangeness){
 690                         // step 6 = TrigVtxRecMcPropNrecStrangeCorr
 691
 692                         // read tracks
 693                         if(fESDEvent)       ntracks = ReadEventESDRecMcProp(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 6);//stagness version not yet implemented
 694                         else if(fAODEvent)  ntracks = ReadEventAODRecMcProp(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 6);
 695                         else AliInfo("Fatal Error");
 696
 697                         if (fCentralityMethod.Length() > 0)
 698                             ntracks = centrality;
 699
 700                         //analyse
 701                         if(pt.size()){//(internally ntracks=fNRecAccept)
 702                             Analyse(pt, eta, phi, charge, strangenessWeight, fNRecAcceptStrangeCorr, nTracksTracklets[1], nTracksTracklets[2], 6); //step 6 = TrigVtxRecMcPropNrecStrangeCorr
 703                         }
 704                     }
 705                     // step 3 = TrigVtxMcNrec
 706
 707                     // read tracks
 708                     if(fESDEvent)       ntracks = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 3);
 709                     else if(fAODEvent)  ntracks = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 3);
 710                     else AliInfo("Fatal Error");
 711
 712                     if (fCentralityMethod.Length() > 0){
 713                         fNRecAccept = centrality;
 714                         fNMcPrimAccept = centrality;
 715                     }
 716
 717                     // analyse
 718                     if(pt.size()){
 719                         Analyse(pt, eta, phi, charge, strangenessWeight, fNRecAccept,    nTracksTracklets[1],nTracksTracklets[2], 3); //step 3 = TrigVtxMcNrec
 720                         Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 2); //step 2 = TrigVtxMcNmc
 721                     }
 722
 723                 }
 724             }
 725
 726         }//check event (true)
 727
 728
 729         if(fUseMC && !fMcOnly){
 730             //reset values
 731             fNMcPrimAccept=0;// number of accepted primaries
 732             fNRecAccept=0;   // number of accepted tracks
 733             fNMcPrimAcceptTracklet=0;// number of accepted primaries (no pt cut)
 734             fNRecAcceptTracklet=0;   // number of accepted tracklets
 735
 736             if(CheckEvent(false)){// all events, with and without reconstucted vertex
 737                 if(fESDEvent) ntracks       = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 1);//read tracks
 738                 else if(fAODEvent) ntracks  = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 1);//read tracks
 739                 else AliInfo("Fatal Error");
 740
 741                 if (fCentralityMethod.Length() > 0)
 742                     fNMcPrimAccept = centrality;
 743
 744
 745                 // analyse
 746                 if(pt.size()){
 747                     Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 1);  // step 1 = TrigAllMcNmc
 748
 749                     Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);  //first part of step 0 // step 0 = AllAllMcNmc
 750                 }
 751
 752
 753             }
 754         }
 755     }// triggered event
 756
 757     else { // not selected by physics selection task = not triggered
 758         if(fUseMC && !fMcOnly){
 759             if(CheckEvent(false)){
 760
 761                 //read tracks
 762                 if(fESDEvent)      ntracks  = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
 763                 else if(fAODEvent) ntracks  = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
 764                 else AliInfo("Fatal Error");
 765
 766                 if (fCentralityMethod.Length() > 0)
 767                     fNMcPrimAccept = centrality;
 768
 769                 //analyse
 770                 if(pt.size()){
 771                     Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);  //second part of step 0 // step 0 = AllAllMcNmc
 772                 }
 773             }
 774         }
 775     }
 776
 777     if(fMcOnly){
 778         // read event
 779         if(fMode==0)       ntracks  = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
 780         else if (fMode==1) ntracks  = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
 781
 782         if (fCentralityMethod.Length() > 0)
 783             fNMcPrimAccept = centrality;
 784
 785         // analyse
 786         if(pt.size()){
 787             Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);
 788         }
 789     }
 790 }
 791
 792
 793 //________________________________________________________________________
 794 Double_t AliAnalysisTaskMinijet::ReadEventESD( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 795                                            vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
 796                                            vector<Float_t> &strangeArray,
 797                                            vector<Double_t>   &nTracksTracklets, const Int_t step)
 798 {
 799     // gives back the number of esd tracks and pointer to arrays with track
 800     // properties (pt, eta, phi)
 801     // Uses TPC tracks with SPD vertex from now on
 802
 803     ptArray.clear();
 804     etaArray.clear();
 805     phiArray.clear();
 806     chargeArray.clear();
 807     strangeArray.clear();
 808     nTracksTracklets.clear();
 809
 810     const AliESDVertex* vtxSPD   = fESDEvent->GetPrimaryVertexSPD(); // uses track or SPD vertexer
 811     fVertexZ[step]->Fill(vtxSPD->GetZ());
 812
 813     // Retreive the number of all tracks for this event.
 814     Int_t ntracks = fESDEvent->GetNumberOfTracks();
 815     if(fDebug>1)  Printf("all ESD tracks: %d", ntracks);
 816
 817     //first loop to check how many tracks are accepted
 818     //------------------
 819     Double_t nAcceptedTracks=0;
 820     //Float_t nAcceptedTracksStrange=0;
 821     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 822
 823         AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
 824         if (!esdTrack) {
 825             Error("ReadEventESD", "Could not receive track %d", iTracks);
 826             continue;
 827         }
 828
 829         // use TPC only tracks with non default SPD vertex
 830         AliESDtrack *track = AliESDtrackCuts::
 831         GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
 832         if(!track) continue;
 833         if(!fCuts->AcceptTrack(track)) {
 834             delete track;
 835             continue;// apply TPC track cuts before constrain to SPD vertex
 836         }
 837         if(track->Pt()>0.){
 838             // only constrain tracks above threshold
 839             AliExternalTrackParam exParam;
 840             // take the B-field from the ESD, no 3D fieldMap available at this point
 841             Bool_t relate = false;
 842             relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
 843                                               kVeryBig,&exParam);
 844             if(!relate){
 845                 delete track;
 846                 continue;
 847             }
 848             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
 849                        exParam.GetCovariance());
 850         }
 851         else{
 852             delete track;
 853             continue;// only if tracks have pt<=0
 854         }
 855
 856         if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax) {
 857             ptArray.push_back(track->Pt());
 858             etaArray.push_back(track->Eta());
 859             phiArray.push_back(track->Phi());
 860             chargeArray.push_back(track->Charge());
 861             strangeArray.push_back(1);
 862             ++nAcceptedTracks;
 863             fHistPt->Fill(track->Pt());
 864         }
 865
 866         // TPC only track memory needs to be cleaned up
 867         if(track)
 868             delete track;
 869
 870     }
 871
 872     //need to be checked
 873     if(nAcceptedTracks==0) return 0;
 874
 875     //tracklet loop
 876     Int_t ntrackletsAccept=0;
 877     AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
 878     Int_t ntracklets = mult->GetNumberOfTracklets();
 879     for(Int_t i=0;i< ntracklets;i++){
 880         if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
 881             ntrackletsAccept++;
 882         }
 883     }
 884     nTracksTracklets.push_back(nAcceptedTracks);
 885     nTracksTracklets.push_back(ntrackletsAccept);
 886     nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
 887     //where here also neutral particles are counted.
 888
 889
 890     fVzEvent=vtxSPD->GetZ(); // needed for correction map
 891     if(step==5 || step ==2) {
 892         fNRecAccept=nAcceptedTracks;
 893         fNRecAcceptTracklet=ntrackletsAccept;
 894     }
 895     return fNRecAccept;
 896
 897
 898 }
 899
 900 //________________________________________________________________________
 901 Double_t AliAnalysisTaskMinijet::ReadEventESDRecMcProp( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 902                                                     vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
 903                                                     vector<Float_t> &strangeArray,
 904                                                     vector<Double_t> &nTracksTracklets, const Int_t step)
 905 {
 906     // gives back the number of esd tracks and pointer to arrays with track
 907     // properties (pt, eta, phi) of mc particles if available
 908     // Uses TPC tracks with SPD vertex from now on
 909
 910     ptArray.clear();
 911     etaArray.clear();
 912     phiArray.clear();
 913     chargeArray.clear();
 914     strangeArray.clear();
 915     nTracksTracklets.clear();
 916
 917
 918     AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
 919     if (!mcEvent) {
 920         Error("ReadEventESDRecMcProp", "Could not retrieve MC event");
 921         return 0;
 922     }
 923     AliStack* stack = MCEvent()->Stack();
 924     if(!stack) return 0;
 925
 926
 927     // Retreive the number of all tracks for this event.
 928     Int_t ntracks = fESDEvent->GetNumberOfTracks();
 929     if(fDebug>1)Printf("all ESD tracks: %d", ntracks);
 930
 931     const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
 932     fVertexZ[step]->Fill(vtxSPD->GetZ());
 933
 934     //track loop
 935     Double_t nAcceptedTracks=0;
 936     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 937
 938         AliVParticle *vtrack = fESDEvent->GetTrack(iTracks);
 939         AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
 940         if (!esdTrack) {
 941             Error("ReadEventESDRecMcProp", "Could not receive track %d", iTracks);
 942             continue;
 943         }
 944
 945         // use TPC only tracks with non default SPD vertex
 946         AliESDtrack *track = AliESDtrackCuts::
 947         GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
 948         if(!track) continue;
 949         if(!fCuts->AcceptTrack(track)) {
 950             delete track;
 951             continue;// apply TPC track cuts before constrain to SPD vertex
 952         }
 953         if(track->Pt()>0.){
 954             // only constrain tracks above threshold
 955             AliExternalTrackParam exParam;
 956             // take the B-field from the ESD, no 3D fieldMap available at this point
 957             Bool_t relate = false;
 958             relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
 959                                               kVeryBig,&exParam);
 960             if(!relate){
 961                 delete track;
 962                 continue;
 963             }
 964             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
 965                        exParam.GetCovariance());
 966         }
 967         else{
 968             delete track;
 969             continue;// only if tracks have pt<=0
 970         }
 971
 972         //count tracks, if available, use mc particle properties
 973         if(vtrack->GetLabel()<=0){
 974             if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax){
 975                 ptArray.push_back(track->Pt());
 976                 etaArray.push_back(track->Eta());
 977                 phiArray.push_back(track->Phi());
 978                 chargeArray.push_back(track->Charge());
 979                 strangeArray.push_back(1);
 980                 ++nAcceptedTracks;
 981             }
 982         }
 983         else{
 984             TParticle *partOfTrack = stack->Particle(vtrack->GetLabel());
 985             if (TMath::Abs(partOfTrack->Eta())<fEtaCut && partOfTrack->Pt()>fPtMin && partOfTrack->Pt()<fPtMax) {
 986                 ptArray.push_back(partOfTrack->Pt());
 987                 etaArray.push_back(partOfTrack->Eta());
 988                 phiArray.push_back(partOfTrack->Phi());
 989                 chargeArray.push_back(vtrack->Charge());
 990                 strangeArray.push_back(1);
 991                 ++nAcceptedTracks;
 992             }
 993         }
 994
 995         // TPC only track memory needs to be cleaned up
 996         if(track)
 997             delete track;
 998
 999     }
1000
1001     if(nAcceptedTracks==0) return 0;
1002
1003     //tracklet loop
1004     Int_t ntrackletsAccept=0;
1005     AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
1006     Int_t ntracklets = mult->GetNumberOfTracklets();
1007     for(Int_t i=0;i< ntracklets;i++){
1008         if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
1009             ntrackletsAccept++;
1010         }
1011     }
1012
1013     nTracksTracklets.push_back(nAcceptedTracks);
1014     nTracksTracklets.push_back(ntrackletsAccept);
1015     nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1016     //where here also neutral particles are counted.
1017
1018
1019     //get mc vertex for correction maps
1020     AliGenEventHeader*  header = MCEvent()->GenEventHeader();
1021     TArrayF mcV;
1022     header->PrimaryVertex(mcV);
1023     fVzEvent= mcV[2];
1024
1025     return fNRecAccept; // give back reconstructed value
1026
1027
1028 }
1029
1030
1031
1032
1033 //________________________________________________________________________
1034 Double_t AliAnalysisTaskMinijet::ReadEventESDMC(vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1035                                              vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1036                                              vector<Float_t> &strangeArray,
1037                                              vector<Double_t> &nTracksTracklets, const Int_t step)
1038 {
1039     // gives back the number of charged prim MC particle and pointer to arrays
1040     // with particle properties (pt, eta, phi)
1041
1042     ptArray.clear();
1043     etaArray.clear();
1044     phiArray.clear();
1045     chargeArray.clear();
1046     strangeArray.clear();
1047     nTracksTracklets.clear();
1048
1049     fNMcPrimAccept=0;
1050
1051     AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
1052     if (!mcEvent) {
1053         Error("ReadEventESDMC", "Could not retrieve MC event");
1054         return 0;
1055     }
1056
1057     AliStack* stack = MCEvent()->Stack();
1058     if(!stack) return 0;
1059
1060     Int_t ntracks = mcEvent->GetNumberOfTracks();
1061     if(fDebug>1)Printf("MC particles: %d", ntracks);
1062
1063     //vertex
1064     AliGenEventHeader*  header = MCEvent()->GenEventHeader();
1065     TArrayF mcV;
1066     Float_t vzMC=0.;
1067     if(header){
1068         header->PrimaryVertex(mcV);
1069         vzMC = mcV[2];
1070         if(step==1){
1071             fVertexZ[0]->Fill(vzMC);
1072         }
1073         fVertexZ[step]->Fill(vzMC);
1074     }
1075
1076     //----------------------------------
1077     //first track loop to check how many chared primary tracks are available
1078     Double_t nChargedPrimaries=0;
1079     Double_t nAllPrimaries=0;
1080
1081     Double_t nPseudoTracklets=0;
1082     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1083         AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
1084         if (!track) {
1085             Error("ReadEventESDMC", "Could not receive track %d", iTracks);
1086             continue;
1087         }
1088
1089
1090         if(//count also charged particles in case of fSelectParticles==2 (only neutral)
1091            !SelectParticlePlusCharged(
1092                                       track->Charge(),
1093                                       track->PdgCode(),
1094                                       stack->IsPhysicalPrimary(track->Label())
1095                                       )
1096            )
1097             continue;
1098
1099         //count number of pseudo tracklets
1100         if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.0) ++nPseudoTracklets; //0.035
1101         //same cuts as on ESDtracks
1102         if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1103
1104         //count all primaries
1105         ++nAllPrimaries;
1106         //count charged primaries
1107         if (track->Charge() != 0) ++nChargedPrimaries;
1108
1109         if(fDebug>2) Printf("PDG=%d, IsPrim=%d",  track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
1110
1111
1112     }
1113     //----------------------------------
1114
1115     if(fDebug>2){
1116         Printf("All in acceptance=%f",nAllPrimaries);
1117         Printf("Charged in acceptance =%f",nChargedPrimaries);
1118     }
1119
1120     fChargedPi0->Fill(nAllPrimaries,nChargedPrimaries);
1121
1122     if(nAllPrimaries==0) return 0;
1123     if(nChargedPrimaries==0) return 0;
1124
1125     //track loop
1126     //Int_t iChargedPiK=0;
1127     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1128         AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
1129         if (!track) {
1130             Error("ReadEventESDMC", "Could not receive track %d", iTracks);
1131             continue;
1132         }
1133
1134         if(!SelectParticle(
1135                            track->Charge(),
1136                            track->PdgCode(),
1137                            stack->IsPhysicalPrimary(track->Label())
1138                            )
1139            )
1140             continue;
1141
1142
1143         //same cuts as on ESDtracks
1144         if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin  || track->Pt()>fPtMax) continue;
1145
1146         if(fDebug>2) Printf("After: PDG=%d, IsPrim=%d",  track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
1147
1148
1149         fHistPtMC->Fill(track->Pt());
1150         //fills arrays with track properties
1151         ptArray.push_back(track->Pt());
1152         etaArray.push_back(track->Eta());
1153         phiArray.push_back(track->Phi());
1154         chargeArray.push_back(track->Charge());
1155         strangeArray.push_back(1);
1156
1157     } //track loop
1158
1159     nTracksTracklets.push_back(nChargedPrimaries);
1160     nTracksTracklets.push_back(nPseudoTracklets);
1161     if(fSelectParticles!=2){
1162         nTracksTracklets.push_back(nAllPrimaries);
1163     }
1164     else{
1165         nTracksTracklets.push_back(nAllPrimaries-nChargedPrimaries); // only neutral
1166     }
1167
1168     fNMcPrimAccept = nChargedPrimaries;
1169     fNMcPrimAcceptTracklet = nPseudoTracklets;
1170
1171     if(step==1){
1172         fNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1173         fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1174         fNmcNchTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1175         fPNmcNchTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1176     }
1177     if(step==3){
1178         fNmcNchVtx->Fill(fNMcPrimAccept,fNRecAccept);
1179         fPNmcNchVtx->Fill(fNMcPrimAccept,fNRecAccept);
1180         fNmcNchVtxTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1181         fPNmcNchVtxTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1182     }
1183
1184     fVzEvent= vzMC;
1185     return fNRecAccept;
1186
1187 }
1188
1189 //________________________________________________________________________
1190 Double_t AliAnalysisTaskMinijet::ReadEventAOD( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1191                                            vector<Float_t> &phiArray,  vector<Short_t> &chargeArray,
1192                                            vector<Float_t> &strangeArray,
1193                                            vector<Double_t> &nTracksTracklets, const Int_t step)
1194 {
1195     // gives back the number of AOD tracks and pointer to arrays with track
1196     // properties (pt, eta, phi)
1197
1198     ptArray.clear();
1199     etaArray.clear();
1200     phiArray.clear();
1201     chargeArray.clear();
1202     strangeArray.clear();
1203     nTracksTracklets.clear();
1204
1205     TClonesArray *mcArray=0x0;
1206     if(fAnalysePrimOnly || (fCorrStrangeness && fUseMC)){
1207         mcArray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
1208     }
1209
1210
1211     AliAODVertex*       vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1212     Double_t vzAOD=vertex->GetZ();
1213     fVertexZ[step]->Fill(vzAOD);
1214
1215     // Retreive the number of tracks for this event.
1216     Int_t ntracks = fAODEvent->GetNumberOfTracks();
1217     if(fDebug>1) Printf("AOD tracks: %d", ntracks);
1218
1219
1220     Double_t nAcceptedTracks=0;
1221     Float_t nAcceptedTracksStrange=0;
1222     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1223         AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
1224         if (!track) {
1225             Error("ReadEventAOD", "Could not receive track %d", iTracks);
1226             continue;
1227         }
1228
1229         AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
1230
1231         //use only tracks from primaries
1232         if(fAnalysePrimOnly){
1233             if(vtrack->GetLabel()<=0)continue;
1234             if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
1235         }
1236
1237         Double_t save= track->Pt();
1238         Double_t d0rphiz[2],covd0[3];
1239         Bool_t isDca= track->PropagateToDCA(fAODEvent->GetPrimaryVertex(),fAODEvent->GetMagneticField(),9999.,d0rphiz,covd0);
1240         fPropagateDca->Fill(Int_t(isDca));
1241         if(TMath::Abs(save - track->Pt())>1e-6) Printf("Before pt=%f, After pt=%f",save, track->Pt());
1242
1243         if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut
1244            && track->Pt()>fPtMin && track->Pt()<fPtMax){
1245
1246             nAcceptedTracks++;
1247
1248             // Printf("dca= %f", track->DCA());
1249             //save track properties in vector
1250             ptArray.push_back(track->Pt());
1251             etaArray.push_back(track->Eta());
1252             phiArray.push_back(track->Phi());
1253             chargeArray.push_back(track->Charge());
1254             fHistPt->Fill(track->Pt());
1255
1256
1257             //correction for underestimation of strangeness in Monte Carlos -> underestimation of contamination
1258             Float_t factor=1.;
1259             if(fUseMC && fCorrStrangeness && step==7){
1260                 if(vtrack->GetLabel()>0){
1261                     AliAODMCParticle* mcprong =(AliAODMCParticle*)mcArray->At(vtrack->GetLabel());
1262                     if(mcprong){
1263                         Int_t labmom = mcprong->GetMother();
1264                         if(labmom>=0){
1265                             AliAODMCParticle* mcmother=(AliAODMCParticle*)mcArray->At(labmom);
1266                             Int_t pdgMother=0;
1267                             if(mcmother) {
1268                                 pdgMother = mcmother->GetPdgCode();
1269                                 if(TMath::Abs(pdgMother)==310 || TMath::Abs(pdgMother)==130 || TMath::Abs(pdgMother)==321){ //K^0_S, K^0_L, K^+-
1270                                     if(track->Pt()<=1)factor=1./0.7; // values from strangeness publication
1271                                     else factor=1./0.6;// values from strangeness publication
1272                                 }
1273                                 if(TMath::Abs(pdgMother)==3122) { //Lambda
1274                                     factor=1./0.25; // values from strangeness publication
1275                                 }
1276                             }
1277                         }
1278                     }
1279                 }
1280             }
1281             nAcceptedTracksStrange+=factor;
1282             strangeArray.push_back(factor);
1283             fDcaXY[step]->Fill(d0rphiz[0], factor);
1284             fDcaZ[step]->Fill(d0rphiz[0], factor);
1285
1286         }
1287     }
1288     //need to check this option for MC
1289     if(nAcceptedTracks==0) return 0;
1290
1291
1292     //tracklet loop
1293     Int_t ntrackletsAccept=0;
1294     AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1295     for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1296         if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut){
1297             ++ntrackletsAccept;
1298         }
1299     }
1300
1301
1302     nTracksTracklets.push_back(nAcceptedTracks);
1303     nTracksTracklets.push_back(ntrackletsAccept);
1304     nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1305     //where here also neutral particles are counted.
1306
1307
1308     fVzEvent= vzAOD;
1309     if(step==5 || step==2){
1310         fNRecAccept = nAcceptedTracks; // needed for MC case //step5 = TrigVtxRecNrec
1311         fNRecAcceptTracklet = ntrackletsAccept; // needed for MC case //step5 = TrigVtxRecNrec
1312     }
1313     if(step==7)fNRecAcceptStrangeCorr = nAcceptedTracksStrange;
1314
1315     return fNRecAccept; // at the moment, always return reconstructed multiplicity
1316
1317 }
1318
1319 //________________________________________________________________________
1320 Double_t AliAnalysisTaskMinijet::ReadEventAODRecMcProp( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1321                                                     vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1322                                                     vector<Float_t> &strangeArray,
1323                                                     vector<Double_t> &nTracksTracklets, const Int_t step)
1324 {
1325     // gives back the number of AOD tracks and pointer to arrays with track
1326     // properties (pt, eta, phi)
1327
1328
1329     ptArray.clear();
1330     etaArray.clear();
1331     phiArray.clear();
1332     chargeArray.clear();
1333     strangeArray.clear();
1334     nTracksTracklets.clear();
1335
1336
1337     // Retreive the number of tracks for this event.
1338     Int_t ntracks = fAODEvent->GetNumberOfTracks();
1339     if(fDebug>1) Printf("AOD tracks: %d", ntracks);
1340
1341
1342     //get array of mc particles
1343     TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1344     FindListObject(AliAODMCParticle::StdBranchName());
1345     if(!mcArray){
1346         AliInfo("No MC particle branch found");
1347         return kFALSE;
1348     }
1349
1350     AliAODVertex*       vtx= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1351     Double_t vzAOD=vtx->GetZ();
1352     fVertexZ[step]->Fill(vzAOD);
1353
1354     Double_t nAcceptedTracks=0;
1355     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1356         AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
1357
1358         AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
1359
1360         if (!track) {
1361             Error("ReadEventAODRecMcProp", "Could not receive track %d", iTracks);
1362             continue;
1363         }
1364
1365         //use only tracks from primaries
1366         if(fAnalysePrimOnly){
1367             if(vtrack->GetLabel()<=0)continue;
1368             if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
1369         }
1370
1371         if(track->TestFilterBit(fFilterBit) &&  TMath::Abs(track->Eta())<fEtaCut &&
1372            track->Pt()>fPtMin && track->Pt()<fPtMax){
1373
1374             nAcceptedTracks++;
1375             Float_t factor=1.;
1376
1377             //save track properties in vector
1378             if(vtrack->GetLabel()<=0){ //fake tracks before "label<0", but crash in AOD079 // what is the meaning of label 0
1379                 //      Printf("Fake track");
1380                 //      continue;
1381                 ptArray.push_back(track->Pt());
1382                 etaArray.push_back(track->Eta());
1383                 phiArray.push_back(track->Phi());
1384                 chargeArray.push_back(track->Charge());
1385
1386             }
1387             else{//mc properties
1388                 AliAODMCParticle *partOfTrack = (AliAODMCParticle*)mcArray->At(vtrack->GetLabel());
1389                 if(!partOfTrack) Printf("label=%d", vtrack->GetLabel());
1390                 if(partOfTrack){
1391                     ptArray.push_back(partOfTrack->Pt());
1392                     etaArray.push_back(partOfTrack->Eta());
1393                     phiArray.push_back(partOfTrack->Phi());
1394                     chargeArray.push_back(vtrack->Charge());//partOfTrack?
1395
1396                     //correction for underestimation of strangeness in Monte Carlos -> underestimation of contamination
1397                     if(fUseMC && fCorrStrangeness && step==6){
1398                         Int_t labmom = partOfTrack->GetMother();
1399                         if(labmom>=0){
1400                             AliAODMCParticle* mcmother=(AliAODMCParticle*)mcArray->At(labmom);
1401                             Int_t pdgMother=0;
1402                             if(mcmother) {
1403                                 pdgMother = mcmother->GetPdgCode();
1404                                 if(TMath::Abs(pdgMother)==310 || TMath::Abs(pdgMother)==130 || TMath::Abs(pdgMother)==321){ //K^0_S, K^0_L, K^+-
1405                                     if(track->Pt()<=1)factor=1./0.7; // values from strangeness publication
1406                                     else factor=1./0.6;// values from strangeness publication
1407                                 }
1408                                 if(TMath::Abs(pdgMother)==3122) { //Lambda
1409                                     factor=1./0.25; // values from strangeness publication
1410                                 }
1411                             }
1412                         }
1413                     }
1414                 }
1415             }
1416             strangeArray.push_back(factor);
1417
1418         }
1419     }
1420     //need to check this option for MC
1421     if(nAcceptedTracks==0) return 0;
1422
1423     //tracklet loop
1424     Int_t ntrackletsAccept=0;
1425     AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1426     for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1427         if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut ){
1428             ++ntrackletsAccept;
1429         }
1430     }
1431
1432
1433     nTracksTracklets.push_back(nAcceptedTracks);
1434     nTracksTracklets.push_back(ntrackletsAccept);
1435     nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1436     //where here also neutral particles are counted.
1437
1438
1439     //check vertex (mc)
1440     AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1441     FindListObject(AliAODMCHeader::StdBranchName());
1442     Float_t vzMC = aodMCheader->GetVtxZ();
1443
1444     fVzEvent= vzMC;
1445     return fNRecAccept;//this is the rec value from step 5
1446
1447 }
1448
1449
1450
1451 //________________________________________________________________________
1452 Double_t AliAnalysisTaskMinijet::ReadEventAODMC( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1453                                              vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1454                                              vector<Float_t> &strangeArray,
1455                                              vector<Double_t> &nTracksTracklets, const Int_t step)
1456 {
1457     // gives back the number of AOD MC particles and pointer to arrays with particle
1458     // properties (pt, eta, phi)
1459
1460     ptArray.clear();
1461     etaArray.clear();
1462     phiArray.clear();
1463     chargeArray.clear();
1464     strangeArray.clear();
1465     nTracksTracklets.clear();
1466
1467     //check vertex
1468     AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1469     FindListObject(AliAODMCHeader::StdBranchName());
1470     Float_t vzMC = aodMCheader->GetVtxZ();
1471     if(step==1){
1472         fVertexZ[0]->Fill(vzMC);
1473     }
1474     fVertexZ[step]->Fill(vzMC);
1475
1476     //retreive MC particles from event
1477     TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1478     FindListObject(AliAODMCParticle::StdBranchName());
1479     if(!mcArray){
1480         AliInfo("No MC particle branch found");
1481         return kFALSE;
1482     }
1483
1484     Int_t ntracks = mcArray->GetEntriesFast();
1485     if(fDebug>1)Printf("MC particles: %d", ntracks);
1486
1487
1488     // Track loop: chek how many particles will be accepted
1489     //Float_t vzMC=0.;
1490     Double_t nChargedPrim=0;
1491     Double_t nAllPrim=0;
1492     Double_t nPseudoTracklets=0;
1493     for (Int_t it = 0; it < ntracks; it++) {
1494         AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1495         if (!track) {
1496             Error("ReadEventAODMC", "Could not receive particle %d", it);
1497             continue;
1498         }
1499
1500         if(!SelectParticlePlusCharged(
1501                                       track->Charge(),
1502                                       track->PdgCode(),
1503                                       track->IsPhysicalPrimary()
1504                                       )
1505            )
1506             continue;
1507
1508         if(TMath::Abs(track->Eta())<fEtaCut && track->Pt()>0.0)++nPseudoTracklets; //0.035
1509         if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1510
1511         nAllPrim++;
1512         if(track->Charge()!=0) nChargedPrim++;
1513
1514     }
1515
1516
1517     if(nAllPrim==0) return 0;
1518     if(nChargedPrim==0) return 0;
1519
1520     //generate array with size of number of accepted tracks
1521     fChargedPi0->Fill(nAllPrim,nChargedPrim);
1522
1523
1524     // Track loop: fill arrays for accepted tracks
1525     // Int_t iChargedPrim=0;
1526     for (Int_t it = 0; it < ntracks; it++) {
1527         AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1528         if (!track) {
1529             Error("ReadEventAODMC", "Could not receive particle %d", it);
1530             continue;
1531         }
1532
1533         if(!SelectParticle(
1534                            track->Charge(),
1535                            track->PdgCode(),
1536                            track->IsPhysicalPrimary()
1537                            )
1538            )
1539             continue;
1540
1541         if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1542
1543         fHistPtMC->Fill(track->Pt());
1544         ptArray.push_back(track->Pt());
1545         etaArray.push_back(track->Eta());
1546         phiArray.push_back(track->Phi());
1547         chargeArray.push_back(track->Charge());
1548         strangeArray.push_back(1);
1549     }
1550
1551     nTracksTracklets.push_back(nChargedPrim);
1552     nTracksTracklets.push_back(nPseudoTracklets);
1553     if(fSelectParticles!=2){
1554         nTracksTracklets.push_back(nAllPrim);
1555     }
1556     else{
1557         nTracksTracklets.push_back(nAllPrim-nChargedPrim); // only neutral
1558     }
1559
1560
1561
1562     fVzEvent= vzMC;
1563     fNMcPrimAccept = nChargedPrim;
1564     fNMcPrimAcceptTracklet = nPseudoTracklets;
1565
1566     if(step==1){ // step 1 = Trig All Mc Nmc
1567         fNmcNch->Fill( fNMcPrimAccept,fNRecAccept);
1568         fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1569         fNmcNchTracklet->Fill( fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1570         fPNmcNchTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1571     }
1572     if(step==3){ // step 3 = Trig vtx Mc
1573         fNmcNchVtx->Fill( fNMcPrimAccept,fNRecAccept);
1574         fNmcNchVtxStrangeCorr->Fill( fNMcPrimAccept,fNRecAcceptStrangeCorr);
1575         fPNmcNchVtx->Fill(fNMcPrimAccept,fNRecAccept);
1576         fNmcNchVtxTracklet->Fill( fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1577         fPNmcNchVtxTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1578     }
1579     return fNRecAccept; // rec value from step 5 or step 2
1580
1581
1582 }
1583
1584 //________________________________________________________________________
1585 void AliAnalysisTaskMinijet::Analyse(const vector<Float_t> &pt,
1586                                      const vector<Float_t> &eta,
1587                                      const vector<Float_t> &phi,
1588                                      const vector<Short_t> &charge,
1589                                      const vector<Float_t> &strangeWeight,
1590                                      const Double_t ntracksCharged,
1591                                      const Int_t ntracklets,
1592                                      const Int_t nAll,
1593                                      const Int_t step)
1594 {
1595
1596     // analyse track properties (comming from either ESDs or AODs) in order to compute
1597     // mini jet activity (chared tracks) as function of charged multiplicity
1598
1599     fStep->Fill(step);
1600
1601     if(fDebug){
1602         Printf("Analysis Step=%d", step);
1603         if(fDebug>2){
1604             Printf("nAll=%d",nAll);
1605             Printf("nCharged=%f",ntracksCharged);
1606             for (Int_t i = 0; i < nAll; i++) {
1607                 Printf("pt[%d]=%f",i,pt[i]);
1608             }
1609         }
1610     }
1611
1612     //calculation of mean pt for all tracks in case of step==0
1613     if(step==5 || step ==2){ // rec step
1614         Double_t meanPt=0.;
1615         Double_t leadingPt=0.;
1616         for (Int_t i = 0; i < nAll; i++) {
1617             if(pt[i]<0.01)continue;
1618             meanPt+=pt[i];
1619             if(leadingPt<pt[i])leadingPt=pt[i];
1620         }
1621         meanPt=meanPt/nAll;
1622         fMeanPtRec=meanPt;
1623         fLeadingPtRec=leadingPt;
1624     }
1625
1626     Double_t propEvent[4] = {ntracksCharged,fVzEvent,fMeanPtRec, fLeadingPtRec}; //vz: {rec, mc, mc}, meanPt and Nrec is always rec value
1627     fMapEvent[step]->Fill(propEvent);
1628
1629     fNcharge[step]->Fill(ntracksCharged);
1630
1631     Float_t ptEventAxis=0;  // pt event axis
1632     Float_t etaEventAxis=0; // eta event axis
1633     Float_t phiEventAxis=0; // phi event axis
1634     Short_t chargeEventAxis=0; // charge event axis
1635     Float_t strangeWeightEventAxis=0;  // strange weight event axis
1636
1637     Float_t ptOthers  = 0; // pt others // for all other tracks around event axis -> see loop
1638     Float_t etaOthers = 0; // eta others
1639     Float_t phiOthers = 0; // phi others
1640     Short_t chargeOthers = 0; // charge others
1641     Float_t strangeWeightOthers  = 0; // strange weight others
1642
1643     Int_t   *pindexInnerEta  = new Int_t  [nAll+1];
1644     Float_t *ptInnerEta      = new Float_t[nAll+1];
1645
1646
1647
1648     for (Int_t i = 0; i < nAll; i++) {
1649
1650         if(pt[i]<0.01)continue;
1651
1652         //fill single particle correction for first step of pair correction
1653         Double_t propAll[3] = {eta[i],pt[i],ntracksCharged };
1654         fMapAll[step]->Fill(propAll, strangeWeight[i]);
1655
1656
1657         //filling of simple check plots
1658         if(pt[i]<0.7) continue;
1659         fPt[step]    -> Fill( pt[i]);
1660         fEta[step]   -> Fill(eta[i]);
1661         fPhi[step]   -> Fill(phi[i]);
1662         fPhiEta[step]-> Fill(phi[i], eta[i]);
1663
1664         pindexInnerEta[i]=0; //set all values to zero
1665         //fill new array for eta check
1666         ptInnerEta[i]= pt[i];
1667
1668     }
1669
1670
1671
1672     // define event axis: leading or random track (with pt>fTriggerPtCut)
1673     // ---------------------------------------
1674     Int_t highPtTracks=0;
1675     Int_t highPtTracksInnerEta=0;
1676     // Double_t highPtTracksInnerEtaCorr=0;
1677     Int_t mult09=0;
1678
1679     //count high pt tracks and high pt tracks in acceptance for seeds
1680     for(Int_t i=0;i<nAll;i++){
1681
1682         if(pt[i]<0.01)continue;
1683
1684         if(TMath::Abs(eta[i])<0.9){
1685             mult09++;
1686         }
1687
1688         if(pt[i]>fTriggerPtCut) {
1689             highPtTracks++;
1690         }
1691
1692         // seed should be place in middle of acceptance, that complete cone is in acceptance
1693         if(pt[i]>fTriggerPtCut && TMath::Abs(eta[i])<fEtaCutSeed && charge[i]!=0){
1694
1695             // Printf("eta=%f", eta[i]);
1696             highPtTracksInnerEta++;
1697
1698         }
1699         else{
1700             ptInnerEta[i]=0;
1701         }
1702     }
1703
1704
1705     //sort array in order to get highest pt tracks first
1706     //index can be computed with pindexInnerEta[number]
1707     if(nAll) TMath::Sort(nAll, ptInnerEta, pindexInnerEta, kTRUE);
1708
1709     //     plot of multiplicity distributions
1710     fNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1711     fPNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1712
1713     if(ntracklets){
1714         fNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1715         fNchTracklet[step]->Fill(ntracklets, ntracksCharged);
1716         fPNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1717     }
1718
1719     //analysis can only be performed with event axis, defined by high pt track
1720
1721
1722     if(highPtTracks>0 && highPtTracksInnerEta>0){
1723
1724         // build pairs in two track loops
1725         // loop over all possible trigger particles (defined by pt_trig and eta_acceptance)
1726         for(Int_t axis=0;(axis<nAll) && (pt[pindexInnerEta[axis]]>=fTriggerPtCut); axis++){
1727
1728             //EventAxisRandom track properties
1729             ptEventAxis  = pt [pindexInnerEta[axis]];
1730             etaEventAxis = eta[pindexInnerEta[axis]];
1731             phiEventAxis = phi[pindexInnerEta[axis]];
1732             chargeEventAxis = charge[pindexInnerEta[axis]];
1733             strangeWeightEventAxis = strangeWeight[pindexInnerEta[axis]];
1734             fPtSeed[step]    -> Fill( ptEventAxis);
1735             fEtaSeed[step]   -> Fill(etaEventAxis);
1736             fPhiSeed[step]   -> Fill(phiEventAxis);
1737
1738
1739             Double_t prop[3] = {etaEventAxis,ptEventAxis,ntracksCharged };
1740             fMapSingleTrig[step]->Fill(prop, strangeWeightEventAxis);
1741
1742             //associated tracks
1743             for (Int_t iTrack = axis+1; iTrack < nAll; iTrack++) {
1744
1745                 if(pt[pindexInnerEta[iTrack]]<fAssociatePtCut) continue;
1746
1747                 if(fSelectParticlesAssoc==1){
1748                     if(charge[pindexInnerEta[iTrack]]==0)continue;
1749                 }
1750                 if(fSelectParticlesAssoc==2){
1751                     if(charge[pindexInnerEta[iTrack]]!=0)continue;
1752                 }
1753
1754
1755                 ptOthers   = pt [pindexInnerEta[iTrack]];
1756                 etaOthers  = eta[pindexInnerEta[iTrack]];
1757                 phiOthers  = phi[pindexInnerEta[iTrack]];
1758                 chargeOthers = charge[pindexInnerEta[iTrack]];
1759                 strangeWeightOthers = strangeWeight[pindexInnerEta[iTrack]];
1760
1761
1762                 //plot only properties of tracks with pt>ptassoc
1763                 fPtOthers[step]    -> Fill( ptOthers);
1764                 fEtaOthers[step]   -> Fill(etaOthers);
1765                 fPhiOthers[step]   -> Fill(phiOthers);
1766                 fPtEtaOthers[step]   -> Fill(ptOthers, etaOthers);
1767
1768                 //      if(fDebug>2)Printf("%f, %f", pt[pindexInnerEta[axis]], pt[pindexInnerEta[iTrack]]);
1769
1770                 Float_t dPhi = phiOthers-phiEventAxis;
1771                 if(dPhi>1.5*TMath::Pi()) dPhi = dPhi-2*TMath::Pi();
1772                 else if(dPhi<-0.5*TMath::Pi())dPhi=dPhi+2*TMath::Pi();
1773                 Float_t dEta=etaOthers-etaEventAxis;
1774
1775
1776                 fDPhiDEtaEventAxis[step]->Fill(dPhi, dEta, strangeWeightEventAxis*strangeWeightOthers);
1777                 fDPhiEventAxis[step]->Fill(dPhi, strangeWeightEventAxis*strangeWeightOthers);
1778
1779                 //check outliers
1780                 if(ptEventAxis< 0.4 || ptEventAxis > 100) AliInfo("particles out of range pt");
1781                 if(ntracksCharged<-1 || ntracksCharged>1500) AliInfo("particles out of range ncharge");
1782                 if(TMath::Abs(dEta)>2*fEtaCut) {
1783                     AliInfo("particles out of range dEta");
1784                     AliInfo(Form("eta1=%f, eta2=%f", etaOthers, etaEventAxis));
1785                     AliInfo(Form("step=%d",step));
1786                 }
1787                 if(dPhi<-0.5*TMath::Pi() || dPhi>1.5*TMath::Pi()){
1788                     AliInfo(Form("particles out of range dPhi"));
1789                     AliInfo(Form("phi1=%f, phi2=%f", phiOthers, phiEventAxis));
1790                 }
1791
1792                 Bool_t isLikeSign = CheckLikeSign(chargeEventAxis, chargeOthers);
1793
1794                 Double_t prop6[6] = {ptEventAxis,ptOthers,dEta,dPhi,ntracksCharged, isLikeSign };
1795                 fMapPair[step]->Fill(prop6, strangeWeightEventAxis*strangeWeightOthers);
1796
1797                 //thrid track loop (Andreas: three particle auto-correlations)
1798                 if(fThreeParticleCorr){
1799                     for (Int_t third = iTrack+1; third < nAll; third++) {
1800                         if(pt[pindexInnerEta[iTrack]]<fTriggerPtCut) continue;
1801                         if(pt[pindexInnerEta[third]]<fTriggerPtCut) continue;
1802
1803                         //dphi1
1804                         Float_t dPhi1 = phiEventAxis - phiOthers;
1805                         if(dPhi1>1.5*TMath::Pi()) dPhi1 = dPhi1-2*TMath::Pi();
1806                         else if(dPhi1<-0.5*TMath::Pi())dPhi1=dPhi1+2*TMath::Pi();
1807
1808                         Float_t phiThird = phi[pindexInnerEta[third]];
1809                         Float_t strangeWeightThird = strangeWeight[pindexInnerEta[third]];
1810
1811                         //dphi2
1812                         Float_t dPhi2 = phiEventAxis - phiThird;
1813                         if(dPhi2>1.5*TMath::Pi()) dPhi2 = dPhi2-2*TMath::Pi();
1814                         else if(dPhi2<-0.5*TMath::Pi())dPhi2=dPhi2+2*TMath::Pi();
1815
1816                         fDPhi1DPhi2[step]-> Fill(dPhi1, dPhi2, strangeWeightEventAxis*strangeWeightOthers*strangeWeightThird);
1817                         Double_t propThree[3] = {dPhi1,dPhi2,ntracksCharged};
1818                         fMapThree[step]->Fill(propThree,strangeWeightEventAxis*strangeWeightOthers*strangeWeightThird );
1819
1820
1821                     }// end of three particle correlation loop
1822
1823                 }// if fThreeParticleCorr is set to true
1824
1825             }// end of inner track loop
1826
1827         } //end of outer track loop
1828
1829         fTriggerNch[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1830         fTriggerNchSeeds[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1831         fTriggerTracklet[step]->Fill(ntracklets);
1832
1833
1834     }//if there is at least one high pt track
1835
1836
1837     if(pindexInnerEta){// clean up array memory used for TMath::Sort
1838         delete[] pindexInnerEta;
1839         pindexInnerEta=0;
1840     }
1841
1842     if(ptInnerEta){// clean up array memory used for TMath::Sort
1843         delete[] ptInnerEta;
1844         ptInnerEta=0;
1845     }
1846
1847 }
1848
1849
1850
1851 //________________________________________________________________________
1852 void AliAnalysisTaskMinijet::Terminate(Option_t*)
1853 {
1854     //terminate function is called at the end
1855     //can be used to draw histograms etc.
1856
1857
1858 }
1859
1860 //________________________________________________________________________
1861 Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(const Short_t charge, const Int_t pdg, Bool_t prim)
1862 {
1863     //selection of mc particle
1864     //fSelectParticles=0: use charged primaries and pi0 and k0
1865     //fSelectParticles=1: use only charged primaries
1866     //fSelectParticles=2: use only pi0 and k0
1867
1868     if(fSelectParticles==0 || fSelectParticles==2){ // in case of 2: need to count also charged particles
1869         if(charge==0){
1870             if(!(pdg==111||pdg==130||pdg==310))
1871                 return false;
1872         }
1873         else{// charge !=0
1874             if(!prim)
1875                 return false;
1876         }
1877     }
1878
1879     else if(fSelectParticles==1){
1880         if (charge==0 || !prim){
1881             return false;
1882         }
1883     }
1884
1885     else{
1886         AliInfo("Error: wrong selection of charged/pi0/k0");
1887         return 0;
1888     }
1889
1890     return true;
1891
1892 }
1893
1894 //________________________________________________________________________
1895 Bool_t AliAnalysisTaskMinijet::SelectParticle(const Short_t charge, const Int_t pdg, const Bool_t prim)
1896 {
1897     //selection of mc particle
1898     //fSelectParticles=0: use charged primaries and pi0 and k0
1899     //fSelectParticles=1: use only charged primaries
1900     //fSelectParticles=2: use only pi0 and k0
1901
1902     if(fSelectParticles==0){
1903         if(charge==0){
1904             if(!(pdg==111||pdg==130||pdg==310))
1905                 return false;
1906         }
1907         else{// charge !=0
1908             if(!prim)
1909                 return false;
1910         }
1911     }
1912     else if (fSelectParticles==1){
1913         if (charge==0 || !prim){
1914             return false;
1915         }
1916     }
1917     else if(fSelectParticles==2){
1918         if(!(pdg==111||pdg==130||pdg==310))
1919             return false;
1920     }
1921
1922     return true;
1923
1924 }
1925
1926 //________________________________________________________________________
1927 Bool_t AliAnalysisTaskMinijet::CheckEvent(const Bool_t recVertex)
1928 {
1929     // This function tests the quality of an event (ESD/AOD) (rec and/or mc part)
1930     // recVertex=false:  check if Mc events and stack is there, Nmc>0
1931     // recVertex=false: " + check if there is a good, reconstructed SPD vertex
1932     // defined by |z|<fVertexCut(10cm), Contributer>0, no PileUpFromSPD(3,0,8)
1933
1934     if(fMode==0){//esd
1935
1936         //mc
1937         if(fUseMC){
1938
1939             //mc event
1940             AliMCEvent *mcEvente = (AliMCEvent*) MCEvent();
1941             if (!mcEvente) {
1942                 Error("CheckEvent", "Could not retrieve MC event");
1943                 return false;
1944             }
1945
1946             //stack
1947             AliStack* stackg = MCEvent()->Stack();
1948             if(!stackg) return false;
1949             Int_t ntracksg = mcEvente->GetNumberOfTracks();
1950             if(ntracksg<0) return false;
1951
1952             //vertex
1953             AliGenEventHeader*  headerg = MCEvent()->GenEventHeader();
1954             TArrayF mcVg;
1955             headerg->PrimaryVertex(mcVg);
1956             //if(TMath::Abs(mcVg[0])<1e-8 && TMath::Abs(mcVg[0])<1e-8 &&
1957             //   TMath::Abs(mcVg[0])<1e-8) return false;
1958             Float_t vzMCg = mcVg[2];
1959             if(TMath::Abs(vzMCg)>fVertexZCut) return false;
1960             //hasVtxMc=true;
1961         }
1962
1963         //rec
1964         if(recVertex==true){
1965             if(fESDEvent->IsPileupFromSPD(3,0,8))return false;
1966
1967             //rec vertex
1968             const AliESDVertex* vertexESDg   = fESDEvent->GetPrimaryVertex(); // uses track or SPD vertexer
1969             if(!vertexESDg) return false;
1970             fVertexCheck->Fill(vertexESDg->GetNContributors());
1971             if(vertexESDg->GetNContributors()<=0)return false;
1972             Float_t fVzg= vertexESDg->GetZ();
1973             if(TMath::Abs(fVzg)>fVertexZCut) return false;
1974
1975             //rec spd vertex
1976             const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
1977             if(!vtxSPD) return false;
1978             if(vtxSPD->GetNContributors()<=0)return false;
1979             Float_t fVzSPD= vtxSPD->GetZ();
1980             if(TMath::Abs(fVzSPD)>fVertexZCut) return false;
1981
1982         }
1983         return true;
1984     }
1985
1986
1987     else if(fMode==1){ //aod
1988
1989         if(fUseMC){
1990
1991             //retreive MC particles from event
1992             TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1993             FindListObject(AliAODMCParticle::StdBranchName());
1994             if(!mcArray){
1995                 AliInfo("No MC particle branch found");
1996                 return false;
1997             }
1998
1999             //mc
2000             AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
2001             FindListObject(AliAODMCHeader::StdBranchName());
2002             Float_t vzMC = aodMCheader->GetVtxZ();
2003             if(TMath::Abs(vzMC)>fVertexZCut) return false;
2004
2005             //hasVtxMc=true;
2006         }
2007
2008         //rec
2009         if(recVertex==true){
2010             //if(fAODEvent->IsPileupFromSPD(3,0.8)) return false;
2011
2012             AliAODVertex*       vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertex();
2013             if(!vertex) return false;
2014             TString vtxTitle(vertex->GetTitle());// only allow vertex from tracks, no vertexer z
2015             // Printf("vtxTitle: %s",vtxTitle.Data());
2016             //if (!(vtxTitle.Contains("VertexerTracksWithConstraint"))) return false;
2017             fVertexCheck->Fill(vertex->GetNContributors());
2018             if(vertex->GetNContributors()<=0) return false;
2019             Double_t vzAOD=vertex->GetZ();
2020             // if(TMath::Abs(vzAOD)<1e-9) return false;
2021             if(TMath::Abs(vzAOD)>fVertexZCut) return false;
2022
2023             AliAODVertex*       vertexSPD= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();
2024             if(!vertexSPD) return false;
2025             if(vertexSPD->GetNContributors()<=0) return false;
2026             Double_t vzSPD=vertexSPD->GetZ();
2027             //if(TMath::Abs(vzSPD)<1e-9) return false;
2028             if(TMath::Abs(vzSPD)>fVertexZCut) return false;
2029
2030             //check TPC reconstruction: check for corrupted chunks
2031             //better: check TPCvertex, but this is not available yet in AODs
2032             Int_t nAcceptedTracksTPC=0;
2033             Int_t nAcceptedTracksITSTPC=0;
2034             for (Int_t iTracks = 0; iTracks < fAODEvent->GetNumberOfTracks(); iTracks++) {
2035                 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
2036                 if (!track) continue;
2037                 if(track->TestFilterBit(128) && TMath::Abs(track->Eta())<fEtaCut &&
2038                    track->Pt()>fPtMin && track->Pt()<fPtMax)
2039                     nAcceptedTracksTPC++;
2040                 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut &&
2041                    track->Pt()>fPtMin && track->Pt()<fPtMax)
2042                     nAcceptedTracksITSTPC++;
2043             }
2044             fCorruptedChunks->Fill(nAcceptedTracksTPC,nAcceptedTracksITSTPC);
2045             if(fRejectChunks){
2046                 if(nAcceptedTracksTPC>fNTPC && nAcceptedTracksITSTPC==0)
2047                     return false;//most likely corrupted chunk. No ITS tracks are reconstructed
2048             }
2049             fCorruptedChunksAfter->Fill(nAcceptedTracksTPC,nAcceptedTracksITSTPC);
2050
2051             //control histograms=================
2052             //tracklet loop
2053             Int_t ntrackletsAccept=0;
2054             AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
2055             for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
2056                 if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 &&
2057                    TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut) ++ntrackletsAccept;
2058             }
2059             Int_t nAcceptedTracks=0;
2060             for (Int_t iTracks = 0; iTracks < fAODEvent->GetNumberOfTracks(); iTracks++) {
2061                 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
2062                 if (!track) continue;
2063                 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut
2064                    && track->Pt()>fPtMin && track->Pt()<fPtMax) nAcceptedTracks++;
2065             }
2066             fNContrNtracklets->Fill(ntrackletsAccept,vertexSPD->GetNContributors());
2067             fNContrNtracks->Fill(nAcceptedTracks,vertexSPD->GetNContributors());
2068             //====================================
2069         }
2070         return true;
2071
2072     }
2073
2074     else {
2075         Printf("Wrong mode!");
2076         return false;
2077     }
2078
2079 }
2080
2081 //_____________________________________________________________________________
2082 const Double_t * AliAnalysisTaskMinijet::CreateLogAxis(const Int_t nbins,
2083                                                        const Double_t xmin,
2084                                                        const Double_t xmax)
2085 {
2086     // returns pointer to an array with can be used to build a logarithmic axis
2087     // it is user responsibility to delete the array
2088
2089     Double_t logxmin = TMath::Log10(xmin);
2090     Double_t logxmax = TMath::Log10(xmax);
2091     Double_t binwidth = (logxmax-logxmin)/nbins;
2092
2093     Double_t *xbins =  new Double_t[nbins+1];
2094
2095     xbins[0] = xmin;
2096     for (Int_t i=1;i<=nbins;i++) {
2097         xbins[i] = xmin + TMath::Power(10,logxmin+i*binwidth);
2098     }
2099
2100     return xbins;
2101 }
2102
2103 //_____________________________________________________________________________
2104 Bool_t AliAnalysisTaskMinijet::CheckLikeSign(const Short_t chargeEventAxis,
2105                                              const Short_t chargeOthers)
2106 {
2107     // compute if charge of two particles/tracks has same sign or different sign
2108
2109     if(fDebug>2) Printf("Charge1=%d, Charge2=%d",chargeEventAxis,chargeOthers);
2110
2111     if(chargeEventAxis<0){
2112         if(chargeOthers<0){
2113             return true;
2114         }
2115         else if(chargeOthers>0){
2116             return false;
2117         }
2118     }
2119
2120     else if(chargeEventAxis>0){
2121         if(chargeOthers>0){
2122             return true;
2123         }
2124         else if(chargeOthers<0){
2125             return false;
2126         }
2127     }
2128
2129     else{
2130         AliInfo("Error: Charge not lower nor higher as zero");
2131         return false;
2132     }
2133
2134     AliInfo("Error: Check values of Charge");
2135     return false;
2136 }
2137
2138