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 = 105;
 180         minbinCentr=0;
 181         maxbinCentr=105;
 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
1240         AliAODTrack* clone = (AliAODTrack*) track->Clone("trk_clone"); //need clone, in order not to change track parameters
1241         Bool_t isDca= clone->PropagateToDCA(fAODEvent->GetPrimaryVertex(),fAODEvent->GetMagneticField(),9999.,d0rphiz,covd0);
1242         delete clone;
1243         fPropagateDca->Fill(Int_t(isDca));
1244         if(TMath::Abs(save - track->Pt())>1e-6) Printf("Before pt=%f, After pt=%f",save, track->Pt());
1245
1246         if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut
1247            && track->Pt()>fPtMin && track->Pt()<fPtMax){
1248
1249             nAcceptedTracks++;
1250
1251             // Printf("dca= %f", track->DCA());
1252             //save track properties in vector
1253             ptArray.push_back(track->Pt());
1254             etaArray.push_back(track->Eta());
1255             phiArray.push_back(track->Phi());
1256             chargeArray.push_back(track->Charge());
1257             fHistPt->Fill(track->Pt());
1258
1259
1260             //correction for underestimation of strangeness in Monte Carlos -> underestimation of contamination
1261             Float_t factor=1.;
1262             if(fUseMC && fCorrStrangeness && step==7){
1263                 if(vtrack->GetLabel()>0){
1264                     AliAODMCParticle* mcprong =(AliAODMCParticle*)mcArray->At(vtrack->GetLabel());
1265                     if(mcprong){
1266                         Int_t labmom = mcprong->GetMother();
1267                         if(labmom>=0){
1268                             AliAODMCParticle* mcmother=(AliAODMCParticle*)mcArray->At(labmom);
1269                             Int_t pdgMother=0;
1270                             if(mcmother) {
1271                                 pdgMother = mcmother->GetPdgCode();
1272                                 if(TMath::Abs(pdgMother)==310 || TMath::Abs(pdgMother)==130 || TMath::Abs(pdgMother)==321){ //K^0_S, K^0_L, K^+-
1273                                     if(track->Pt()<=1)factor=1./0.7; // values from strangeness publication
1274                                     else factor=1./0.6;// values from strangeness publication
1275                                 }
1276                                 if(TMath::Abs(pdgMother)==3122) { //Lambda
1277                                     factor=1./0.25; // values from strangeness publication
1278                                 }
1279                             }
1280                         }
1281                     }
1282                 }
1283             }
1284             nAcceptedTracksStrange+=factor;
1285             strangeArray.push_back(factor);
1286             fDcaXY[step]->Fill(d0rphiz[0], factor);
1287             fDcaZ[step]->Fill(d0rphiz[0], factor);
1288
1289         }
1290     }
1291     //need to check this option for MC
1292     if(nAcceptedTracks==0) return 0;
1293
1294
1295     //tracklet loop
1296     Int_t ntrackletsAccept=0;
1297     AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1298     for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1299         if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut){
1300             ++ntrackletsAccept;
1301         }
1302     }
1303
1304
1305     nTracksTracklets.push_back(nAcceptedTracks);
1306     nTracksTracklets.push_back(ntrackletsAccept);
1307     nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1308     //where here also neutral particles are counted.
1309
1310
1311     fVzEvent= vzAOD;
1312     if(step==5 || step==2){
1313         fNRecAccept = nAcceptedTracks; // needed for MC case //step5 = TrigVtxRecNrec
1314         fNRecAcceptTracklet = ntrackletsAccept; // needed for MC case //step5 = TrigVtxRecNrec
1315     }
1316     if(step==7)fNRecAcceptStrangeCorr = nAcceptedTracksStrange;
1317
1318     return fNRecAccept; // at the moment, always return reconstructed multiplicity
1319
1320 }
1321
1322 //________________________________________________________________________
1323 Double_t AliAnalysisTaskMinijet::ReadEventAODRecMcProp( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1324                                                     vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1325                                                     vector<Float_t> &strangeArray,
1326                                                     vector<Double_t> &nTracksTracklets, const Int_t step)
1327 {
1328     // gives back the number of AOD tracks and pointer to arrays with track
1329     // properties (pt, eta, phi)
1330
1331
1332     ptArray.clear();
1333     etaArray.clear();
1334     phiArray.clear();
1335     chargeArray.clear();
1336     strangeArray.clear();
1337     nTracksTracklets.clear();
1338
1339
1340     // Retreive the number of tracks for this event.
1341     Int_t ntracks = fAODEvent->GetNumberOfTracks();
1342     if(fDebug>1) Printf("AOD tracks: %d", ntracks);
1343
1344
1345     //get array of mc particles
1346     TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1347     FindListObject(AliAODMCParticle::StdBranchName());
1348     if(!mcArray){
1349         AliInfo("No MC particle branch found");
1350         return kFALSE;
1351     }
1352
1353     AliAODVertex*       vtx= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1354     Double_t vzAOD=vtx->GetZ();
1355     fVertexZ[step]->Fill(vzAOD);
1356
1357     Double_t nAcceptedTracks=0;
1358     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1359         AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
1360
1361         AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
1362
1363         if (!track) {
1364             Error("ReadEventAODRecMcProp", "Could not receive track %d", iTracks);
1365             continue;
1366         }
1367
1368         //use only tracks from primaries
1369         if(fAnalysePrimOnly){
1370             if(vtrack->GetLabel()<=0)continue;
1371             if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
1372         }
1373
1374         if(track->TestFilterBit(fFilterBit) &&  TMath::Abs(track->Eta())<fEtaCut &&
1375            track->Pt()>fPtMin && track->Pt()<fPtMax){
1376
1377             nAcceptedTracks++;
1378             Float_t factor=1.;
1379
1380             //save track properties in vector
1381             if(vtrack->GetLabel()<=0){ //fake tracks before "label<0", but crash in AOD079 // what is the meaning of label 0
1382                 //      Printf("Fake track");
1383                 //      continue;
1384                 ptArray.push_back(track->Pt());
1385                 etaArray.push_back(track->Eta());
1386                 phiArray.push_back(track->Phi());
1387                 chargeArray.push_back(track->Charge());
1388
1389             }
1390             else{//mc properties
1391                 AliAODMCParticle *partOfTrack = (AliAODMCParticle*)mcArray->At(vtrack->GetLabel());
1392                 if(!partOfTrack) Printf("label=%d", vtrack->GetLabel());
1393                 if(partOfTrack){
1394                     ptArray.push_back(partOfTrack->Pt());
1395                     etaArray.push_back(partOfTrack->Eta());
1396                     phiArray.push_back(partOfTrack->Phi());
1397                     chargeArray.push_back(vtrack->Charge());//partOfTrack?
1398
1399                     //correction for underestimation of strangeness in Monte Carlos -> underestimation of contamination
1400                     if(fUseMC && fCorrStrangeness && step==6){
1401                         Int_t labmom = partOfTrack->GetMother();
1402                         if(labmom>=0){
1403                             AliAODMCParticle* mcmother=(AliAODMCParticle*)mcArray->At(labmom);
1404                             Int_t pdgMother=0;
1405                             if(mcmother) {
1406                                 pdgMother = mcmother->GetPdgCode();
1407                                 if(TMath::Abs(pdgMother)==310 || TMath::Abs(pdgMother)==130 || TMath::Abs(pdgMother)==321){ //K^0_S, K^0_L, K^+-
1408                                     if(track->Pt()<=1)factor=1./0.7; // values from strangeness publication
1409                                     else factor=1./0.6;// values from strangeness publication
1410                                 }
1411                                 if(TMath::Abs(pdgMother)==3122) { //Lambda
1412                                     factor=1./0.25; // values from strangeness publication
1413                                 }
1414                             }
1415                         }
1416                     }
1417                 }
1418             }
1419             strangeArray.push_back(factor);
1420
1421         }
1422     }
1423     //need to check this option for MC
1424     if(nAcceptedTracks==0) return 0;
1425
1426     //tracklet loop
1427     Int_t ntrackletsAccept=0;
1428     AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1429     for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1430         if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut ){
1431             ++ntrackletsAccept;
1432         }
1433     }
1434
1435
1436     nTracksTracklets.push_back(nAcceptedTracks);
1437     nTracksTracklets.push_back(ntrackletsAccept);
1438     nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1439     //where here also neutral particles are counted.
1440
1441
1442     //check vertex (mc)
1443     AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1444     FindListObject(AliAODMCHeader::StdBranchName());
1445     Float_t vzMC = aodMCheader->GetVtxZ();
1446
1447     fVzEvent= vzMC;
1448     return fNRecAccept;//this is the rec value from step 5
1449
1450 }
1451
1452
1453
1454 //________________________________________________________________________
1455 Double_t AliAnalysisTaskMinijet::ReadEventAODMC( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1456                                              vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1457                                              vector<Float_t> &strangeArray,
1458                                              vector<Double_t> &nTracksTracklets, const Int_t step)
1459 {
1460     // gives back the number of AOD MC particles and pointer to arrays with particle
1461     // properties (pt, eta, phi)
1462
1463     ptArray.clear();
1464     etaArray.clear();
1465     phiArray.clear();
1466     chargeArray.clear();
1467     strangeArray.clear();
1468     nTracksTracklets.clear();
1469
1470     //check vertex
1471     AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1472     FindListObject(AliAODMCHeader::StdBranchName());
1473     Float_t vzMC = aodMCheader->GetVtxZ();
1474     if(step==1){
1475         fVertexZ[0]->Fill(vzMC);
1476     }
1477     fVertexZ[step]->Fill(vzMC);
1478
1479     //retreive MC particles from event
1480     TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1481     FindListObject(AliAODMCParticle::StdBranchName());
1482     if(!mcArray){
1483         AliInfo("No MC particle branch found");
1484         return kFALSE;
1485     }
1486
1487     Int_t ntracks = mcArray->GetEntriesFast();
1488     if(fDebug>1)Printf("MC particles: %d", ntracks);
1489
1490
1491     // Track loop: chek how many particles will be accepted
1492     //Float_t vzMC=0.;
1493     Double_t nChargedPrim=0;
1494     Double_t nAllPrim=0;
1495     Double_t nPseudoTracklets=0;
1496     for (Int_t it = 0; it < ntracks; it++) {
1497         AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1498         if (!track) {
1499             Error("ReadEventAODMC", "Could not receive particle %d", it);
1500             continue;
1501         }
1502
1503         if(!SelectParticlePlusCharged(
1504                                       track->Charge(),
1505                                       track->PdgCode(),
1506                                       track->IsPhysicalPrimary()
1507                                       )
1508            )
1509             continue;
1510
1511         if(TMath::Abs(track->Eta())<fEtaCut && track->Pt()>0.0)++nPseudoTracklets; //0.035
1512         if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1513
1514         nAllPrim++;
1515         if(track->Charge()!=0) nChargedPrim++;
1516
1517     }
1518
1519
1520     if(nAllPrim==0) return 0;
1521     if(nChargedPrim==0) return 0;
1522
1523     //generate array with size of number of accepted tracks
1524     fChargedPi0->Fill(nAllPrim,nChargedPrim);
1525
1526
1527     // Track loop: fill arrays for accepted tracks
1528     // Int_t iChargedPrim=0;
1529     for (Int_t it = 0; it < ntracks; it++) {
1530         AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1531         if (!track) {
1532             Error("ReadEventAODMC", "Could not receive particle %d", it);
1533             continue;
1534         }
1535
1536         if(!SelectParticle(
1537                            track->Charge(),
1538                            track->PdgCode(),
1539                            track->IsPhysicalPrimary()
1540                            )
1541            )
1542             continue;
1543
1544         if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1545
1546         fHistPtMC->Fill(track->Pt());
1547         ptArray.push_back(track->Pt());
1548         etaArray.push_back(track->Eta());
1549         phiArray.push_back(track->Phi());
1550         chargeArray.push_back(track->Charge());
1551         strangeArray.push_back(1);
1552     }
1553
1554     nTracksTracklets.push_back(nChargedPrim);
1555     nTracksTracklets.push_back(nPseudoTracklets);
1556     if(fSelectParticles!=2){
1557         nTracksTracklets.push_back(nAllPrim);
1558     }
1559     else{
1560         nTracksTracklets.push_back(nAllPrim-nChargedPrim); // only neutral
1561     }
1562
1563
1564
1565     fVzEvent= vzMC;
1566     fNMcPrimAccept = nChargedPrim;
1567     fNMcPrimAcceptTracklet = nPseudoTracklets;
1568
1569     if(step==1){ // step 1 = Trig All Mc Nmc
1570         fNmcNch->Fill( fNMcPrimAccept,fNRecAccept);
1571         fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1572         fNmcNchTracklet->Fill( fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1573         fPNmcNchTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1574     }
1575     if(step==3){ // step 3 = Trig vtx Mc
1576         fNmcNchVtx->Fill( fNMcPrimAccept,fNRecAccept);
1577         fNmcNchVtxStrangeCorr->Fill( fNMcPrimAccept,fNRecAcceptStrangeCorr);
1578         fPNmcNchVtx->Fill(fNMcPrimAccept,fNRecAccept);
1579         fNmcNchVtxTracklet->Fill( fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1580         fPNmcNchVtxTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1581     }
1582     return fNRecAccept; // rec value from step 5 or step 2
1583
1584
1585 }
1586
1587 //________________________________________________________________________
1588 void AliAnalysisTaskMinijet::Analyse(const vector<Float_t> &pt,
1589                                      const vector<Float_t> &eta,
1590                                      const vector<Float_t> &phi,
1591                                      const vector<Short_t> &charge,
1592                                      const vector<Float_t> &strangeWeight,
1593                                      const Double_t ntracksCharged,
1594                                      const Int_t ntracklets,
1595                                      const Int_t nAll,
1596                                      const Int_t step)
1597 {
1598
1599     // analyse track properties (comming from either ESDs or AODs) in order to compute
1600     // mini jet activity (chared tracks) as function of charged multiplicity
1601
1602     fStep->Fill(step);
1603
1604     if(fDebug){
1605         Printf("Analysis Step=%d", step);
1606         if(fDebug>2){
1607             Printf("nAll=%d",nAll);
1608             Printf("nCharged=%f",ntracksCharged);
1609             for (Int_t i = 0; i < nAll; i++) {
1610                 Printf("pt[%d]=%f",i,pt[i]);
1611             }
1612         }
1613     }
1614
1615     //calculation of mean pt for all tracks in case of step==0
1616     if(step==5 || step ==2){ // rec step
1617         Double_t meanPt=0.;
1618         Double_t leadingPt=0.;
1619         for (Int_t i = 0; i < nAll; i++) {
1620             if(pt[i]<0.01)continue;
1621             meanPt+=pt[i];
1622             if(leadingPt<pt[i])leadingPt=pt[i];
1623         }
1624         meanPt=meanPt/nAll;
1625         fMeanPtRec=meanPt;
1626         fLeadingPtRec=leadingPt;
1627     }
1628
1629     Double_t propEvent[4] = {ntracksCharged,fVzEvent,fMeanPtRec, fLeadingPtRec}; //vz: {rec, mc, mc}, meanPt and Nrec is always rec value
1630     fMapEvent[step]->Fill(propEvent);
1631
1632     fNcharge[step]->Fill(ntracksCharged);
1633
1634     Float_t ptEventAxis=0;  // pt event axis
1635     Float_t etaEventAxis=0; // eta event axis
1636     Float_t phiEventAxis=0; // phi event axis
1637     Short_t chargeEventAxis=0; // charge event axis
1638     Float_t strangeWeightEventAxis=0;  // strange weight event axis
1639
1640     Float_t ptOthers  = 0; // pt others // for all other tracks around event axis -> see loop
1641     Float_t etaOthers = 0; // eta others
1642     Float_t phiOthers = 0; // phi others
1643     Short_t chargeOthers = 0; // charge others
1644     Float_t strangeWeightOthers  = 0; // strange weight others
1645
1646     Int_t   *pindexInnerEta  = new Int_t  [nAll+1];
1647     Float_t *ptInnerEta      = new Float_t[nAll+1];
1648
1649
1650
1651     for (Int_t i = 0; i < nAll; i++) {
1652
1653         if(pt[i]<0.01)continue;
1654
1655         //fill single particle correction for first step of pair correction
1656         Double_t propAll[3] = {eta[i],pt[i],ntracksCharged };
1657         fMapAll[step]->Fill(propAll, strangeWeight[i]);
1658
1659
1660         //filling of simple check plots
1661         if(pt[i]<0.7) continue;
1662         fPt[step]    -> Fill( pt[i]);
1663         fEta[step]   -> Fill(eta[i]);
1664         fPhi[step]   -> Fill(phi[i]);
1665         fPhiEta[step]-> Fill(phi[i], eta[i]);
1666
1667         pindexInnerEta[i]=0; //set all values to zero
1668         //fill new array for eta check
1669         ptInnerEta[i]= pt[i];
1670
1671     }
1672
1673
1674
1675     // define event axis: leading or random track (with pt>fTriggerPtCut)
1676     // ---------------------------------------
1677     Int_t highPtTracks=0;
1678     Int_t highPtTracksInnerEta=0;
1679     // Double_t highPtTracksInnerEtaCorr=0;
1680     Int_t mult09=0;
1681
1682     //count high pt tracks and high pt tracks in acceptance for seeds
1683     for(Int_t i=0;i<nAll;i++){
1684
1685         if(pt[i]<0.01)continue;
1686
1687         if(TMath::Abs(eta[i])<0.9){
1688             mult09++;
1689         }
1690
1691         if(pt[i]>fTriggerPtCut) {
1692             highPtTracks++;
1693         }
1694
1695         // seed should be place in middle of acceptance, that complete cone is in acceptance
1696         if(pt[i]>fTriggerPtCut && TMath::Abs(eta[i])<fEtaCutSeed && charge[i]!=0){
1697
1698             // Printf("eta=%f", eta[i]);
1699             highPtTracksInnerEta++;
1700
1701         }
1702         else{
1703             ptInnerEta[i]=0;
1704         }
1705     }
1706
1707
1708     //sort array in order to get highest pt tracks first
1709     //index can be computed with pindexInnerEta[number]
1710     if(nAll) TMath::Sort(nAll, ptInnerEta, pindexInnerEta, kTRUE);
1711
1712     //     plot of multiplicity distributions
1713     fNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1714     fPNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1715
1716     if(ntracklets){
1717         fNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1718         fNchTracklet[step]->Fill(ntracklets, ntracksCharged);
1719         fPNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1720     }
1721
1722     //analysis can only be performed with event axis, defined by high pt track
1723
1724
1725     if(highPtTracks>0 && highPtTracksInnerEta>0){
1726
1727         // build pairs in two track loops
1728         // loop over all possible trigger particles (defined by pt_trig and eta_acceptance)
1729         for(Int_t axis=0;(axis<nAll) && (pt[pindexInnerEta[axis]]>=fTriggerPtCut); axis++){
1730
1731             //EventAxisRandom track properties
1732             ptEventAxis  = pt [pindexInnerEta[axis]];
1733             etaEventAxis = eta[pindexInnerEta[axis]];
1734             phiEventAxis = phi[pindexInnerEta[axis]];
1735             chargeEventAxis = charge[pindexInnerEta[axis]];
1736             strangeWeightEventAxis = strangeWeight[pindexInnerEta[axis]];
1737             fPtSeed[step]    -> Fill( ptEventAxis);
1738             fEtaSeed[step]   -> Fill(etaEventAxis);
1739             fPhiSeed[step]   -> Fill(phiEventAxis);
1740
1741
1742             Double_t prop[3] = {etaEventAxis,ptEventAxis,ntracksCharged };
1743             fMapSingleTrig[step]->Fill(prop, strangeWeightEventAxis);
1744
1745             //associated tracks
1746             for (Int_t iTrack = axis+1; iTrack < nAll; iTrack++) {
1747
1748                 if(pt[pindexInnerEta[iTrack]]<fAssociatePtCut) continue;
1749
1750                 if(fSelectParticlesAssoc==1){
1751                     if(charge[pindexInnerEta[iTrack]]==0)continue;
1752                 }
1753                 if(fSelectParticlesAssoc==2){
1754                     if(charge[pindexInnerEta[iTrack]]!=0)continue;
1755                 }
1756
1757
1758                 ptOthers   = pt [pindexInnerEta[iTrack]];
1759                 etaOthers  = eta[pindexInnerEta[iTrack]];
1760                 phiOthers  = phi[pindexInnerEta[iTrack]];
1761                 chargeOthers = charge[pindexInnerEta[iTrack]];
1762                 strangeWeightOthers = strangeWeight[pindexInnerEta[iTrack]];
1763
1764
1765                 //plot only properties of tracks with pt>ptassoc
1766                 fPtOthers[step]    -> Fill( ptOthers);
1767                 fEtaOthers[step]   -> Fill(etaOthers);
1768                 fPhiOthers[step]   -> Fill(phiOthers);
1769                 fPtEtaOthers[step]   -> Fill(ptOthers, etaOthers);
1770
1771                 //      if(fDebug>2)Printf("%f, %f", pt[pindexInnerEta[axis]], pt[pindexInnerEta[iTrack]]);
1772
1773                 Float_t dPhi = phiOthers-phiEventAxis;
1774                 if(dPhi>1.5*TMath::Pi()) dPhi = dPhi-2*TMath::Pi();
1775                 else if(dPhi<-0.5*TMath::Pi())dPhi=dPhi+2*TMath::Pi();
1776                 Float_t dEta=etaOthers-etaEventAxis;
1777
1778
1779                 fDPhiDEtaEventAxis[step]->Fill(dPhi, dEta, strangeWeightEventAxis*strangeWeightOthers);
1780                 fDPhiEventAxis[step]->Fill(dPhi, strangeWeightEventAxis*strangeWeightOthers);
1781
1782                 //check outliers
1783                 if(ptEventAxis< 0.4 || ptEventAxis > 100) AliInfo("particles out of range pt");
1784                 if(ntracksCharged<-1 || ntracksCharged>1500) AliInfo("particles out of range ncharge");
1785                 if(TMath::Abs(dEta)>2*fEtaCut) {
1786                     AliInfo("particles out of range dEta");
1787                     AliInfo(Form("eta1=%f, eta2=%f", etaOthers, etaEventAxis));
1788                     AliInfo(Form("step=%d",step));
1789                 }
1790                 if(dPhi<-0.5*TMath::Pi() || dPhi>1.5*TMath::Pi()){
1791                     AliInfo(Form("particles out of range dPhi"));
1792                     AliInfo(Form("phi1=%f, phi2=%f", phiOthers, phiEventAxis));
1793                 }
1794
1795                 Bool_t isLikeSign = CheckLikeSign(chargeEventAxis, chargeOthers);
1796
1797                 Double_t prop6[6] = {ptEventAxis,ptOthers,dEta,dPhi,ntracksCharged, static_cast<Double_t>(isLikeSign) };
1798                 fMapPair[step]->Fill(prop6, strangeWeightEventAxis*strangeWeightOthers);
1799
1800                 //thrid track loop (Andreas: three particle auto-correlations)
1801                 if(fThreeParticleCorr){
1802                     for (Int_t third = iTrack+1; third < nAll; third++) {
1803                         if(pt[pindexInnerEta[iTrack]]<fTriggerPtCut) continue;
1804                         if(pt[pindexInnerEta[third]]<fTriggerPtCut) continue;
1805
1806                         //dphi1
1807                         Float_t dPhi1 = phiEventAxis - phiOthers;
1808                         if(dPhi1>1.5*TMath::Pi()) dPhi1 = dPhi1-2*TMath::Pi();
1809                         else if(dPhi1<-0.5*TMath::Pi())dPhi1=dPhi1+2*TMath::Pi();
1810
1811                         Float_t phiThird = phi[pindexInnerEta[third]];
1812                         Float_t strangeWeightThird = strangeWeight[pindexInnerEta[third]];
1813
1814                         //dphi2
1815                         Float_t dPhi2 = phiEventAxis - phiThird;
1816                         if(dPhi2>1.5*TMath::Pi()) dPhi2 = dPhi2-2*TMath::Pi();
1817                         else if(dPhi2<-0.5*TMath::Pi())dPhi2=dPhi2+2*TMath::Pi();
1818
1819                         fDPhi1DPhi2[step]-> Fill(dPhi1, dPhi2, strangeWeightEventAxis*strangeWeightOthers*strangeWeightThird);
1820                         Double_t propThree[3] = {dPhi1,dPhi2,ntracksCharged};
1821                         fMapThree[step]->Fill(propThree,strangeWeightEventAxis*strangeWeightOthers*strangeWeightThird );
1822
1823
1824                     }// end of three particle correlation loop
1825
1826                 }// if fThreeParticleCorr is set to true
1827
1828             }// end of inner track loop
1829
1830         } //end of outer track loop
1831
1832         fTriggerNch[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1833         fTriggerNchSeeds[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1834         fTriggerTracklet[step]->Fill(ntracklets);
1835
1836
1837     }//if there is at least one high pt track
1838
1839
1840     if(pindexInnerEta){// clean up array memory used for TMath::Sort
1841         delete[] pindexInnerEta;
1842         pindexInnerEta=0;
1843     }
1844
1845     if(ptInnerEta){// clean up array memory used for TMath::Sort
1846         delete[] ptInnerEta;
1847         ptInnerEta=0;
1848     }
1849
1850 }
1851
1852
1853
1854 //________________________________________________________________________
1855 void AliAnalysisTaskMinijet::Terminate(Option_t*)
1856 {
1857     //terminate function is called at the end
1858     //can be used to draw histograms etc.
1859
1860
1861 }
1862
1863 //________________________________________________________________________
1864 Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(const Short_t charge, const Int_t pdg, Bool_t prim)
1865 {
1866     //selection of mc particle
1867     //fSelectParticles=0: use charged primaries and pi0 and k0
1868     //fSelectParticles=1: use only charged primaries
1869     //fSelectParticles=2: use only pi0 and k0
1870
1871     if(fSelectParticles==0 || fSelectParticles==2){ // in case of 2: need to count also charged particles
1872         if(charge==0){
1873             if(!(pdg==111||pdg==130||pdg==310))
1874                 return false;
1875         }
1876         else{// charge !=0
1877             if(!prim)
1878                 return false;
1879         }
1880     }
1881
1882     else if(fSelectParticles==1){
1883         if (charge==0 || !prim){
1884             return false;
1885         }
1886     }
1887
1888     else{
1889         AliInfo("Error: wrong selection of charged/pi0/k0");
1890         return 0;
1891     }
1892
1893     return true;
1894
1895 }
1896
1897 //________________________________________________________________________
1898 Bool_t AliAnalysisTaskMinijet::SelectParticle(const Short_t charge, const Int_t pdg, const Bool_t prim)
1899 {
1900     //selection of mc particle
1901     //fSelectParticles=0: use charged primaries and pi0 and k0
1902     //fSelectParticles=1: use only charged primaries
1903     //fSelectParticles=2: use only pi0 and k0
1904
1905     if(fSelectParticles==0){
1906         if(charge==0){
1907             if(!(pdg==111||pdg==130||pdg==310))
1908                 return false;
1909         }
1910         else{// charge !=0
1911             if(!prim)
1912                 return false;
1913         }
1914     }
1915     else if (fSelectParticles==1){
1916         if (charge==0 || !prim){
1917             return false;
1918         }
1919     }
1920     else if(fSelectParticles==2){
1921         if(!(pdg==111||pdg==130||pdg==310))
1922             return false;
1923     }
1924
1925     return true;
1926
1927 }
1928
1929 //________________________________________________________________________
1930 Bool_t AliAnalysisTaskMinijet::CheckEvent(const Bool_t recVertex)
1931 {
1932     // This function tests the quality of an event (ESD/AOD) (rec and/or mc part)
1933     // recVertex=false:  check if Mc events and stack is there, Nmc>0
1934     // recVertex=false: " + check if there is a good, reconstructed SPD vertex
1935     // defined by |z|<fVertexCut(10cm), Contributer>0, no PileUpFromSPD(3,0,8)
1936
1937     if(fMode==0){//esd
1938
1939         //mc
1940         if(fUseMC){
1941
1942             //mc event
1943             AliMCEvent *mcEvente = (AliMCEvent*) MCEvent();
1944             if (!mcEvente) {
1945                 Error("CheckEvent", "Could not retrieve MC event");
1946                 return false;
1947             }
1948
1949             //stack
1950             AliStack* stackg = MCEvent()->Stack();
1951             if(!stackg) return false;
1952             Int_t ntracksg = mcEvente->GetNumberOfTracks();
1953             if(ntracksg<0) return false;
1954
1955             //vertex
1956             AliGenEventHeader*  headerg = MCEvent()->GenEventHeader();
1957             TArrayF mcVg;
1958             headerg->PrimaryVertex(mcVg);
1959             //if(TMath::Abs(mcVg[0])<1e-8 && TMath::Abs(mcVg[0])<1e-8 &&
1960             //   TMath::Abs(mcVg[0])<1e-8) return false;
1961             Float_t vzMCg = mcVg[2];
1962             if(TMath::Abs(vzMCg)>fVertexZCut) return false;
1963             //hasVtxMc=true;
1964         }
1965
1966         //rec
1967         if(recVertex==true){
1968             if(fESDEvent->IsPileupFromSPD(3,0,8))return false;
1969
1970             //rec vertex
1971             const AliESDVertex* vertexESDg   = fESDEvent->GetPrimaryVertex(); // uses track or SPD vertexer
1972             if(!vertexESDg) return false;
1973             fVertexCheck->Fill(vertexESDg->GetNContributors());
1974             if(vertexESDg->GetNContributors()<=0)return false;
1975             Float_t fVzg= vertexESDg->GetZ();
1976             if(TMath::Abs(fVzg)>fVertexZCut) return false;
1977
1978             //rec spd vertex
1979             const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
1980             if(!vtxSPD) return false;
1981             if(vtxSPD->GetNContributors()<=0)return false;
1982             Float_t fVzSPD= vtxSPD->GetZ();
1983             if(TMath::Abs(fVzSPD)>fVertexZCut) return false;
1984
1985         }
1986         return true;
1987     }
1988
1989
1990     else if(fMode==1){ //aod
1991
1992         if(fUseMC){
1993
1994             //retreive MC particles from event
1995             TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1996             FindListObject(AliAODMCParticle::StdBranchName());
1997             if(!mcArray){
1998                 AliInfo("No MC particle branch found");
1999                 return false;
2000             }
2001
2002             //mc
2003             AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
2004             FindListObject(AliAODMCHeader::StdBranchName());
2005             Float_t vzMC = aodMCheader->GetVtxZ();
2006             if(TMath::Abs(vzMC)>fVertexZCut) return false;
2007
2008             //hasVtxMc=true;
2009         }
2010
2011         //rec
2012         if(recVertex==true){
2013             //if(fAODEvent->IsPileupFromSPD(3,0.8)) return false;
2014
2015             AliAODVertex*       vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertex();
2016             if(!vertex) return false;
2017             TString vtxTitle(vertex->GetTitle());// only allow vertex from tracks, no vertexer z
2018             // Printf("vtxTitle: %s",vtxTitle.Data());
2019             //if (!(vtxTitle.Contains("VertexerTracksWithConstraint"))) return false;
2020             fVertexCheck->Fill(vertex->GetNContributors());
2021             if(vertex->GetNContributors()<=0) return false;
2022             Double_t vzAOD=vertex->GetZ();
2023             // if(TMath::Abs(vzAOD)<1e-9) return false;
2024             if(TMath::Abs(vzAOD)>fVertexZCut) return false;
2025
2026             AliAODVertex*       vertexSPD= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();
2027             if(!vertexSPD) return false;
2028             if(vertexSPD->GetNContributors()<=0) return false;
2029             Double_t vzSPD=vertexSPD->GetZ();
2030             //if(TMath::Abs(vzSPD)<1e-9) return false;
2031             if(TMath::Abs(vzSPD)>fVertexZCut) return false;
2032
2033             //check TPC reconstruction: check for corrupted chunks
2034             //better: check TPCvertex, but this is not available yet in AODs
2035             Int_t nAcceptedTracksTPC=0;
2036             Int_t nAcceptedTracksITSTPC=0;
2037             for (Int_t iTracks = 0; iTracks < fAODEvent->GetNumberOfTracks(); iTracks++) {
2038                 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
2039                 if (!track) continue;
2040                 if(track->TestFilterBit(128) && TMath::Abs(track->Eta())<fEtaCut &&
2041                    track->Pt()>fPtMin && track->Pt()<fPtMax)
2042                     nAcceptedTracksTPC++;
2043                 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut &&
2044                    track->Pt()>fPtMin && track->Pt()<fPtMax)
2045                     nAcceptedTracksITSTPC++;
2046             }
2047             fCorruptedChunks->Fill(nAcceptedTracksTPC,nAcceptedTracksITSTPC);
2048             if(fRejectChunks){
2049                 if(nAcceptedTracksTPC>fNTPC && nAcceptedTracksITSTPC==0)
2050                     return false;//most likely corrupted chunk. No ITS tracks are reconstructed
2051             }
2052             fCorruptedChunksAfter->Fill(nAcceptedTracksTPC,nAcceptedTracksITSTPC);
2053
2054             //control histograms=================
2055             //tracklet loop
2056             Int_t ntrackletsAccept=0;
2057             AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
2058             for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
2059                 if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 &&
2060                    TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut) ++ntrackletsAccept;
2061             }
2062             Int_t nAcceptedTracks=0;
2063             for (Int_t iTracks = 0; iTracks < fAODEvent->GetNumberOfTracks(); iTracks++) {
2064                 AliAODTrack *track = (AliAODTrack *)fAODEvent->GetTrack(iTracks);
2065                 if (!track) continue;
2066                 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut
2067                    && track->Pt()>fPtMin && track->Pt()<fPtMax) nAcceptedTracks++;
2068             }
2069             fNContrNtracklets->Fill(ntrackletsAccept,vertexSPD->GetNContributors());
2070             fNContrNtracks->Fill(nAcceptedTracks,vertexSPD->GetNContributors());
2071             //====================================
2072         }
2073         return true;
2074
2075     }
2076
2077     else {
2078         Printf("Wrong mode!");
2079         return false;
2080     }
2081
2082 }
2083
2084 //_____________________________________________________________________________
2085 const Double_t * AliAnalysisTaskMinijet::CreateLogAxis(const Int_t nbins,
2086                                                        const Double_t xmin,
2087                                                        const Double_t xmax)
2088 {
2089     // returns pointer to an array with can be used to build a logarithmic axis
2090     // it is user responsibility to delete the array
2091
2092     Double_t logxmin = TMath::Log10(xmin);
2093     Double_t logxmax = TMath::Log10(xmax);
2094     Double_t binwidth = (logxmax-logxmin)/nbins;
2095
2096     Double_t *xbins =  new Double_t[nbins+1];
2097
2098     xbins[0] = xmin;
2099     for (Int_t i=1;i<=nbins;i++) {
2100         xbins[i] = xmin + TMath::Power(10,logxmin+i*binwidth);
2101     }
2102
2103     return xbins;
2104 }
2105
2106 //_____________________________________________________________________________
2107 Bool_t AliAnalysisTaskMinijet::CheckLikeSign(const Short_t chargeEventAxis,
2108                                              const Short_t chargeOthers)
2109 {
2110     // compute if charge of two particles/tracks has same sign or different sign
2111
2112     if(fDebug>2) Printf("Charge1=%d, Charge2=%d",chargeEventAxis,chargeOthers);
2113
2114     if(chargeEventAxis<0){
2115         if(chargeOthers<0){
2116             return true;
2117         }
2118         else if(chargeOthers>0){
2119             return false;
2120         }
2121     }
2122
2123     else if(chargeEventAxis>0){
2124         if(chargeOthers>0){
2125             return true;
2126         }
2127         else if(chargeOthers<0){
2128             return false;
2129         }
2130     }
2131
2132     else{
2133         AliInfo("Error: Charge not lower nor higher as zero");
2134         return false;
2135     }
2136
2137     AliInfo("Error: Check values of Charge");
2138     return false;
2139 }
2140
2141