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 = ((AliVAODHeader*)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->GetNumberOfTracks();
 592                 for(Int_t itrack=0; itrack<nTracks; itrack++) {
 593                     AliAODTrack * track = dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(itrack));
 594                     if(!track) AliFatal("Not a standard AOD");
 595                     if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
 596                         useEvent=true;
 597                         break;
 598                     }
 599                 }
 600                 if (!useEvent) return;
 601             }
 602             else if(fSelOption==1){
 603                 TString trcl = fAODEvent->GetFiredTriggerClasses().Data();
 604                 if (!(trcl.Contains("CINT1-B-NOPF-ALLNOTRD"))) return;
 605             }
 606         }
 607     }
 608
 609     //reset values
 610     fNMcPrimAccept=0;// number of accepted primaries
 611     fNRecAccept=0;   // number of accepted tracks
 612     fNRecAcceptStrangeCorr=0;   // number of accepted tracks + strangeness correction
 613     fNMcPrimAcceptTracklet=0;// number of accepted primaries (no pt cut)
 614     fNRecAcceptTracklet=0;   // number of accepted tracklets
 615
 616     // instead of task->SelectCollisionCandidate(mask) in AddTask macro
 617     Bool_t isSelected = (((AliInputEventHandler*)
 618                           (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 619                          ->IsEventSelected() &  fTriggerType);
 620
 621
 622     if(fDebug){
 623         Printf("IsSelected = %d", isSelected);
 624         Printf("CheckEvent(true)= %d", CheckEvent(true));
 625         Printf("CheckEvent(false)= %d", CheckEvent(false));
 626     }
 627
 628     fEventStat->Fill(0);//all events
 629
 630     //check trigger
 631     if(isSelected){ // has offline trigger
 632
 633         fEventStat->Fill(1);//triggered event
 634
 635         if(CheckEvent(true)){//step 5 = TrigVtxRecNrec, step 4 = TrigVtxRecMcPropNrec ,step 3 = TrigVtxMcNrec
 636
 637             fEventStat->Fill(2);//triggered event with vertex
 638
 639             if(!fMcOnly){
 640                 //step 5 = TrigVtxRecNrec
 641
 642                 // read tracks
 643                 if(fESDEvent)     ntracks = ReadEventESD(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 5);
 644                 else if(fAODEvent)ntracks = ReadEventAOD(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 5);
 645                 else AliInfo("Fatal Error");
 646
 647                 if (fCentralityMethod.Length() > 0)
 648                     ntracks = centrality;
 649
 650                 // analyse
 651                 if(pt.size()){ //(internally ntracks=fNRecAccept)
 652                     fEventStat->Fill(3);//triggered event with vertex and one reconstructed track in acceptance
 653                         Analyse(pt, eta, phi, charge, strangenessWeight, ntracks, nTracksTracklets[1], nTracksTracklets[2], 5);//step 5 = TrigVtxRecNrec
 654                     }
 655
 656                 if(fCorrStrangeness){
 657                     //step 7 = TrigVtxRecNrecStrangeCorr
 658
 659                     // read tracks
 660                     if(fESDEvent)     ntracks = ReadEventESD(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 7);//stagness version not yet implemented
 661                     else if(fAODEvent)ntracks = ReadEventAOD(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 7);
 662                     else AliInfo("Fatal Error");
 663
 664                     if (fCentralityMethod.Length() > 0)
 665                         ntracks = centrality;
 666
 667                     // analyse
 668                     if(pt.size()){ //(internally ntracks=fNRecAccept)
 669                         Analyse(pt, eta, phi, charge, strangenessWeight, fNRecAcceptStrangeCorr, nTracksTracklets[1], nTracksTracklets[2], 7);//step 7 = TrigVtxRecNrecStrangeCorr
 670                     }
 671                 }
 672
 673                 if(fUseMC){
 674                     // step 4 = TrigVtxRecMcPropNrec
 675
 676                     // read tracks
 677                     if(fESDEvent)       ntracks = ReadEventESDRecMcProp(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 4);
 678                     else if(fAODEvent)  ntracks = ReadEventAODRecMcProp(pt, eta, phi, charge,strangenessWeight, nTracksTracklets, 4);
 679                     else AliInfo("Fatal Error");
 680
 681                     if (fCentralityMethod.Length() > 0)
 682                         ntracks = centrality;
 683
 684                     //analyse
 685                     if(pt.size()){//(internally ntracks=fNRecAccept)
 686                         Analyse(pt, eta, phi, charge,strangenessWeight, ntracks, nTracksTracklets[1], nTracksTracklets[2], 4); //step 4 = TrigVtxRecMcPropNrec
 687                     }
 688
 689
 690                     if(fCorrStrangeness){
 691                         // step 6 = TrigVtxRecMcPropNrecStrangeCorr
 692
 693                         // read tracks
 694                         if(fESDEvent)       ntracks = ReadEventESDRecMcProp(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 6);//stagness version not yet implemented
 695                         else if(fAODEvent)  ntracks = ReadEventAODRecMcProp(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 6);
 696                         else AliInfo("Fatal Error");
 697
 698                         if (fCentralityMethod.Length() > 0)
 699                             ntracks = centrality;
 700
 701                         //analyse
 702                         if(pt.size()){//(internally ntracks=fNRecAccept)
 703                             Analyse(pt, eta, phi, charge, strangenessWeight, fNRecAcceptStrangeCorr, nTracksTracklets[1], nTracksTracklets[2], 6); //step 6 = TrigVtxRecMcPropNrecStrangeCorr
 704                         }
 705                     }
 706                     // step 3 = TrigVtxMcNrec
 707
 708                     // read tracks
 709                     if(fESDEvent)       ntracks = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 3);
 710                     else if(fAODEvent)  ntracks = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 3);
 711                     else AliInfo("Fatal Error");
 712
 713                     if (fCentralityMethod.Length() > 0){
 714                         fNRecAccept = centrality;
 715                         fNMcPrimAccept = centrality;
 716                     }
 717
 718                     // analyse
 719                     if(pt.size()){
 720                         Analyse(pt, eta, phi, charge, strangenessWeight, fNRecAccept,    nTracksTracklets[1],nTracksTracklets[2], 3); //step 3 = TrigVtxMcNrec
 721                         Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 2); //step 2 = TrigVtxMcNmc
 722                     }
 723
 724                 }
 725             }
 726
 727         }//check event (true)
 728
 729
 730         if(fUseMC && !fMcOnly){
 731             //reset values
 732             fNMcPrimAccept=0;// number of accepted primaries
 733             fNRecAccept=0;   // number of accepted tracks
 734             fNMcPrimAcceptTracklet=0;// number of accepted primaries (no pt cut)
 735             fNRecAcceptTracklet=0;   // number of accepted tracklets
 736
 737             if(CheckEvent(false)){// all events, with and without reconstucted vertex
 738                 if(fESDEvent) ntracks       = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 1);//read tracks
 739                 else if(fAODEvent) ntracks  = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 1);//read tracks
 740                 else AliInfo("Fatal Error");
 741
 742                 if (fCentralityMethod.Length() > 0)
 743                     fNMcPrimAccept = centrality;
 744
 745
 746                 // analyse
 747                 if(pt.size()){
 748                     Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 1);  // step 1 = TrigAllMcNmc
 749
 750                     Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);  //first part of step 0 // step 0 = AllAllMcNmc
 751                 }
 752
 753
 754             }
 755         }
 756     }// triggered event
 757
 758     else { // not selected by physics selection task = not triggered
 759         if(fUseMC && !fMcOnly){
 760             if(CheckEvent(false)){
 761
 762                 //read tracks
 763                 if(fESDEvent)      ntracks  = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
 764                 else if(fAODEvent) ntracks  = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
 765                 else AliInfo("Fatal Error");
 766
 767                 if (fCentralityMethod.Length() > 0)
 768                     fNMcPrimAccept = centrality;
 769
 770                 //analyse
 771                 if(pt.size()){
 772                     Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);  //second part of step 0 // step 0 = AllAllMcNmc
 773                 }
 774             }
 775         }
 776     }
 777
 778     if(fMcOnly){
 779         // read event
 780         if(fMode==0)       ntracks  = ReadEventESDMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
 781         else if (fMode==1) ntracks  = ReadEventAODMC(pt, eta, phi, charge, strangenessWeight, nTracksTracklets, 0);
 782
 783         if (fCentralityMethod.Length() > 0)
 784             fNMcPrimAccept = centrality;
 785
 786         // analyse
 787         if(pt.size()){
 788             Analyse(pt, eta, phi, charge, strangenessWeight, fNMcPrimAccept, nTracksTracklets[1],nTracksTracklets[2], 0);
 789         }
 790     }
 791 }
 792
 793
 794 //________________________________________________________________________
 795 Double_t AliAnalysisTaskMinijet::ReadEventESD( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 796                                            vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
 797                                            vector<Float_t> &strangeArray,
 798                                            vector<Double_t>   &nTracksTracklets, const Int_t step)
 799 {
 800     // gives back the number of esd tracks and pointer to arrays with track
 801     // properties (pt, eta, phi)
 802     // Uses TPC tracks with SPD vertex from now on
 803
 804     ptArray.clear();
 805     etaArray.clear();
 806     phiArray.clear();
 807     chargeArray.clear();
 808     strangeArray.clear();
 809     nTracksTracklets.clear();
 810
 811     const AliESDVertex* vtxSPD   = fESDEvent->GetPrimaryVertexSPD(); // uses track or SPD vertexer
 812     fVertexZ[step]->Fill(vtxSPD->GetZ());
 813
 814     // Retreive the number of all tracks for this event.
 815     Int_t ntracks = fESDEvent->GetNumberOfTracks();
 816     if(fDebug>1)  Printf("all ESD tracks: %d", ntracks);
 817
 818     //first loop to check how many tracks are accepted
 819     //------------------
 820     Double_t nAcceptedTracks=0;
 821     //Float_t nAcceptedTracksStrange=0;
 822     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 823
 824         AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
 825         if (!esdTrack) {
 826             Error("ReadEventESD", "Could not receive track %d", iTracks);
 827             continue;
 828         }
 829
 830         // use TPC only tracks with non default SPD vertex
 831         AliESDtrack *track = AliESDtrackCuts::
 832         GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
 833         if(!track) continue;
 834         if(!fCuts->AcceptTrack(track)) {
 835             delete track;
 836             continue;// apply TPC track cuts before constrain to SPD vertex
 837         }
 838         if(track->Pt()>0.){
 839             // only constrain tracks above threshold
 840             AliExternalTrackParam exParam;
 841             // take the B-field from the ESD, no 3D fieldMap available at this point
 842             Bool_t relate = false;
 843             relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
 844                                               kVeryBig,&exParam);
 845             if(!relate){
 846                 delete track;
 847                 continue;
 848             }
 849             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
 850                        exParam.GetCovariance());
 851         }
 852         else{
 853             delete track;
 854             continue;// only if tracks have pt<=0
 855         }
 856
 857         if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax) {
 858             ptArray.push_back(track->Pt());
 859             etaArray.push_back(track->Eta());
 860             phiArray.push_back(track->Phi());
 861             chargeArray.push_back(track->Charge());
 862             strangeArray.push_back(1);
 863             ++nAcceptedTracks;
 864             fHistPt->Fill(track->Pt());
 865         }
 866
 867         // TPC only track memory needs to be cleaned up
 868         if(track)
 869             delete track;
 870
 871     }
 872
 873     //need to be checked
 874     if(nAcceptedTracks==0) return 0;
 875
 876     //tracklet loop
 877     Int_t ntrackletsAccept=0;
 878     AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
 879     Int_t ntracklets = mult->GetNumberOfTracklets();
 880     for(Int_t i=0;i< ntracklets;i++){
 881         if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
 882             ntrackletsAccept++;
 883         }
 884     }
 885     nTracksTracklets.push_back(nAcceptedTracks);
 886     nTracksTracklets.push_back(ntrackletsAccept);
 887     nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
 888     //where here also neutral particles are counted.
 889
 890
 891     fVzEvent=vtxSPD->GetZ(); // needed for correction map
 892     if(step==5 || step ==2) {
 893         fNRecAccept=nAcceptedTracks;
 894         fNRecAcceptTracklet=ntrackletsAccept;
 895     }
 896     return fNRecAccept;
 897
 898
 899 }
 900
 901 //________________________________________________________________________
 902 Double_t AliAnalysisTaskMinijet::ReadEventESDRecMcProp( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
 903                                                     vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
 904                                                     vector<Float_t> &strangeArray,
 905                                                     vector<Double_t> &nTracksTracklets, const Int_t step)
 906 {
 907     // gives back the number of esd tracks and pointer to arrays with track
 908     // properties (pt, eta, phi) of mc particles if available
 909     // Uses TPC tracks with SPD vertex from now on
 910
 911     ptArray.clear();
 912     etaArray.clear();
 913     phiArray.clear();
 914     chargeArray.clear();
 915     strangeArray.clear();
 916     nTracksTracklets.clear();
 917
 918
 919     AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
 920     if (!mcEvent) {
 921         Error("ReadEventESDRecMcProp", "Could not retrieve MC event");
 922         return 0;
 923     }
 924     AliStack* stack = MCEvent()->Stack();
 925     if(!stack) return 0;
 926
 927
 928     // Retreive the number of all tracks for this event.
 929     Int_t ntracks = fESDEvent->GetNumberOfTracks();
 930     if(fDebug>1)Printf("all ESD tracks: %d", ntracks);
 931
 932     const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
 933     fVertexZ[step]->Fill(vtxSPD->GetZ());
 934
 935     //track loop
 936     Double_t nAcceptedTracks=0;
 937     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 938
 939         AliVParticle *vtrack = fESDEvent->GetTrack(iTracks);
 940         AliESDtrack *esdTrack = (AliESDtrack *)fESDEvent->GetTrack(iTracks);
 941         if (!esdTrack) {
 942             Error("ReadEventESDRecMcProp", "Could not receive track %d", iTracks);
 943             continue;
 944         }
 945
 946         // use TPC only tracks with non default SPD vertex
 947         AliESDtrack *track = AliESDtrackCuts::
 948         GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(fESDEvent),esdTrack->GetID());
 949         if(!track) continue;
 950         if(!fCuts->AcceptTrack(track)) {
 951             delete track;
 952             continue;// apply TPC track cuts before constrain to SPD vertex
 953         }
 954         if(track->Pt()>0.){
 955             // only constrain tracks above threshold
 956             AliExternalTrackParam exParam;
 957             // take the B-field from the ESD, no 3D fieldMap available at this point
 958             Bool_t relate = false;
 959             relate = track->RelateToVertexTPC(vtxSPD,fESDEvent->GetMagneticField(),
 960                                               kVeryBig,&exParam);
 961             if(!relate){
 962                 delete track;
 963                 continue;
 964             }
 965             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
 966                        exParam.GetCovariance());
 967         }
 968         else{
 969             delete track;
 970             continue;// only if tracks have pt<=0
 971         }
 972
 973         //count tracks, if available, use mc particle properties
 974         if(vtrack->GetLabel()<=0){
 975             if (TMath::Abs(track->Eta())<fEtaCut && track->Pt()>fPtMin && track->Pt()<fPtMax){
 976                 ptArray.push_back(track->Pt());
 977                 etaArray.push_back(track->Eta());
 978                 phiArray.push_back(track->Phi());
 979                 chargeArray.push_back(track->Charge());
 980                 strangeArray.push_back(1);
 981                 ++nAcceptedTracks;
 982             }
 983         }
 984         else{
 985             TParticle *partOfTrack = stack->Particle(vtrack->GetLabel());
 986             if (TMath::Abs(partOfTrack->Eta())<fEtaCut && partOfTrack->Pt()>fPtMin && partOfTrack->Pt()<fPtMax) {
 987                 ptArray.push_back(partOfTrack->Pt());
 988                 etaArray.push_back(partOfTrack->Eta());
 989                 phiArray.push_back(partOfTrack->Phi());
 990                 chargeArray.push_back(vtrack->Charge());
 991                 strangeArray.push_back(1);
 992                 ++nAcceptedTracks;
 993             }
 994         }
 995
 996         // TPC only track memory needs to be cleaned up
 997         if(track)
 998             delete track;
 999
