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