PWGCF/EBYE/Fluctuations/AliEbyEHigherMomentsTask.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Satyajit Jena.                                                 *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 //=========================================================================//
  18 //                                                                         //
  19 //           Analysis Task for Net-Charge Higher Moment Analysis           //
  20 //              Author: Satyajit Jena || Nirbhay K. Behera                 //
  21 //                      sjena@cern.ch || nbehera@cern.ch                   //
  22 //                                                                         //
  23 //=========================================================================//
  24
  25 #include "TChain.h"
  26 #include "TList.h"
  27 #include "TFile.h"
  28 #include "TTree.h"
  29 #include "TH1D.h"
  30 #include "TH2D.h"
  31 #include "TH3D.h"
  32 #include "THnSparse.h"
  33 #include "TCanvas.h"
  34
  35 #include "AliAnalysisTaskSE.h"
  36 #include "AliAnalysisManager.h"
  37
  38 #include "AliESDVertex.h"
  39 #include "AliESDEvent.h"
  40 #include "AliESDInputHandler.h"
  41 #include "AliAODEvent.h"
  42 #include "AliAODTrack.h"
  43 #include "AliAODInputHandler.h"
  44 #include "AliESD.h"
  45 #include "AliESDEvent.h"
  46 #include "AliAODEvent.h"
  47 #include "AliStack.h"
  48 #include "AliESDtrackCuts.h"
  49 #include "AliAODMCHeader.h"
  50 #include "AliAODMCParticle.h"
  51 #include "AliMCEventHandler.h"
  52 #include "AliMCEvent.h"
  53 #include "AliStack.h"
  54 #include "AliGenHijingEventHeader.h"
  55 #include "AliGenEventHeader.h"
  56 #include "AliPID.h"
  57 #include "AliAODPid.h"
  58 #include "AliPIDResponse.h"
  59 #include "AliAODpidUtil.h"
  60 #include "AliPIDCombined.h"
  61
  62 #include "AliEbyEHigherMomentsTask.h"
  63
  64 using std::cout;
  65 using std::endl;
  66
  67 ClassImp(AliEbyEHigherMomentsTask)
  68
  69
  70 //-----------------------------------------------------------------------
  71 AliEbyEHigherMomentsTask::AliEbyEHigherMomentsTask( const char *name )
  72 : AliAnalysisTaskSE( name ),
  73   fListOfHistos(0),
  74   fAOD(0),
  75   fArrayMC(0),
  76   fPIDResponse(0),
  77   fParticleSpecies(AliPID::kProton),
  78   fAnalysisType(0),
  79   fCentralityEstimator("V0M"),
  80   fCentrality(0),
  81   fVxMax(3.),
  82   fVyMax(3.),
  83   fVzMax(10.),
  84   fDCAxy(2.4),
  85   fDCAz(3.),
  86   fPtLowerLimit(0.3),
  87   fPtHigherLimit(1.5),
  88   fEtaLowerLimit(-0.8),
  89   fEtaHigherLimit(0.8),
  90   fRapidityCut(0.5),
  91   fNSigmaCut(3.),
  92   fTPCNClus(80),
  93   fChi2perNDF(4.),
  94   fAODtrackCutBit(128),
  95   fUsePid(kFALSE),
  96   fEventCounter(0),
  97   fHistDCA(0),
  98   fTPCSig(0),
  99   fTPCSigA(0),
 100   fTHnCentNplusNminusCh(0),
 101   fTHnCentNplusNminusChTruth(0),
 102   fTHnCentNplusNminus(0)
 103 {
 104
 105   for ( Int_t i = 0; i < 13; i++) {
 106     fHistQA[i] = NULL;
 107   }
 108
 109   for ( Int_t i = 0; i < 5; i++ ){
 110     fTHnCentNplusNminusPid[i] = NULL;
 111     fTHnCentNplusNminusPidTruth[i] = NULL;
 112   }
 113
 114   DefineOutput(1, TList::Class()); // Outpput....
 115   //DefineOutput(2, TList::Class());
 116 }
 117
 118 AliEbyEHigherMomentsTask::~AliEbyEHigherMomentsTask() {
 119   //if(fListOfHistosQA) delete fListOfHistosQA;
 120   if(fListOfHistos) delete fListOfHistos;
 121
 122 }
 123
 124 //---------------------------------------------------------------------------------
 125 void AliEbyEHigherMomentsTask::UserCreateOutputObjects() {
 126
 127   // fListOfHistosQA = new TList();
 128   //fListOfHistosQA->SetOwner(kTRUE);
 129   fListOfHistos = new TList();
 130   fListOfHistos->SetOwner(kTRUE);
 131
 132   fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
 133   fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
 134   fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-90% centrality");
 135   fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
 136   fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
 137   fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
 138   fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
 139   fListOfHistos->Add(fEventCounter);
 140
 141   //For QA-Histograms
 142   fHistQA[0] = new TH1D("fHistQAvx", "Histo Vx After Cut", 400, -4., 4.);
 143   fHistQA[1] = new TH1D("fHistQAvy", "Histo Vy After Cut", 400, -4., 4.);
 144   fHistQA[2] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
 145   fHistQA[3] = new TH1D("fHistQAvxA", "Histo Vx All", 500, -5., 5.);
 146   fHistQA[4] = new TH1D("fHistQAvyA", "Histo Vy All", 500, -5., 5.);
 147   fHistQA[5] = new TH1D("fHistQAvzA", "Histo Vz All", 500, -25., 25.);
 148   fHistQA[6] = new TH1D("fHistQADcaXyC", "Histo DCAxy after cut", 500, -5., 5.);
 149   fHistQA[7] = new TH1D("fHistQADcaZC", "Histo DCAz after cut", 500, -5., 5.);
 150   fHistQA[8] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
 151   fHistQA[9] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
 152   fHistQA[10] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
 153   fHistQA[11] = new TH1D("fHistQANCls","Number of TPC cluster",200,0,200);
 154   fHistQA[12] = new TH1D("fHistQAChi2","Chi2 per NDF",100,0,10);
 155   for(Int_t i = 0; i < 13; i++)
 156     {
 157       fListOfHistos->Add(fHistQA[i]);
 158     }
 159
 160   fHistDCA = new TH2D("fHistDCA","DCAxy Vs DCAz", 500, -5., 5., 500, -5., 5.);
 161   fTPCSig = new TH2D("fTPCSig","TPC signal",200, 0.0, 10. ,1000,0.,1000);
 162   fTPCSig->SetMarkerColor(kRed);
 163   fTPCSigA = new TH2D("fTPCSigA","TPC signal all ",200, 0.0, 10. ,1000,0.,1000);
 164   fListOfHistos->Add(fHistDCA);
 165   fListOfHistos->Add(fTPCSig);
 166   fListOfHistos->Add(fTPCSigA);
 167
 168   const Int_t nDim = 3;
 169   const Int_t nPid = 5;
 170   TString Species[nPid] = {"Electron","Muon","Pion","Kaon","Proton"};
 171
 172   Int_t fBins[nPid][nDim];
 173   Double_t fMin[nPid][nDim];
 174   Double_t fMax[nPid][nDim];
 175
 176   for( Int_t i = 0; i < nPid; i++ ){
 177     for( Int_t j = 0; j < nDim; j++ ){
 178       fBins[i][j] = 0;
 179       fMin[i][j] = 0.;
 180       fMax[i][j] = 0.;
 181     }
 182   }
 183
 184
 185   Int_t fBinsCh[nDim] = {100, 1500, 1500};
 186   Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 };
 187   Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5};
 188
 189   fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
 190   fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality");
 191   fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus");
 192   fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus");
 193   fListOfHistos->Add(fTHnCentNplusNminusCh);
 194
 195   if( fAnalysisType == "MCAOD" ){
 196
 197     fTHnCentNplusNminusChTruth = new THnSparseD("fTHnCentNplusNminusChTruth","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
 198     fTHnCentNplusNminusChTruth->GetAxis(0)->SetTitle("Centrality");
 199     fTHnCentNplusNminusChTruth->GetAxis(1)->SetTitle("Nplus");
 200     fTHnCentNplusNminusChTruth->GetAxis(2)->SetTitle("Nminus");
 201     fListOfHistos->Add(fTHnCentNplusNminusChTruth);
 202
 203
 204   }//MCAOD---
 205
 206   TString hname1, hname11;
 207   TString htitle1, axisTitle1, axisTitle2;
 208
 209
 210   if( fUsePid ){
 211
 212     Int_t gPid = (Int_t)fParticleSpecies;
 213
 214     if( gPid > 1 && gPid < 5 ){
 215       //Pion----
 216       fBins[2][0] = 100; fBins[2][1] = 1000; fBins[2][2] = 1000;
 217       fMin[2][0] = -0.5; fMin[2][1] = -0.5; fMin[2][2] = -0.5;
 218       fMax[2][0] = 99.5; fMax[2][1] = 999.5; fMax[2][2] = 999.5;
 219       //Kaon------
 220       fBins[3][0] = 100; fBins[3][1] = 600; fBins[3][2] = 600;
 221       fMin[3][0] = -0.5; fMin[3][1] = -0.5; fMin[3][2] = -0.5;
 222       fMax[3][0] = 99.5; fMax[3][1] = 599.5; fMax[3][2] = 599.5;
 223       //Proton-----
 224       fBins[4][0] = 100; fBins[4][1] = 400; fBins[4][2] = 400;
 225       fMin[4][0] = -0.5; fMin[4][1] = -0.5; fMin[4][2] = -0.5;
 226       fMax[4][0] = 99.5; fMax[4][1] = 399.5; fMax[4][2] = 399.5;
 227
 228       hname1 = "fCentNplusNminusPid"; hname1 +=gPid;
 229       htitle1 = Species[gPid]; htitle1 +=" And Neg-"; htitle1 +=Species[gPid];
 230       axisTitle1 ="Pos-"; axisTitle1 += Species[gPid];
 231       axisTitle2 ="Neg-"; axisTitle2 += Species[gPid];
 232
 233       fTHnCentNplusNminusPid[gPid] = new THnSparseD(hname1.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
 234
 235       fTHnCentNplusNminusPid[gPid]->GetAxis(0)->SetTitle("Centrality");
 236       fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
 237       fTHnCentNplusNminusPid[gPid]->GetAxis(2)->SetTitle(axisTitle2.Data());
 238
 239       fListOfHistos->Add(fTHnCentNplusNminusPid[gPid]);
 240
 241       if( fAnalysisType == "MCAOD" ){
 242
 243         hname11 = "fCentNplusNminusPidTruth"; hname11 +=gPid;
 244         fTHnCentNplusNminusPidTruth[gPid] = new THnSparseD(hname11.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
 245
 246         fTHnCentNplusNminusPidTruth[gPid]->GetAxis(0)->SetTitle("Centrality");
 247         fTHnCentNplusNminusPidTruth[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
 248         fTHnCentNplusNminusPidTruth[gPid]->GetAxis(2)->SetTitle(axisTitle2.Data());
 249         fListOfHistos->Add(fTHnCentNplusNminusPidTruth[gPid]);
 250
 251       }//MCAOD-----
 252     }//Pion, Koan and Proton-------
 253     else{
 254
 255       Int_t fBinsX[nDim] = {100, 1500, 1500};
 256       Double_t fMinX[nDim] = { -0.5, -0.5, -0.5 };
 257       Double_t fMaxX[nDim] = { 99.5, 1499.5, 1499.5};
 258
 259       fTHnCentNplusNminus = new THnSparseD("fTHnCentNplusNminus","Cent-NplusChrg-NminusChrg", nDim, fBinsX, fMinX, fMaxX);
 260       fTHnCentNplusNminus->GetAxis(0)->SetTitle("Centrality");
 261       fTHnCentNplusNminus->GetAxis(1)->SetTitle("UnwantedPlus");
 262       fTHnCentNplusNminus->GetAxis(2)->SetTitle("UnwantedMinus");
 263       fListOfHistos->Add(fTHnCentNplusNminus);
 264
 265     }//Unwanted particle -(other than pi, K or proton)
 266
 267   }//fUsePid-------
 268
 269   //PostData(1, fListOfHistosQA);
 270   PostData(1, fListOfHistos);
 271
 272
 273 }
 274
 275
 276 //----------------------------------------------------------------------------------
 277 void AliEbyEHigherMomentsTask::UserExec( Option_t * ){
 278
 279   fEventCounter->Fill(1);
 280
 281   if(fAnalysisType == "AOD") {
 282
 283     doAODEvent();
 284
 285   }//AOD--analysis-----
 286
 287   else if(fAnalysisType == "MCAOD") {
 288
 289     doMCAODEvent();
 290
 291   }
 292
 293   else return;
 294
 295   fEventCounter->Fill(8);
 296
 297
 298 }
 299
 300 //--------------------------------------------------------------------------------------
 301 void AliEbyEHigherMomentsTask::doAODEvent(){
 302
 303   Double_t positiveSum = 0.;
 304   Double_t negativeSum = 0.;
 305   Double_t posPidSum = 0.;
 306   Double_t negPidSum = 0.;
 307   Int_t gPid = 0;
 308
 309   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
 310   if (!manager) {
 311     cout<<"ERROR: Analysis manager not found."<<endl;
 312     return;
 313   }
 314   //coneect to the inputHandler------------
 315   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
 316   if (!inputHandler) {
 317     cout<<"ERROR: Input handler not found."<<endl;
 318     return;
 319   }
 320
 321   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 322   if (!fAOD)
 323     {
 324       cout<< "ERROR 01: AOD not found " <<endl;
 325       return;
 326     }
 327
 328   fPIDResponse =inputHandler->GetPIDResponse();
 329
 330
 331   if (!fPIDResponse){
 332     AliFatal("This Task needs the PID response attached to the inputHandler");
 333     return;
 334   }
 335
 336   AliAODHeader *aodHeader = fAOD->GetHeader();
 337
 338   fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 339
 340   if(fCentrality < 0 || fCentrality >= 91) return;
 341
 342
 343
 344   fEventCounter->Fill(2);
 345
 346   if(!ProperVertex()) return;
 347
 348   Int_t nTracks = fAOD->GetNumberOfTracks();
 349
 350   TExMap *trackMap = new TExMap();//Mapping matrix----
 351
 352   for(Int_t i = 0; i < nTracks; i++) {
 353
 354     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
 355
 356     if(!aodTrack) {
 357       AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
 358       continue;
 359     }
 360
 361     Double_t tpcSignalAll = aodTrack->GetTPCsignal();
 362     fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
 363
 364     Int_t gID = aodTrack->GetID();
 365
 366     if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
 367   }//1st track loop----
 368
 369   AliAODTrack* newAodTrack;
 370
 371   for( Int_t j = 0; j < nTracks; j++ )
 372     {
 373
 374       AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
 375
 376       if(!aodTrack1) {
 377         AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
 378         continue;
 379       }
 380
 381
 382       if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
 383
 384       Int_t gID = aodTrack1->GetID();
 385
 386       //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
 387       newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
 388
 389
 390       Double_t pt = aodTrack1->Pt();
 391       Double_t eta = aodTrack1->Eta();
 392       Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
 393       Double_t chi2ndf = aodTrack1->Chi2perNDF();
 394
 395       if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
 396       if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
 397       if( nclus < fTPCNClus ) continue;
 398       if( chi2ndf > fChi2perNDF ) continue;
 399
 400       if( fAODtrackCutBit == 128 ){
 401         //TPC only tracks----
 402         Float_t dxy = 0., dz = 0.;
 403         dxy = aodTrack1->DCA();
 404         dz  = aodTrack1->ZAtDCA();
 405
 406         if( fabs(dxy) > fDCAxy ) continue;
 407         if( fabs(dz) > fDCAz ) continue;
 408         //Extra cut on DCA---( Similar to BF Task.. )
 409         if( fDCAxy !=0 && fDCAz !=0 ){
 410           if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
 411         }
 412
 413         fHistQA[6]->Fill(dxy);
 414         fHistQA[7]->Fill(dz);
 415         fHistDCA->Fill(dxy,dz);
 416
 417       }
 418
 419       fHistQA[8]->Fill(pt);
 420       fHistQA[9]->Fill(eta);
 421       fHistQA[10]->Fill(aodTrack1->Phi());
 422       fHistQA[11]->Fill(nclus);
 423       fHistQA[12]->Fill(chi2ndf);
 424
 425       Short_t gCharge = aodTrack1->Charge();
 426
 427       if(gCharge > 0) positiveSum += 1.;
 428       if(gCharge < 0) negativeSum += 1.;
 429
 430       if( fUsePid ) {
 431
 432         gPid = (Int_t)fParticleSpecies;
 433
 434         Double_t rap =  newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
 435         Double_t tpcSignal = newAodTrack->GetTPCsignal();
 436         //Double_t rap =  aodTrack1->Y(AliPID::ParticleMass(fParticleSpecies));
 437         //Double_t tpcSignal = aodTrack1->GetTPCsignal();
 438
 439         if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
 440
 441         fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
 442
 443         Float_t nsigmaTPCPID = -999.;
 444         Float_t nsigmaTOFPID = -999.;
 445         //Float_t nsigmaTPCTOFPID = -999.;
 446
 447         nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
 448         nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
 449
 450         if ( nsigmaTPCPID < fNSigmaCut  ){
 451
 452           if (gCharge > 0) posPidSum +=1.;
 453           if( gCharge < 0 ) negPidSum +=1.;
 454
 455         }
 456       }//fUsepid-----
 457
 458     }//--------- Track Loop to select with filterbit
 459
 460
 461   //cout << fCentrality <<" "<< nPlusCharge << " " << nMinusCharge << endl;
 462   //cout << fCentrality <<" "<< nParticle << " " << nAntiParticle << endl;
 463
 464   Double_t fContainerCh[3] = { fCentrality, positiveSum, negativeSum};
 465
 466
 467
 468   fTHnCentNplusNminusCh->Fill(fContainerCh);
 469
 470   if( fUsePid ){
 471     gPid = (Int_t)fParticleSpecies;
 472     Double_t fContainerPid[3] = { fCentrality, posPidSum, negPidSum};
 473     fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
 474   }
 475
 476   fEventCounter->Fill(7);
 477   PostData(1, fListOfHistos);
 478   return;
 479
 480 }
 481 //--------------------------------------------------------------------------------------
 482 //--------------------------------------------------------------------------------------
 483 void AliEbyEHigherMomentsTask::doMCAODEvent(){
 484
 485
 486   //---------
 487   Double_t positiveSumMCRec = 0.;
 488   Double_t negativeSumMCRec = 0.;
 489   Double_t posPidSumMCRec = 0.;
 490   Double_t negPidSumMCRec = 0.;
 491
 492   Double_t positiveSumMCTruth = 0.;
 493   Double_t negativeSumMCTruth = 0.;
 494   Double_t posPidSumMCTruth = 0.;
 495   Double_t negPidSumMCTruth = 0.;
 496
 497   Int_t gPid = 0;
 498   Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies);
 499
 500   //Connect to Anlaysis manager------
 501   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
 502   if (!manager) {
 503     cout<<"ERROR: Analysis manager not found."<<endl;
 504     return;
 505   }
 506
 507   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
 508   if (!inputHandler) {
 509     cout<<"ERROR: Input handler not found."<<endl;
 510     return;
 511   }
 512
 513   //AOD event------
 514   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 515   if (!fAOD)
 516     {
 517       cout<< "ERROR 01: AOD not found " <<endl;
 518       return;
 519     }
 520
 521   fPIDResponse =inputHandler->GetPIDResponse();
 522
 523
 524   if (!fPIDResponse){
 525     AliFatal("This Task needs the PID response attached to the inputHandler");
 526     return;
 527   }
 528
 529   // -- Get MC header
 530   // ------------------------------------------------------------------
 531
 532   fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 533   if (!fArrayMC) {
 534     AliFatal("Error: MC particles branch not found!\n");
 535     return;
 536   }
 537
 538   AliAODMCHeader *mcHdr=NULL;
 539   mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 540   if(!mcHdr) {
 541     printf("MC header branch not found!\n");
 542     return;
 543   }
 544
 545   AliAODHeader *aodHeader = fAOD->GetHeader();
 546
 547   fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 548
 549
 550
 551   if( fCentrality < 0 || fCentrality >= 91) return;
 552
 553   fEventCounter->Fill(2);
 554
 555   if(!ProperVertex()) return;
 556
 557   Int_t nTracks = fAOD->GetNumberOfTracks();
 558
 559   TExMap *trackMap = new TExMap();//Mapping matrix----
 560
 561   for(Int_t i = 0; i < nTracks; i++) {
 562
 563     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
 564
 565     if(!aodTrack) {
 566       AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
 567       continue;
 568     }
 569
 570     Double_t tpcSignalAll = aodTrack->GetTPCsignal();
 571     fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
 572
 573     Int_t gID = aodTrack->GetID();
 574
 575     if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global(primary) tracks-----
 576
 577   }//1st track loop----
 578
 579   AliAODTrack* newAodTrack;
 580
 581   for( Int_t j = 0; j < nTracks; j++ ){
 582
 583     AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
 584
 585     if(!aodTrack1) {
 586       AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
 587       continue;
 588     }
 589
 590
 591     if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
 592
 593     Int_t gID = aodTrack1->GetID();
 594
 595     //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
 596
 597     newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
 598
 599
 600     //cout << dxy << endl;
 601     Double_t pt = aodTrack1->Pt();
 602     Double_t eta = aodTrack1->Eta();
 603     Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
 604     Double_t chi2ndf = aodTrack1->Chi2perNDF();
 605
 606     if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
 607     if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
 608     if( nclus < fTPCNClus ) continue;
 609     if( chi2ndf > fChi2perNDF ) continue;
 610
 611     if( fAODtrackCutBit == 128 ){
 612       Float_t dxy = 0., dz = 0.;
 613       dxy = aodTrack1->DCA();
 614       dz  = aodTrack1->ZAtDCA();
 615
 616       if( fabs(dxy) > fDCAxy ) continue;
 617       if( fabs(dz) > fDCAz ) continue;
 618       //Extra cut on DCA---( Similar to BF Task.. )
 619       if( fDCAxy !=0 && fDCAz !=0 ){
 620         if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
 621       }
 622
 623       fHistQA[6]->Fill(dxy);
 624       fHistQA[7]->Fill(dz);
 625       fHistDCA->Fill(dxy,dz);
 626
 627     }
 628
 629
 630
 631     fHistQA[8]->Fill(pt);
 632     fHistQA[9]->Fill(eta);
 633     fHistQA[10]->Fill(aodTrack1->Phi());
 634     fHistQA[11]->Fill(nclus);
 635     fHistQA[12]->Fill(chi2ndf);
 636
 637
 638     Short_t gCharge = aodTrack1->Charge();
 639
 640     if( gCharge == 0 ) continue;
 641
 642     if(gCharge > 0) positiveSumMCRec += 1.;
 643     if(gCharge < 0) negativeSumMCRec += 1.;
 644
 645     if( fUsePid ) {
 646
 647       gPid = (Int_t)fParticleSpecies;
 648       Double_t rap =  newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
 649       Double_t tpcSignal = newAodTrack->GetTPCsignal();
 650
 651       if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
 652
 653       fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
 654
 655       Float_t nsigmaTPCPID = -999.;
 656       Float_t nsigmaTOFPID = -999.;
 657       //Float_t nsigmaTPCTOFPID = -999.;
 658
 659       nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
 660       nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
 661
 662       if( nsigmaTPCPID < fNSigmaCut  ){
 663
 664         if (gCharge > 0) posPidSumMCRec +=1;
 665         if( gCharge < 0 ) negPidSumMCRec +=1.;
 666       }//nSigmaCut-----
 667     }//fUsepid-----
 668
 669   }//--------- Track Loop to select with filterbit
 670
 671
 672   Double_t fContainerCh[3] = { fCentrality, positiveSumMCRec, negativeSumMCRec};//Reco. values ch. hadrons
 673   fTHnCentNplusNminusCh->Fill(fContainerCh);//Fill the rec. ch. particles---
 674
 675   if( fUsePid ){
 676     gPid = (Int_t)fParticleSpecies;
 677     Double_t fContainerPid[3] = { fCentrality, posPidSumMCRec, negPidSumMCRec};//Reco. values pid.
 678     fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);//Fill the rec. pid tracks
 679   }
 680
 681
 682   //cout << fCentrality << " "<< nPlusCharge << " "<< nMinusCharge << endl;
 683   //cout << nParticle <<" And " << nAntiParticle << endl;
 684   //===========================================
 685   //--------MC Truth---------------------------
 686   //===========================================
 687
 688   Int_t nMCTrack = fArrayMC->GetEntriesFast();
 689
 690   for (Int_t iMC = 0; iMC < nMCTrack; iMC++) {
 691
 692     AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
 693
 694     if(!partMC){
 695       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
 696       continue;
 697     }
 698
 699     if(partMC->Charge() == 0) continue;
 700     if(!partMC->IsPhysicalPrimary()) continue;
 701
 702     if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
 703     if (partMC->Pt() < fPtLowerLimit ||  partMC->Pt() > fPtHigherLimit) continue;
 704
 705     Short_t gCharge = partMC->Charge();
 706
 707     if(gCharge > 0) positiveSumMCTruth += 1.;
 708     if(gCharge < 0) negativeSumMCTruth += 1.;
 709
 710     if(fUsePid){
 711
 712       Double_t rap =  partMC->Y();
 713       if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
 714       if(TMath::Abs(partMC->GetPdgCode()) != gPdgCode) continue;
 715
 716       if(gCharge > 0) posPidSumMCTruth += 1.;
 717       if(gCharge < 0) negPidSumMCTruth += 1.;
 718
 719     }//if(fUsePid) ----
 720
 721
 722   }//MC-Truth Track loop--
 723
 724   Double_t fContainerChTruth[3] = { fCentrality, positiveSumMCTruth, negativeSumMCTruth };
 725   fTHnCentNplusNminusChTruth->Fill(fContainerChTruth);//MC -Truth ch. particles
 726
 727   if( fUsePid ){
 728     gPid = (Int_t)fParticleSpecies;
 729     Double_t fContainerPidTruth[3] = { fCentrality, posPidSumMCTruth, negPidSumMCTruth};
 730     fTHnCentNplusNminusPidTruth[gPid]->Fill(fContainerPidTruth);//MC-Truth pid
 731   }
 732
 733   fEventCounter->Fill(7);
 734   PostData(1, fListOfHistos);
 735   return;
 736
 737 }
 738
 739 //---------------------------------------------------------------------------------------
 740 Bool_t AliEbyEHigherMomentsTask::ProperVertex(){
 741
 742   Bool_t IsVtx = kFALSE;
 743
 744   const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
 745
 746   if(vertex) {
 747     Double32_t fCov[6];
 748     vertex->GetCovarianceMatrix(fCov);
 749     if(vertex->GetNContributors() > 0) {
 750       if(fCov[5] != 0) {
 751
 752         Double_t lvx = vertex->GetX();
 753         Double_t lvy = vertex->GetY();
 754         Double_t lvz = vertex->GetZ();
 755
 756         fEventCounter->Fill(5);
 757
 758         fHistQA[3]->Fill(lvx);fHistQA[4]->Fill(lvy);fHistQA[5]->Fill(lvz);
 759
 760         if(TMath::Abs(lvx) < fVxMax) {
 761           if(TMath::Abs(lvy) < fVyMax) {
 762             if(TMath::Abs(lvz) < fVzMax) {
 763
 764               fEventCounter->Fill(6);
 765               fHistQA[0]->Fill(lvx);fHistQA[1]->Fill(lvy);fHistQA[2]->Fill(lvz);
 766
 767               IsVtx = kTRUE;
 768
 769             }//Z-Vertex cut---
 770           }//Y-vertex cut--
 771         }//X-vertex cut---
 772       }//Covariance------
 773     }//Contributors check---
 774   }//If vertex-----------
 775
 776   return IsVtx;
 777 }
 778 //---------------------------------------------------------------------------------------
 779
 780 //------------------------------------------------------------------------
 781 void AliEbyEHigherMomentsTask::Terminate( Option_t * ){
 782
 783   Info("AliEbyEHigherMomentTask"," Task Successfully finished");
 784
 785 }
 786
 787 //------------------------------------------------------------------------