1000     }
1001
1002     if(nAcceptedTracks==0) return 0;
1003
1004     //tracklet loop
1005     Int_t ntrackletsAccept=0;
1006     AliMultiplicity * mult = (AliMultiplicity*)(fESDEvent->GetMultiplicity());
1007     Int_t ntracklets = mult->GetNumberOfTracklets();
1008     for(Int_t i=0;i< ntracklets;i++){
1009         if(mult->GetDeltaPhi(i)<0.05 && TMath::Abs(mult->GetEta(i))<fEtaCut){
1010             ntrackletsAccept++;
1011         }
1012     }
1013
1014     nTracksTracklets.push_back(nAcceptedTracks);
1015     nTracksTracklets.push_back(ntrackletsAccept);
1016     nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1017     //where here also neutral particles are counted.
1018
1019
1020     //get mc vertex for correction maps
1021     AliGenEventHeader*  header = MCEvent()->GenEventHeader();
1022     TArrayF mcV;
1023     header->PrimaryVertex(mcV);
1024     fVzEvent= mcV[2];
1025
1026     return fNRecAccept; // give back reconstructed value
1027
1028
1029 }
1030
1031
1032
1033
1034 //________________________________________________________________________
1035 Double_t AliAnalysisTaskMinijet::ReadEventESDMC(vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1036                                              vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1037                                              vector<Float_t> &strangeArray,
1038                                              vector<Double_t> &nTracksTracklets, const Int_t step)
1039 {
1040     // gives back the number of charged prim MC particle and pointer to arrays
1041     // with particle properties (pt, eta, phi)
1042
1043     ptArray.clear();
1044     etaArray.clear();
1045     phiArray.clear();
1046     chargeArray.clear();
1047     strangeArray.clear();
1048     nTracksTracklets.clear();
1049
1050     fNMcPrimAccept=0;
1051
1052     AliMCEvent *mcEvent = (AliMCEvent*) MCEvent();
1053     if (!mcEvent) {
1054         Error("ReadEventESDMC", "Could not retrieve MC event");
1055         return 0;
1056     }
1057
1058     AliStack* stack = MCEvent()->Stack();
1059     if(!stack) return 0;
1060
1061     Int_t ntracks = mcEvent->GetNumberOfTracks();
1062     if(fDebug>1)Printf("MC particles: %d", ntracks);
1063
1064     //vertex
1065     AliGenEventHeader*  header = MCEvent()->GenEventHeader();
1066     TArrayF mcV;
1067     Float_t vzMC=0.;
1068     if(header){
1069         header->PrimaryVertex(mcV);
1070         vzMC = mcV[2];
1071         if(step==1){
1072             fVertexZ[0]->Fill(vzMC);
1073         }
1074         fVertexZ[step]->Fill(vzMC);
1075     }
1076
1077     //----------------------------------
1078     //first track loop to check how many chared primary tracks are available
1079     Double_t nChargedPrimaries=0;
1080     Double_t nAllPrimaries=0;
1081
1082     Double_t nPseudoTracklets=0;
1083     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1084         AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
1085         if (!track) {
1086             Error("ReadEventESDMC", "Could not receive track %d", iTracks);
1087             continue;
1088         }
1089
1090
1091         if(//count also charged particles in case of fSelectParticles==2 (only neutral)
1092            !SelectParticlePlusCharged(
1093                                       track->Charge(),
1094                                       track->PdgCode(),
1095                                       stack->IsPhysicalPrimary(track->Label())
1096                                       )
1097            )
1098             continue;
1099
1100         //count number of pseudo tracklets
1101         if(TMath::Abs(track->Eta())<=fEtaCut && track->Pt()>0.0) ++nPseudoTracklets; //0.035
1102         //same cuts as on ESDtracks
1103         if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1104
1105         //count all primaries
1106         ++nAllPrimaries;
1107         //count charged primaries
1108         if (track->Charge() != 0) ++nChargedPrimaries;
1109
1110         if(fDebug>2) Printf("PDG=%d, IsPrim=%d",  track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
1111
1112
1113     }
1114     //----------------------------------
1115
1116     if(fDebug>2){
1117         Printf("All in acceptance=%f",nAllPrimaries);
1118         Printf("Charged in acceptance =%f",nChargedPrimaries);
1119     }
1120
1121     fChargedPi0->Fill(nAllPrimaries,nChargedPrimaries);
1122
1123     if(nAllPrimaries==0) return 0;
1124     if(nChargedPrimaries==0) return 0;
1125
1126     //track loop
1127     //Int_t iChargedPiK=0;
1128     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1129         AliMCParticle *track = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(iTracks));
1130         if (!track) {
1131             Error("ReadEventESDMC", "Could not receive track %d", iTracks);
1132             continue;
1133         }
1134
1135         if(!SelectParticle(
1136                            track->Charge(),
1137                            track->PdgCode(),
1138                            stack->IsPhysicalPrimary(track->Label())
1139                            )
1140            )
1141             continue;
1142
1143
1144         //same cuts as on ESDtracks
1145         if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin  || track->Pt()>fPtMax) continue;
1146
1147         if(fDebug>2) Printf("After: PDG=%d, IsPrim=%d",  track->PdgCode(),stack->IsPhysicalPrimary(track->Label()) );
1148
1149
1150         fHistPtMC->Fill(track->Pt());
1151         //fills arrays with track properties
1152         ptArray.push_back(track->Pt());
1153         etaArray.push_back(track->Eta());
1154         phiArray.push_back(track->Phi());
1155         chargeArray.push_back(track->Charge());
1156         strangeArray.push_back(1);
1157
1158     } //track loop
1159
1160     nTracksTracklets.push_back(nChargedPrimaries);
1161     nTracksTracklets.push_back(nPseudoTracklets);
1162     if(fSelectParticles!=2){
1163         nTracksTracklets.push_back(nAllPrimaries);
1164     }
1165     else{
1166         nTracksTracklets.push_back(nAllPrimaries-nChargedPrimaries); // only neutral
1167     }
1168
1169     fNMcPrimAccept = nChargedPrimaries;
1170     fNMcPrimAcceptTracklet = nPseudoTracklets;
1171
1172     if(step==1){
1173         fNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1174         fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1175         fNmcNchTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1176         fPNmcNchTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1177     }
1178     if(step==3){
1179         fNmcNchVtx->Fill(fNMcPrimAccept,fNRecAccept);
1180         fPNmcNchVtx->Fill(fNMcPrimAccept,fNRecAccept);
1181         fNmcNchVtxTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1182         fPNmcNchVtxTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1183     }
1184
1185     fVzEvent= vzMC;
1186     return fNRecAccept;
1187
1188 }
1189
1190 //________________________________________________________________________
1191 Double_t AliAnalysisTaskMinijet::ReadEventAOD( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1192                                            vector<Float_t> &phiArray,  vector<Short_t> &chargeArray,
1193                                            vector<Float_t> &strangeArray,
1194                                            vector<Double_t> &nTracksTracklets, const Int_t step)
1195 {
1196     // gives back the number of AOD tracks and pointer to arrays with track
1197     // properties (pt, eta, phi)
1198
1199     ptArray.clear();
1200     etaArray.clear();
1201     phiArray.clear();
1202     chargeArray.clear();
1203     strangeArray.clear();
1204     nTracksTracklets.clear();
1205
1206     TClonesArray *mcArray=0x0;
1207     if(fAnalysePrimOnly || (fCorrStrangeness && fUseMC)){
1208         mcArray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
1209     }
1210
1211
1212     AliAODVertex*       vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1213     Double_t vzAOD=vertex->GetZ();
1214     fVertexZ[step]->Fill(vzAOD);
1215
1216     // Retreive the number of tracks for this event.
1217     Int_t ntracks = fAODEvent->GetNumberOfTracks();
1218     if(fDebug>1) Printf("AOD tracks: %d", ntracks);
1219
1220
1221     Double_t nAcceptedTracks=0;
1222     Float_t nAcceptedTracksStrange=0;
1223     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1224         AliAODTrack *track = dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(iTracks));
1225         if(!track) AliFatal("Not a standard AOD");
1226         if (!track) {
1227             Error("ReadEventAOD", "Could not receive track %d", iTracks);
1228             continue;
1229         }
1230
1231         AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
1232
1233         //use only tracks from primaries
1234         if(fAnalysePrimOnly){
1235             if(vtrack->GetLabel()<=0)continue;
1236             if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
1237         }
1238
1239         Double_t save= track->Pt();
1240         Double_t d0rphiz[2],covd0[3];
1241
1242         AliAODTrack* clone = (AliAODTrack*) track->Clone("trk_clone"); //need clone, in order not to change track parameters
1243         Bool_t isDca= clone->PropagateToDCA(fAODEvent->GetPrimaryVertex(),fAODEvent->GetMagneticField(),9999.,d0rphiz,covd0);
1244         delete clone;
1245         fPropagateDca->Fill(Int_t(isDca));
1246         if(TMath::Abs(save - track->Pt())>1e-6) Printf("Before pt=%f, After pt=%f",save, track->Pt());
1247
1248         if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut
1249            && track->Pt()>fPtMin && track->Pt()<fPtMax){
1250
1251             nAcceptedTracks++;
1252
1253             // Printf("dca= %f", track->DCA());
1254             //save track properties in vector
1255             ptArray.push_back(track->Pt());
1256             etaArray.push_back(track->Eta());
1257             phiArray.push_back(track->Phi());
1258             chargeArray.push_back(track->Charge());
1259             fHistPt->Fill(track->Pt());
1260
1261
1262             //correction for underestimation of strangeness in Monte Carlos -> underestimation of contamination
1263             Float_t factor=1.;
1264             if(fUseMC && fCorrStrangeness && step==7){
1265                 if(vtrack->GetLabel()>0){
1266                     AliAODMCParticle* mcprong =(AliAODMCParticle*)mcArray->At(vtrack->GetLabel());
1267                     if(mcprong){
1268                         Int_t labmom = mcprong->GetMother();
1269                         if(labmom>=0){
1270                             AliAODMCParticle* mcmother=(AliAODMCParticle*)mcArray->At(labmom);
1271                             Int_t pdgMother=0;
1272                             if(mcmother) {
1273                                 pdgMother = mcmother->GetPdgCode();
1274                                 if(TMath::Abs(pdgMother)==310 || TMath::Abs(pdgMother)==130 || TMath::Abs(pdgMother)==321){ //K^0_S, K^0_L, K^+-
1275                                     if(track->Pt()<=1)factor=1./0.7; // values from strangeness publication
1276                                     else factor=1./0.6;// values from strangeness publication
1277                                 }
1278                                 if(TMath::Abs(pdgMother)==3122) { //Lambda
1279                                     factor=1./0.25; // values from strangeness publication
1280                                 }
1281                             }
1282                         }
1283                     }
1284                 }
1285             }
1286             nAcceptedTracksStrange+=factor;
1287             strangeArray.push_back(factor);
1288             fDcaXY[step]->Fill(d0rphiz[0], factor);
1289             fDcaZ[step]->Fill(d0rphiz[0], factor);
1290
1291         }
1292     }
1293     //need to check this option for MC
1294     if(nAcceptedTracks==0) return 0;
1295
1296
1297     //tracklet loop
1298     Int_t ntrackletsAccept=0;
1299     AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1300     for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1301         if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut){
1302             ++ntrackletsAccept;
1303         }
1304     }
1305
1306
1307     nTracksTracklets.push_back(nAcceptedTracks);
1308     nTracksTracklets.push_back(ntrackletsAccept);
1309     nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1310     //where here also neutral particles are counted.
1311
1312
1313     fVzEvent= vzAOD;
1314     if(step==5 || step==2){
1315         fNRecAccept = nAcceptedTracks; // needed for MC case //step5 = TrigVtxRecNrec
1316         fNRecAcceptTracklet = ntrackletsAccept; // needed for MC case //step5 = TrigVtxRecNrec
1317     }
1318     if(step==7)fNRecAcceptStrangeCorr = nAcceptedTracksStrange;
1319
1320     return fNRecAccept; // at the moment, always return reconstructed multiplicity
1321
1322 }
1323
1324 //________________________________________________________________________
1325 Double_t AliAnalysisTaskMinijet::ReadEventAODRecMcProp( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1326                                                     vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1327                                                     vector<Float_t> &strangeArray,
1328                                                     vector<Double_t> &nTracksTracklets, const Int_t step)
1329 {
1330     // gives back the number of AOD tracks and pointer to arrays with track
1331     // properties (pt, eta, phi)
1332
1333
1334     ptArray.clear();
1335     etaArray.clear();
1336     phiArray.clear();
1337     chargeArray.clear();
1338     strangeArray.clear();
1339     nTracksTracklets.clear();
1340
1341
1342     // Retreive the number of tracks for this event.
1343     Int_t ntracks = fAODEvent->GetNumberOfTracks();
1344     if(fDebug>1) Printf("AOD tracks: %d", ntracks);
1345
1346
1347     //get array of mc particles
1348     TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1349     FindListObject(AliAODMCParticle::StdBranchName());
1350     if(!mcArray){
1351         AliInfo("No MC particle branch found");
1352         return kFALSE;
1353     }
1354
1355     AliAODVertex*       vtx= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1356     Double_t vzAOD=vtx->GetZ();
1357     fVertexZ[step]->Fill(vzAOD);
1358
1359     Double_t nAcceptedTracks=0;
1360     for (Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1361         AliAODTrack *track = dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(iTracks));
1362         if(!track) AliFatal("Not a standard AOD");
1363
1364         AliVParticle *vtrack = fAODEvent->GetTrack(iTracks);
1365
1366         if (!track) {
1367             Error("ReadEventAODRecMcProp", "Could not receive track %d", iTracks);
1368             continue;
1369         }
1370
1371         //use only tracks from primaries
1372         if(fAnalysePrimOnly){
1373             if(vtrack->GetLabel()<=0)continue;
1374             if(!(static_cast<AliAODMCParticle*>(mcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))continue;
1375         }
1376
1377         if(track->TestFilterBit(fFilterBit) &&  TMath::Abs(track->Eta())<fEtaCut &&
1378            track->Pt()>fPtMin && track->Pt()<fPtMax){
1379
1380             nAcceptedTracks++;
1381             Float_t factor=1.;
1382
1383             //save track properties in vector
1384             if(vtrack->GetLabel()<=0){ //fake tracks before "label<0", but crash in AOD079 // what is the meaning of label 0
1385                 //      Printf("Fake track");
1386                 //      continue;
1387                 ptArray.push_back(track->Pt());
1388                 etaArray.push_back(track->Eta());
1389                 phiArray.push_back(track->Phi());
1390                 chargeArray.push_back(track->Charge());
1391
1392             }
1393             else{//mc properties
1394                 AliAODMCParticle *partOfTrack = (AliAODMCParticle*)mcArray->At(vtrack->GetLabel());
1395                 if(!partOfTrack) Printf("label=%d", vtrack->GetLabel());
1396                 if(partOfTrack){
1397                     ptArray.push_back(partOfTrack->Pt());
1398                     etaArray.push_back(partOfTrack->Eta());
1399                     phiArray.push_back(partOfTrack->Phi());
1400                     chargeArray.push_back(vtrack->Charge());//partOfTrack?
1401
1402                     //correction for underestimation of strangeness in Monte Carlos -> underestimation of contamination
1403                     if(fUseMC && fCorrStrangeness && step==6){
1404                         Int_t labmom = partOfTrack->GetMother();
1405                         if(labmom>=0){
1406                             AliAODMCParticle* mcmother=(AliAODMCParticle*)mcArray->At(labmom);
1407                             Int_t pdgMother=0;
1408                             if(mcmother) {
1409                                 pdgMother = mcmother->GetPdgCode();
1410                                 if(TMath::Abs(pdgMother)==310 || TMath::Abs(pdgMother)==130 || TMath::Abs(pdgMother)==321){ //K^0_S, K^0_L, K^+-
1411                                     if(track->Pt()<=1)factor=1./0.7; // values from strangeness publication
1412                                     else factor=1./0.6;// values from strangeness publication
1413                                 }
1414                                 if(TMath::Abs(pdgMother)==3122) { //Lambda
1415                                     factor=1./0.25; // values from strangeness publication
1416                                 }
1417                             }
1418                         }
1419                     }
1420                 }
1421             }
1422             strangeArray.push_back(factor);
1423
1424         }
1425     }
1426     //need to check this option for MC
1427     if(nAcceptedTracks==0) return 0;
1428
1429     //tracklet loop
1430     Int_t ntrackletsAccept=0;
1431     AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
1432     for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
1433         if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 && TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut ){
1434             ++ntrackletsAccept;
1435         }
1436     }
1437
1438
1439     nTracksTracklets.push_back(nAcceptedTracks);
1440     nTracksTracklets.push_back(ntrackletsAccept);
1441     nTracksTracklets.push_back(nAcceptedTracks);//in order to be compatible with mc analysis
1442     //where here also neutral particles are counted.
1443
1444
1445     //check vertex (mc)
1446     AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1447     FindListObject(AliAODMCHeader::StdBranchName());
1448     Float_t vzMC = aodMCheader->GetVtxZ();
1449
1450     fVzEvent= vzMC;
1451     return fNRecAccept;//this is the rec value from step 5
1452
1453 }
1454
1455
1456
1457 //________________________________________________________________________
1458 Double_t AliAnalysisTaskMinijet::ReadEventAODMC( vector<Float_t> &ptArray,  vector<Float_t> &etaArray,
1459                                              vector<Float_t> &phiArray, vector<Short_t> &chargeArray,
1460                                              vector<Float_t> &strangeArray,
1461                                              vector<Double_t> &nTracksTracklets, const Int_t step)
1462 {
1463     // gives back the number of AOD MC particles and pointer to arrays with particle
1464     // properties (pt, eta, phi)
1465
1466     ptArray.clear();
1467     etaArray.clear();
1468     phiArray.clear();
1469     chargeArray.clear();
1470     strangeArray.clear();
1471     nTracksTracklets.clear();
1472
1473     //check vertex
1474     AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
1475     FindListObject(AliAODMCHeader::StdBranchName());
1476     Float_t vzMC = aodMCheader->GetVtxZ();
1477     if(step==1){
1478         fVertexZ[0]->Fill(vzMC);
1479     }
1480     fVertexZ[step]->Fill(vzMC);
1481
1482     //retreive MC particles from event
1483     TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1484     FindListObject(AliAODMCParticle::StdBranchName());
1485     if(!mcArray){
1486         AliInfo("No MC particle branch found");
1487         return kFALSE;
1488     }
1489
1490     Int_t ntracks = mcArray->GetEntriesFast();
1491     if(fDebug>1)Printf("MC particles: %d", ntracks);
1492
1493
1494     // Track loop: chek how many particles will be accepted
1495     //Float_t vzMC=0.;
1496     Double_t nChargedPrim=0;
1497     Double_t nAllPrim=0;
1498     Double_t nPseudoTracklets=0;
1499     for (Int_t it = 0; it < ntracks; it++) {
1500         AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1501         if (!track) {
1502             Error("ReadEventAODMC", "Could not receive particle %d", it);
1503             continue;
1504         }
1505
1506         if(!SelectParticlePlusCharged(
1507                                       track->Charge(),
1508                                       track->PdgCode(),
1509                                       track->IsPhysicalPrimary()
1510                                       )
1511            )
1512             continue;
1513
1514         if(TMath::Abs(track->Eta())<fEtaCut && track->Pt()>0.0)++nPseudoTracklets; //0.035
1515         if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1516
1517         nAllPrim++;
1518         if(track->Charge()!=0) nChargedPrim++;
1519
1520     }
1521
1522
1523     if(nAllPrim==0) return 0;
1524     if(nChargedPrim==0) return 0;
1525
1526     //generate array with size of number of accepted tracks
1527     fChargedPi0->Fill(nAllPrim,nChargedPrim);
1528
1529
1530     // Track loop: fill arrays for accepted tracks
1531     // Int_t iChargedPrim=0;
1532     for (Int_t it = 0; it < ntracks; it++) {
1533         AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
1534         if (!track) {
1535             Error("ReadEventAODMC", "Could not receive particle %d", it);
1536             continue;
1537         }
1538
1539         if(!SelectParticle(
1540                            track->Charge(),
1541                            track->PdgCode(),
1542                            track->IsPhysicalPrimary()
1543                            )
1544            )
1545             continue;
1546
1547         if(TMath::Abs(track->Eta())>fEtaCut || track->Pt()<fPtMin || track->Pt()>fPtMax) continue;
1548
1549         fHistPtMC->Fill(track->Pt());
1550         ptArray.push_back(track->Pt());
1551         etaArray.push_back(track->Eta());
1552         phiArray.push_back(track->Phi());
1553         chargeArray.push_back(track->Charge());
1554         strangeArray.push_back(1);
1555     }
1556
1557     nTracksTracklets.push_back(nChargedPrim);
1558     nTracksTracklets.push_back(nPseudoTracklets);
1559     if(fSelectParticles!=2){
1560         nTracksTracklets.push_back(nAllPrim);
1561     }
1562     else{
1563         nTracksTracklets.push_back(nAllPrim-nChargedPrim); // only neutral
1564     }
1565
1566
1567
1568     fVzEvent= vzMC;
1569     fNMcPrimAccept = nChargedPrim;
1570     fNMcPrimAcceptTracklet = nPseudoTracklets;
1571
1572     if(step==1){ // step 1 = Trig All Mc Nmc
1573         fNmcNch->Fill( fNMcPrimAccept,fNRecAccept);
1574         fPNmcNch->Fill(fNMcPrimAccept,fNRecAccept);
1575         fNmcNchTracklet->Fill( fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1576         fPNmcNchTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1577     }
1578     if(step==3){ // step 3 = Trig vtx Mc
1579         fNmcNchVtx->Fill( fNMcPrimAccept,fNRecAccept);
1580         fNmcNchVtxStrangeCorr->Fill( fNMcPrimAccept,fNRecAcceptStrangeCorr);
1581         fPNmcNchVtx->Fill(fNMcPrimAccept,fNRecAccept);
1582         fNmcNchVtxTracklet->Fill( fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1583         fPNmcNchVtxTracklet->Fill(fNMcPrimAcceptTracklet,fNRecAcceptTracklet);
1584     }
1585     return fNRecAccept; // rec value from step 5 or step 2
1586
1587
1588 }
1589
1590 //________________________________________________________________________
1591 void AliAnalysisTaskMinijet::Analyse(const vector<Float_t> &pt,
1592                                      const vector<Float_t> &eta,
1593                                      const vector<Float_t> &phi,
1594                                      const vector<Short_t> &charge,
1595                                      const vector<Float_t> &strangeWeight,
1596                                      const Double_t ntracksCharged,
1597                                      const Int_t ntracklets,
1598                                      const Int_t nAll,
1599                                      const Int_t step)
1600 {
1601
1602     // analyse track properties (comming from either ESDs or AODs) in order to compute
1603     // mini jet activity (chared tracks) as function of charged multiplicity
1604
1605     fStep->Fill(step);
1606
1607     if(fDebug){
1608         Printf("Analysis Step=%d", step);
1609         if(fDebug>2){
1610             Printf("nAll=%d",nAll);
1611             Printf("nCharged=%f",ntracksCharged);
1612             for (Int_t i = 0; i < nAll; i++) {
1613                 Printf("pt[%d]=%f",i,pt[i]);
1614             }
1615         }
1616     }
1617
1618     //calculation of mean pt for all tracks in case of step==0
1619     if(step==5 || step ==2){ // rec step
1620         Double_t meanPt=0.;
1621         Double_t leadingPt=0.;
1622         for (Int_t i = 0; i < nAll; i++) {
1623             if(pt[i]<0.01)continue;
1624             meanPt+=pt[i];
1625             if(leadingPt<pt[i])leadingPt=pt[i];
1626         }
1627         meanPt=meanPt/nAll;
1628         fMeanPtRec=meanPt;
1629         fLeadingPtRec=leadingPt;
1630     }
1631
1632     Double_t propEvent[4] = {ntracksCharged,fVzEvent,fMeanPtRec, fLeadingPtRec}; //vz: {rec, mc, mc}, meanPt and Nrec is always rec value
1633     fMapEvent[step]->Fill(propEvent);
1634
1635     fNcharge[step]->Fill(ntracksCharged);
1636
1637     Float_t ptEventAxis=0;  // pt event axis
1638     Float_t etaEventAxis=0; // eta event axis
1639     Float_t phiEventAxis=0; // phi event axis
1640     Short_t chargeEventAxis=0; // charge event axis
1641     Float_t strangeWeightEventAxis=0;  // strange weight event axis
1642
1643     Float_t ptOthers  = 0; // pt others // for all other tracks around event axis -> see loop
1644     Float_t etaOthers = 0; // eta others
1645     Float_t phiOthers = 0; // phi others
1646     Short_t chargeOthers = 0; // charge others
1647     Float_t strangeWeightOthers  = 0; // strange weight others
1648
1649     Int_t   *pindexInnerEta  = new Int_t  [nAll+1];
1650     Float_t *ptInnerEta      = new Float_t[nAll+1];
1651
1652
1653
1654     for (Int_t i = 0; i < nAll; i++) {
1655
1656         if(pt[i]<0.01)continue;
1657
1658         //fill single particle correction for first step of pair correction
1659         Double_t propAll[3] = {eta[i],pt[i],ntracksCharged };
1660         fMapAll[step]->Fill(propAll, strangeWeight[i]);
1661
1662
1663         //filling of simple check plots
1664         if(pt[i]<0.7) continue;
1665         fPt[step]    -> Fill( pt[i]);
1666         fEta[step]   -> Fill(eta[i]);
1667         fPhi[step]   -> Fill(phi[i]);
1668         fPhiEta[step]-> Fill(phi[i], eta[i]);
1669
1670         pindexInnerEta[i]=0; //set all values to zero
1671         //fill new array for eta check
1672         ptInnerEta[i]= pt[i];
1673
1674     }
1675
1676
1677
1678     // define event axis: leading or random track (with pt>fTriggerPtCut)
1679     // ---------------------------------------
1680     Int_t highPtTracks=0;
1681     Int_t highPtTracksInnerEta=0;
1682     // Double_t highPtTracksInnerEtaCorr=0;
1683     Int_t mult09=0;
1684
1685     //count high pt tracks and high pt tracks in acceptance for seeds
1686     for(Int_t i=0;i<nAll;i++){
1687
1688         if(pt[i]<0.01)continue;
1689
1690         if(TMath::Abs(eta[i])<0.9){
1691             mult09++;
1692         }
1693
1694         if(pt[i]>fTriggerPtCut) {
1695             highPtTracks++;
1696         }
1697
1698         // seed should be place in middle of acceptance, that complete cone is in acceptance
1699         if(pt[i]>fTriggerPtCut && TMath::Abs(eta[i])<fEtaCutSeed && charge[i]!=0){
1700
1701             // Printf("eta=%f", eta[i]);
1702             highPtTracksInnerEta++;
1703
1704         }
1705         else{
1706             ptInnerEta[i]=0;
1707         }
1708     }
1709
1710
1711     //sort array in order to get highest pt tracks first
1712     //index can be computed with pindexInnerEta[number]
1713     if(nAll) TMath::Sort(nAll, ptInnerEta, pindexInnerEta, kTRUE);
1714
1715     //     plot of multiplicity distributions
1716     fNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1717     fPNch07Nch[step]->Fill(ntracksCharged, highPtTracksInnerEta);
1718
1719     if(ntracklets){
1720         fNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1721         fNchTracklet[step]->Fill(ntracklets, ntracksCharged);
1722         fPNch07Tracklet[step]->Fill(ntracklets, highPtTracksInnerEta);//only counts tracks which can be used as seeds
1723     }
1724
1725     //analysis can only be performed with event axis, defined by high pt track
1726
1727
1728     if(highPtTracks>0 && highPtTracksInnerEta>0){
1729
1730         // build pairs in two track loops
1731         // loop over all possible trigger particles (defined by pt_trig and eta_acceptance)
1732         for(Int_t axis=0;(axis<nAll) && (pt[pindexInnerEta[axis]]>=fTriggerPtCut); axis++){
1733
1734             //EventAxisRandom track properties
1735             ptEventAxis  = pt [pindexInnerEta[axis]];
1736             etaEventAxis = eta[pindexInnerEta[axis]];
1737             phiEventAxis = phi[pindexInnerEta[axis]];
1738             chargeEventAxis = charge[pindexInnerEta[axis]];
1739             strangeWeightEventAxis = strangeWeight[pindexInnerEta[axis]];
1740             fPtSeed[step]    -> Fill( ptEventAxis);
1741             fEtaSeed[step]   -> Fill(etaEventAxis);
1742             fPhiSeed[step]   -> Fill(phiEventAxis);
1743
1744
1745             Double_t prop[3] = {etaEventAxis,ptEventAxis,ntracksCharged };
1746             fMapSingleTrig[step]->Fill(prop, strangeWeightEventAxis);
1747
1748             //associated tracks
1749             for (Int_t iTrack = axis+1; iTrack < nAll; iTrack++) {
1750
1751                 if(pt[pindexInnerEta[iTrack]]<fAssociatePtCut) continue;
1752
1753                 if(fSelectParticlesAssoc==1){
1754                     if(charge[pindexInnerEta[iTrack]]==0)continue;
1755                 }
1756                 if(fSelectParticlesAssoc==2){
1757                     if(charge[pindexInnerEta[iTrack]]!=0)continue;
1758                 }
1759
1760
1761                 ptOthers   = pt [pindexInnerEta[iTrack]];
1762                 etaOthers  = eta[pindexInnerEta[iTrack]];
1763                 phiOthers  = phi[pindexInnerEta[iTrack]];
1764                 chargeOthers = charge[pindexInnerEta[iTrack]];
1765                 strangeWeightOthers = strangeWeight[pindexInnerEta[iTrack]];
1766
1767
1768                 //plot only properties of tracks with pt>ptassoc
1769                 fPtOthers[step]    -> Fill( ptOthers);
1770                 fEtaOthers[step]   -> Fill(etaOthers);
1771                 fPhiOthers[step]   -> Fill(phiOthers);
1772                 fPtEtaOthers[step]   -> Fill(ptOthers, etaOthers);
1773
1774                 //      if(fDebug>2)Printf("%f, %f", pt[pindexInnerEta[axis]], pt[pindexInnerEta[iTrack]]);
1775
1776                 Float_t dPhi = phiOthers-phiEventAxis;
1777                 if(dPhi>1.5*TMath::Pi()) dPhi = dPhi-2*TMath::Pi();
1778                 else if(dPhi<-0.5*TMath::Pi())dPhi=dPhi+2*TMath::Pi();
1779                 Float_t dEta=etaOthers-etaEventAxis;
1780
1781
1782                 fDPhiDEtaEventAxis[step]->Fill(dPhi, dEta, strangeWeightEventAxis*strangeWeightOthers);
1783                 fDPhiEventAxis[step]->Fill(dPhi, strangeWeightEventAxis*strangeWeightOthers);
1784
1785                 //check outliers
1786                 if(ptEventAxis< 0.4 || ptEventAxis > 100) AliInfo("particles out of range pt");
1787                 if(ntracksCharged<-1 || ntracksCharged>1500) AliInfo("particles out of range ncharge");
1788                 if(TMath::Abs(dEta)>2*fEtaCut) {
1789                     AliInfo("particles out of range dEta");
1790                     AliInfo(Form("eta1=%f, eta2=%f", etaOthers, etaEventAxis));
1791                     AliInfo(Form("step=%d",step));
1792                 }
1793                 if(dPhi<-0.5*TMath::Pi() || dPhi>1.5*TMath::Pi()){
1794                     AliInfo(Form("particles out of range dPhi"));
1795                     AliInfo(Form("phi1=%f, phi2=%f", phiOthers, phiEventAxis));
1796                 }
1797
1798                 Bool_t isLikeSign = CheckLikeSign(chargeEventAxis, chargeOthers);
1799
1800                 Double_t prop6[6] = {ptEventAxis,ptOthers,dEta,dPhi,ntracksCharged, static_cast<Double_t>(isLikeSign) };
1801                 fMapPair[step]->Fill(prop6, strangeWeightEventAxis*strangeWeightOthers);
1802
1803                 //thrid track loop (Andreas: three particle auto-correlations)
1804                 if(fThreeParticleCorr){
1805                     for (Int_t third = iTrack+1; third < nAll; third++) {
1806                         if(pt[pindexInnerEta[iTrack]]<fTriggerPtCut) continue;
1807                         if(pt[pindexInnerEta[third]]<fTriggerPtCut) continue;
1808
1809                         //dphi1
1810                         Float_t dPhi1 = phiEventAxis - phiOthers;
1811                         if(dPhi1>1.5*TMath::Pi()) dPhi1 = dPhi1-2*TMath::Pi();
1812                         else if(dPhi1<-0.5*TMath::Pi())dPhi1=dPhi1+2*TMath::Pi();
1813
1814                         Float_t phiThird = phi[pindexInnerEta[third]];
1815                         Float_t strangeWeightThird = strangeWeight[pindexInnerEta[third]];
1816
1817                         //dphi2
1818                         Float_t dPhi2 = phiEventAxis - phiThird;
1819                         if(dPhi2>1.5*TMath::Pi()) dPhi2 = dPhi2-2*TMath::Pi();
1820                         else if(dPhi2<-0.5*TMath::Pi())dPhi2=dPhi2+2*TMath::Pi();
1821
1822                         fDPhi1DPhi2[step]-> Fill(dPhi1, dPhi2, strangeWeightEventAxis*strangeWeightOthers*strangeWeightThird);
1823                         Double_t propThree[3] = {dPhi1,dPhi2,ntracksCharged};
1824                         fMapThree[step]->Fill(propThree,strangeWeightEventAxis*strangeWeightOthers*strangeWeightThird );
1825
1826
1827                     }// end of three particle correlation loop
1828
1829                 }// if fThreeParticleCorr is set to true
1830
1831             }// end of inner track loop
1832
1833         } //end of outer track loop
1834
1835         fTriggerNch[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1836         fTriggerNchSeeds[step]->Fill(ntracksCharged,highPtTracksInnerEta);
1837         fTriggerTracklet[step]->Fill(ntracklets);
1838
1839
1840     }//if there is at least one high pt track
1841
1842
1843     if(pindexInnerEta){// clean up array memory used for TMath::Sort
1844         delete[] pindexInnerEta;
1845         pindexInnerEta=0;
1846     }
1847
1848     if(ptInnerEta){// clean up array memory used for TMath::Sort
1849         delete[] ptInnerEta;
1850         ptInnerEta=0;
1851     }
1852
1853 }
1854
1855
1856
1857 //________________________________________________________________________
1858 void AliAnalysisTaskMinijet::Terminate(Option_t*)
1859 {
1860     //terminate function is called at the end
1861     //can be used to draw histograms etc.
1862
1863
1864 }
1865
1866 //________________________________________________________________________
1867 Bool_t AliAnalysisTaskMinijet::SelectParticlePlusCharged(const Short_t charge, const Int_t pdg, Bool_t prim)
1868 {
1869     //selection of mc particle
1870     //fSelectParticles=0: use charged primaries and pi0 and k0
1871     //fSelectParticles=1: use only charged primaries
1872     //fSelectParticles=2: use only pi0 and k0
1873
1874     if(fSelectParticles==0 || fSelectParticles==2){ // in case of 2: need to count also charged particles
1875         if(charge==0){
1876             if(!(pdg==111||pdg==130||pdg==310))
1877                 return false;
1878         }
1879         else{// charge !=0
1880             if(!prim)
1881                 return false;
1882         }
1883     }
1884
1885     else if(fSelectParticles==1){
1886         if (charge==0 || !prim){
1887             return false;
1888         }
1889     }
1890
1891     else{
1892         AliInfo("Error: wrong selection of charged/pi0/k0");
1893         return 0;
1894     }
1895
1896     return true;
1897
1898 }
1899
1900 //________________________________________________________________________
1901 Bool_t AliAnalysisTaskMinijet::SelectParticle(const Short_t charge, const Int_t pdg, const Bool_t prim)
1902 {
1903     //selection of mc particle
1904     //fSelectParticles=0: use charged primaries and pi0 and k0
1905     //fSelectParticles=1: use only charged primaries
1906     //fSelectParticles=2: use only pi0 and k0
1907
1908     if(fSelectParticles==0){
1909         if(charge==0){
1910             if(!(pdg==111||pdg==130||pdg==310))
1911                 return false;
1912         }
1913         else{// charge !=0
1914             if(!prim)
1915                 return false;
1916         }
1917     }
1918     else if (fSelectParticles==1){
1919         if (charge==0 || !prim){
1920             return false;
1921         }
1922     }
1923     else if(fSelectParticles==2){
1924         if(!(pdg==111||pdg==130||pdg==310))
1925             return false;
1926     }
1927
1928     return true;
1929
1930 }
1931
1932 //________________________________________________________________________
1933 Bool_t AliAnalysisTaskMinijet::CheckEvent(const Bool_t recVertex)
1934 {
1935     // This function tests the quality of an event (ESD/AOD) (rec and/or mc part)
1936     // recVertex=false:  check if Mc events and stack is there, Nmc>0
1937     // recVertex=false: " + check if there is a good, reconstructed SPD vertex
1938     // defined by |z|<fVertexCut(10cm), Contributer>0, no PileUpFromSPD(3,0,8)
1939
1940     if(fMode==0){//esd
1941
1942         //mc
1943         if(fUseMC){
1944
1945             //mc event
1946             AliMCEvent *mcEvente = (AliMCEvent*) MCEvent();
1947             if (!mcEvente) {
1948                 Error("CheckEvent", "Could not retrieve MC event");
1949                 return false;
1950             }
1951
1952             //stack
1953             AliStack* stackg = MCEvent()->Stack();
1954             if(!stackg) return false;
1955             Int_t ntracksg = mcEvente->GetNumberOfTracks();
1956             if(ntracksg<0) return false;
1957
1958             //vertex
1959             AliGenEventHeader*  headerg = MCEvent()->GenEventHeader();
1960             TArrayF mcVg;
1961             headerg->PrimaryVertex(mcVg);
1962             //if(TMath::Abs(mcVg[0])<1e-8 && TMath::Abs(mcVg[0])<1e-8 &&
1963             //   TMath::Abs(mcVg[0])<1e-8) return false;
1964             Float_t vzMCg = mcVg[2];
1965             if(TMath::Abs(vzMCg)>fVertexZCut) return false;
1966             //hasVtxMc=true;
1967         }
1968
1969         //rec
1970         if(recVertex==true){
1971             if(fESDEvent->IsPileupFromSPD(3,0,8))return false;
1972
1973             //rec vertex
1974             const AliESDVertex* vertexESDg   = fESDEvent->GetPrimaryVertex(); // uses track or SPD vertexer
1975             if(!vertexESDg) return false;
1976             fVertexCheck->Fill(vertexESDg->GetNContributors());
1977             if(vertexESDg->GetNContributors()<=0)return false;
1978             Float_t fVzg= vertexESDg->GetZ();
1979             if(TMath::Abs(fVzg)>fVertexZCut) return false;
1980
1981             //rec spd vertex
1982             const AliESDVertex *vtxSPD = fESDEvent->GetPrimaryVertexSPD();
1983             if(!vtxSPD) return false;
1984             if(vtxSPD->GetNContributors()<=0)return false;
1985             Float_t fVzSPD= vtxSPD->GetZ();
1986             if(TMath::Abs(fVzSPD)>fVertexZCut) return false;
1987
1988         }
1989         return true;
1990     }
1991
1992
1993     else if(fMode==1){ //aod
1994
1995         if(fUseMC){
1996
1997             //retreive MC particles from event
1998             TClonesArray *mcArray = (TClonesArray*)fAODEvent->
1999             FindListObject(AliAODMCParticle::StdBranchName());
2000             if(!mcArray){
2001                 AliInfo("No MC particle branch found");
2002                 return false;
2003             }
2004
2005             //mc
2006             AliAODMCHeader *aodMCheader = (AliAODMCHeader *) fAODEvent->
2007             FindListObject(AliAODMCHeader::StdBranchName());
2008             Float_t vzMC = aodMCheader->GetVtxZ();
2009             if(TMath::Abs(vzMC)>fVertexZCut) return false;
2010
2011             //hasVtxMc=true;
2012         }
2013
2014         //rec
2015         if(recVertex==true){
2016             //if(fAODEvent->IsPileupFromSPD(3,0.8)) return false;
2017
2018             AliAODVertex*       vertex= (AliAODVertex*)fAODEvent->GetPrimaryVertex();
2019             if(!vertex) return false;
2020             TString vtxTitle(vertex->GetTitle());// only allow vertex from tracks, no vertexer z
2021             // Printf("vtxTitle: %s",vtxTitle.Data());
2022             //if (!(vtxTitle.Contains("VertexerTracksWithConstraint"))) return false;
2023             fVertexCheck->Fill(vertex->GetNContributors());
2024             if(vertex->GetNContributors()<=0) return false;
2025             Double_t vzAOD=vertex->GetZ();
2026             // if(TMath::Abs(vzAOD)<1e-9) return false;
2027             if(TMath::Abs(vzAOD)>fVertexZCut) return false;
2028
2029             AliAODVertex*       vertexSPD= (AliAODVertex*)fAODEvent->GetPrimaryVertexSPD();
2030             if(!vertexSPD) return false;
2031             if(vertexSPD->GetNContributors()<=0) return false;
2032             Double_t vzSPD=vertexSPD->GetZ();
2033             //if(TMath::Abs(vzSPD)<1e-9) return false;
2034             if(TMath::Abs(vzSPD)>fVertexZCut) return false;
2035
2036             //check TPC reconstruction: check for corrupted chunks
2037             //better: check TPCvertex, but this is not available yet in AODs
2038             Int_t nAcceptedTracksTPC=0;
2039             Int_t nAcceptedTracksITSTPC=0;
2040             for (Int_t iTracks = 0; iTracks < fAODEvent->GetNumberOfTracks(); iTracks++) {
2041                 AliAODTrack *track = dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(iTracks));
2042                 if(!track) AliFatal("Not a standard AOD");
2043                 if (!track) continue;
2044                 if(track->TestFilterBit(128) && TMath::Abs(track->Eta())<fEtaCut &&
2045                    track->Pt()>fPtMin && track->Pt()<fPtMax)
2046                     nAcceptedTracksTPC++;
2047                 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut &&
2048                    track->Pt()>fPtMin && track->Pt()<fPtMax)
2049                     nAcceptedTracksITSTPC++;
2050             }
2051             fCorruptedChunks->Fill(nAcceptedTracksTPC,nAcceptedTracksITSTPC);
2052             if(fRejectChunks){
2053                 if(nAcceptedTracksTPC>fNTPC && nAcceptedTracksITSTPC==0)
2054                     return false;//most likely corrupted chunk. No ITS tracks are reconstructed
2055             }
2056             fCorruptedChunksAfter->Fill(nAcceptedTracksTPC,nAcceptedTracksITSTPC);
2057
2058             //control histograms=================
2059             //tracklet loop
2060             Int_t ntrackletsAccept=0;
2061             AliAODTracklets * mult= (AliAODTracklets*)fAODEvent->GetTracklets();
2062             for(Int_t i=0;i<mult->GetNumberOfTracklets();++i){
2063                 if(TMath::Abs(mult->GetDeltaPhi(i))<0.05 &&
2064                    TMath::Abs(TMath::Log(TMath::Tan(0.5 * mult->GetTheta(i))))<fEtaCut) ++ntrackletsAccept;
2065             }
2066             Int_t nAcceptedTracks=0;
2067             for (Int_t iTracks = 0; iTracks < fAODEvent->GetNumberOfTracks(); iTracks++) {
2068                 AliAODTrack *track = dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(iTracks));
2069                 if(!track) AliFatal("Not a standard AOD");
2070                 if (!track) continue;
2071                 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta())<fEtaCut
2072                    && track->Pt()>fPtMin && track->Pt()<fPtMax) nAcceptedTracks++;
2073             }
2074             fNContrNtracklets->Fill(ntrackletsAccept,vertexSPD->GetNContributors());
2075             fNContrNtracks->Fill(nAcceptedTracks,vertexSPD->GetNContributors());
2076             //====================================
2077         }
2078         return true;
2079
2080     }
2081
2082     else {
2083         Printf("Wrong mode!");
2084         return false;
2085     }
2086
2087 }
2088
2089 //_____________________________________________________________________________
2090 const Double_t * AliAnalysisTaskMinijet::CreateLogAxis(const Int_t nbins,
2091                                                        const Double_t xmin,
2092                                                        const Double_t xmax)
2093 {
2094     // returns pointer to an array with can be used to build a logarithmic axis
2095     // it is user responsibility to delete the array
2096
2097     Double_t logxmin = TMath::Log10(xmin);
2098     Double_t logxmax = TMath::Log10(xmax);
2099     Double_t binwidth = (logxmax-logxmin)/nbins;
2100
2101     Double_t *xbins =  new Double_t[nbins+1];
2102
2103     xbins[0] = xmin;
2104     for (Int_t i=1;i<=nbins;i++) {
2105         xbins[i] = xmin + TMath::Power(10,logxmin+i*binwidth);
2106     }
2107
2108     return xbins;
2109 }
2110
2111 //_____________________________________________________________________________
2112 Bool_t AliAnalysisTaskMinijet::CheckLikeSign(const Short_t chargeEventAxis,
2113                                              const Short_t chargeOthers)
2114 {
2115     // compute if charge of two particles/tracks has same sign or different sign
2116
2117     if(fDebug>2) Printf("Charge1=%d, Charge2=%d",chargeEventAxis,chargeOthers);
2118
2119     if(chargeEventAxis<0){
2120         if(chargeOthers<0){
2121             return true;
2122         }
2123         else if(chargeOthers>0){
2124             return false;
2125         }
2126     }
2127
2128     else if(chargeEventAxis>0){
2129         if(chargeOthers>0){
2130             return true;
2131         }
2132         else if(chargeOthers<0){
2133             return false;
2134         }
2135     }
2136
2137     else{
2138         AliInfo("Error: Charge not lower nor higher as zero");
2139         return false;
2140     }
2141
2142     AliInfo("Error: Check values of Charge");
2143     return false;
2144 }
2145
2146