PWGCF/EBYE/Fluctuations/AliEbyEHigherMomentsEffContTask.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 "TH2F.h"
  31 #include "THnSparse.h"
  32 #include "TCanvas.h"
  33 #include "TParticle.h"
  34 #include <TDatabasePDG.h>
  35 #include <TParticlePDG.h>
  36
  37 #include "AliAnalysisTaskSE.h"
  38 #include "AliAnalysisManager.h"
  39
  40 #include "AliESDVertex.h"
  41 #include "AliESDEvent.h"
  42 #include "AliESDInputHandler.h"
  43 #include "AliAODEvent.h"
  44 #include "AliAODTrack.h"
  45 #include "AliAODInputHandler.h"
  46 #include "AliESD.h"
  47 #include "AliESDEvent.h"
  48 #include "AliAODEvent.h"
  49 #include "AliStack.h"
  50 #include "AliESDtrackCuts.h"
  51 #include "AliAODMCHeader.h"
  52 #include "AliAODMCParticle.h"
  53 #include "AliMCEventHandler.h"
  54 #include "AliMCEvent.h"
  55 #include "AliMCParticle.h"
  56 #include "AliStack.h"
  57 #include "AliGenHijingEventHeader.h"
  58 #include "AliGenEventHeader.h"
  59 #include "AliPID.h"
  60 #include "AliAODPid.h"
  61 #include "AliPIDResponse.h"
  62 #include "AliAODpidUtil.h"
  63 #include "AliPIDCombined.h"
  64
  65 #include "AliEbyEHigherMomentsEffContTask.h"
  66
  67 using std::cout;
  68 using std::endl;
  69
  70
  71 ClassImp(AliEbyEHigherMomentsEffContTask)
  72
  73
  74 //-----------------------------------------------------------------------
  75 AliEbyEHigherMomentsEffContTask::AliEbyEHigherMomentsEffContTask( const char *name )
  76 : AliAnalysisTaskSE( name ),
  77   fListOfHistosQA(0),
  78   fListOfHistos(0),
  79   fAOD(0),
  80   fArrayMC(0),
  81   fPIDResponse(0),
  82   fParticleSpecies(AliPID::kProton),
  83   fAnalysisType("AOD"),
  84   fCentralityEstimator("V0M"),
  85   fCentrality(0),
  86   fVxMax(3.),
  87   fVyMax(3.),
  88   fVzMax(10.),
  89   fDCAxy(2.4),
  90   fDCAz(3.),
  91   fPtLowerLimit(0.3),
  92   fPtHigherLimit(1.5),
  93   fEtaLowerLimit(-0.8),
  94   fEtaHigherLimit(0.8),
  95   fRapidityCut(0.5),
  96   fNSigmaCut(3.),
  97   fTPCNClus(80),
  98   fChi2perNDF(4.),
  99   fAODtrackCutBit(1024),
 100   fLabel(NULL),
 101   fUsePid(kFALSE),
 102   fCheckCont(kFALSE),
 103   fEff(kFALSE),
 104   fEventCounter(0),
 105   fHistDCA(0),
 106   fTPCSig(0),
 107   fTPCSigA(0),
 108   fTHnCentNplusNminusCh(0),
 109   fTHnCentNplusNminus(0),
 110   fTHnEff(0),
 111   fTHnCont(0)
 112 {
 113
 114   for ( Int_t i = 0; i < 13; i++) {
 115     fHistQA[i] = NULL;
 116   }
 117
 118   for ( Int_t i = 0; i < 5; i++ ){
 119     fTHnCentNplusNminusPid[i] = NULL;
 120   }
 121
 122   DefineOutput(1, TList::Class()); // Outpput....
 123   DefineOutput(2, TList::Class());
 124
 125 }
 126
 127 AliEbyEHigherMomentsEffContTask::~AliEbyEHigherMomentsEffContTask() {
 128
 129   if(fListOfHistosQA) delete fListOfHistosQA;
 130   if(fListOfHistos) delete fListOfHistos;
 131   if(fLabel[0] ) delete [] (fLabel[0]);
 132   if(fLabel[1] ) delete [] (fLabel[1]);
 133
 134 }
 135
 136 //---------------------------------------------------------------------------------
 137 void AliEbyEHigherMomentsEffContTask::UserCreateOutputObjects() {
 138   //For QA-Histograms
 139   fListOfHistosQA = new TList();
 140   fListOfHistosQA->SetOwner(kTRUE);
 141
 142   fListOfHistos = new TList();
 143   fListOfHistos->SetOwner(kTRUE);
 144
 145   fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
 146   fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
 147   fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-90% centrality");
 148   fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
 149   fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
 150   fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
 151   fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
 152   fListOfHistosQA->Add(fEventCounter);
 153
 154   fHistQA[0] = new TH1D("fHistQAvx", "Histo Vx After Cut", 400, -4., 4.);
 155   fHistQA[1] = new TH1D("fHistQAvy", "Histo Vy After Cut", 400, -4., 4.);
 156   fHistQA[2] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
 157   fHistQA[3] = new TH1D("fHistQAvxA", "Histo Vx All", 500, -5., 5.);
 158   fHistQA[4] = new TH1D("fHistQAvyA", "Histo Vy All", 500, -5., 5.);
 159   fHistQA[5] = new TH1D("fHistQAvzA", "Histo Vz All", 500, -25., 25.);
 160   fHistQA[6] = new TH1D("fHistQADcaXyC", "Histo DCAxy after cut", 500, -5., 5.);
 161   fHistQA[7] = new TH1D("fHistQADcaZC", "Histo DCAz after cut", 500, -5., 5.);
 162   fHistQA[8] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
 163   fHistQA[9] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
 164   fHistQA[10] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
 165   fHistQA[11] = new TH1D("fHistQANCls","Number of TPC cluster",200,0,200);
 166   fHistQA[12] = new TH1D("fHistQAChi2","Chi2 per NDF",100,0,10);
 167
 168   for(Int_t i = 0; i < 13; i++)
 169     {
 170       fListOfHistosQA->Add(fHistQA[i]);
 171     }
 172
 173   fHistDCA = new TH2D("fHistDCA","DCAxy Vs DCAz", 500, -5., 5., 500, -5., 5.);
 174   fTPCSig = new TH2D("fTPCSig","TPC signal",200, 0.0, 10. ,1000,0.,1000);
 175   fTPCSig->SetMarkerColor(kRed);
 176   fTPCSigA = new TH2D("fTPCSigA","TPC signal all ",200, 0.0, 10. ,1000,0.,1000);
 177   fListOfHistosQA->Add(fHistDCA);
 178   fListOfHistosQA->Add(fTPCSig);
 179   fListOfHistosQA->Add(fTPCSigA);
 180
 181   const Int_t nDim = 3;
 182   const Int_t nPid = 5;
 183
 184   TString Species[nPid] = {"Electron","Muon","Pion","Kaon","Proton"};
 185
 186   Int_t fBins[nPid][nDim];
 187   Double_t fMin[nPid][nDim];
 188   Double_t fMax[nPid][nDim];
 189
 190   for( Int_t i = 0; i < nPid; i++ ){
 191     for( Int_t j = 0; j < nDim; j++ ){
 192       fBins[i][j] = 0;
 193       fMin[i][j] = 0.;
 194       fMax[i][j] = 0.;
 195     }
 196   }
 197
 198
 199   // if(fAnalysisType == "AOD" || fAnalysisType == "MCAOD")
 200
 201   //{
 202       TString hname1;
 203       TString htitle1, axisTitle1, axisTitle2;
 204
 205
 206
 207       if( fUsePid ){
 208
 209         Int_t gPid = (Int_t)fParticleSpecies;
 210
 211         if( gPid > 1 && gPid < 5 ){
 212           //Pion----
 213           fBins[2][0] = 100; fBins[2][1] = 1000; fBins[2][2] = 1000;
 214           fMin[2][0] = -0.5; fMin[2][1] = -0.5; fMin[2][2] = -0.5;
 215           fMax[2][0] = 99.5; fMax[2][1] = 999.5; fMax[2][2] = 999.5;
 216           //Kaon------
 217           fBins[3][0] = 100; fBins[3][1] = 600; fBins[3][2] = 600;
 218           fMin[3][0] = -0.5; fMin[3][1] = -0.5; fMin[3][2] = -0.5;
 219           fMax[3][0] = 99.5; fMax[3][1] = 599.5; fMax[3][2] = 599.5;
 220           //Proton-----
 221           fBins[4][0] = 100; fBins[4][1] = 400; fBins[4][2] = 400;
 222           fMin[4][0] = -0.5; fMin[4][1] = -0.5; fMin[4][2] = -0.5;
 223           fMax[4][0] = 99.5; fMax[4][1] = 399.5; fMax[4][2] = 399.5;
 224
 225           hname1 = "fCentNplusNminusPid"; hname1 +=gPid;
 226           htitle1 = Species[gPid]; htitle1 +=" And Neg-"; htitle1 +=Species[gPid];
 227           axisTitle1 = Species[gPid];
 228           axisTitle2 ="Neg-"; axisTitle2 += Species[gPid];
 229
 230           fTHnCentNplusNminusPid[gPid] = new THnSparseD(hname1.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
 231
 232           fTHnCentNplusNminusPid[gPid]->GetAxis(0)->SetTitle("Centrality");
 233           fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
 234           fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle2.Data());
 235
 236           fListOfHistos->Add(fTHnCentNplusNminusPid[gPid]);
 237
 238         }//Pion, Koan and Proton-------
 239
 240         else{
 241
 242           Int_t fBinsX[nDim] = {100, 1500, 1500};
 243           Double_t fMinX[nDim] = { -0.5, -0.5, -0.5 };
 244           Double_t fMaxX[nDim] = { 99.5, 1499.5, 1499.5};
 245
 246           fTHnCentNplusNminus = new THnSparseD("fTHnCentNplusNminus","Cent-NplusChrg-NminusChrg", nDim, fBinsX, fMinX, fMaxX);
 247           fTHnCentNplusNminus->GetAxis(0)->SetTitle("Centrality");
 248           fTHnCentNplusNminus->GetAxis(1)->SetTitle("UnwantedPlus");
 249           fTHnCentNplusNminus->GetAxis(2)->SetTitle("UnwantedMinus");
 250           fListOfHistos->Add(fTHnCentNplusNminus);
 251
 252         }//Unwanted particle -------
 253
 254       }//fUsePid-------
 255
 256       else{
 257
 258         Int_t fBinsCh[nDim] = {100, 1500, 1500};
 259         Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 };
 260         Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5};
 261
 262         fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
 263         fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality");
 264         fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus");
 265         fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus");
 266         fListOfHistos->Add(fTHnCentNplusNminusCh);
 267       }//All charge hadrons---------
 268
 269
 270
 271   if( fAnalysisType == "MCAOD" ){
 272
 273     Double_t dCent[2] = {-0.5, 9.5};
 274     Int_t iCent       = 10;
 275
 276     Double_t dEta[2]  = {-0.9, 0.9};
 277     Int_t iEta        = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ;
 278
 279     Double_t dRap[2]  = {-0.5, 0.5};
 280     Int_t iRap        = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ;
 281
 282     Double_t dPhi[2]  = {0.0, TMath::TwoPi()};
 283     Int_t iPhi        = 42;
 284
 285     Double_t dPt[2]   = {0.1, 3.0};
 286     Int_t iPt         = Int_t((dPt[1]-dPt[0]) / 0.075);
 287
 288     Double_t dSign[2] = {-1.5, 1.5};
 289     Int_t iSign       = 3;
 290
 291     if(fEff){
 292
 293       Int_t nBinEff[15] = { iCent, iEta, iRap, iPhi, iPt, iSign, 2, 2, 2, iEta, iPhi, iPt, iEta, 2*iPhi+1, 2*iPt+1};
 294       Double_t nMinEff[15] = { dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0], -0.5, -0.5, -0.5, dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1] };
 295       Double_t nMaxEff[15] = { dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1], 1.5, 1.5, 1.5, dEta[1], dPhi[1], dPt[1], dEta[1], dPhi[1], dPt[1] };
 296
 297       fTHnEff = new THnSparseF("fHnEff", "cent:etaMC:yMC:phiMC:ptMC:sign:findable:recStatus:pidStatus:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt", 15, nBinEff, nMinEff, nMaxEff);
 298
 299       fTHnEff->Sumw2();
 300
 301       fTHnEff->GetAxis(0)->SetTitle("centrality");
 302       fTHnEff->GetAxis(1)->SetTitle("#eta_{MC}");
 303       fTHnEff->GetAxis(2)->SetTitle("#it{y}_{MC}");
 304       fTHnEff->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");
 305       fTHnEff->GetAxis(4)->SetTitle("p_{T,MC} (GeV/#it{c})");
 306       fTHnEff->GetAxis(5)->SetTitle("sign");
 307       fTHnEff->GetAxis(6)->SetTitle("findable");
 308       fTHnEff->GetAxis(7)->SetTitle("recStatus");
 309       fTHnEff->GetAxis(8)->SetTitle("recPid");
 310       fTHnEff->GetAxis(9)->SetTitle("#eta_{Rec}");
 311       fTHnEff->GetAxis(10)->SetTitle("#varphi_{Rec} (rad)");
 312       fTHnEff->GetAxis(11)->SetTitle("p_{T,Rec} (GeV/#it{c})");
 313       fTHnEff->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}");
 314       fTHnEff->GetAxis(13)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");
 315       fTHnEff->GetAxis(14)->SetTitle("p_{T,MC}-p_{T,Rec} (GeV/#it{c})");
 316
 317       fListOfHistos->Add(fTHnEff);
 318     }
 319
 320     if(fCheckCont){
 321
 322       Int_t nBinCont[14] = { iCent,    iEta,    iRap,    iPhi,    iPt,    iSign,       8,    iSign,    iEta,    iPhi,    iPt,    iEta, 2*iPhi+1, 2*iPt+1 };
 323       Double_t nMinCont[14] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0],     0.5, dSign[0], dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
 324       Double_t nMaxCont[14] = { dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1],     8.5, dSign[1], dEta[1], dPhi[1], dPt[1], dEta[1],  dPhi[1],  dPt[1] };
 325
 326       fTHnCont = new THnSparseF("fHnCont", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",14, nBinCont, nMinCont, nMaxCont);
 327       fTHnCont->Sumw2();
 328
 329       fTHnCont->GetAxis(0)->SetTitle("centrality");
 330       fTHnCont->GetAxis(1)->SetTitle("#eta_{MC}");
 331       fTHnCont->GetAxis(2)->SetTitle("#it{y}_{MC}");
 332       fTHnCont->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");
 333       fTHnCont->GetAxis(4)->SetTitle("p_{T,MC} (GeV/#it{c})");
 334       fTHnCont->GetAxis(5)->SetTitle("sign");
 335       fTHnCont->GetAxis(6)->SetTitle("contPart");
 336       fTHnCont->GetAxis(7)->SetTitle("contSign");
 337       fTHnCont->GetAxis(8)->SetTitle("#eta_{Rec}");
 338       fTHnCont->GetAxis(9)->SetTitle("#varphi_{Rec} (rad)");
 339       fTHnCont->GetAxis(10)->SetTitle("p_{T,Rec} (GeV/#it{c})");
 340       fTHnCont->GetAxis(11)->SetTitle("#eta_{MC}-#eta_{Rec}");
 341       fTHnCont->GetAxis(12)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");
 342       fTHnCont->GetAxis(13)->SetTitle("p_{T,MC}-p_{T,Rec} (GeV/#it{c})");
 343
 344       fListOfHistos->Add(fTHnCont);
 345     }
 346
 347   }//MCAOD------
 348
 349   PostData(1, fListOfHistosQA);
 350   PostData(2, fListOfHistos);
 351
 352 }
 353
 354
 355 //----------------------------------------------------------------------------------
 356 void AliEbyEHigherMomentsEffContTask::UserExec( Option_t * ){
 357
 358   fEventCounter->Fill(1);
 359
 360
 361   if(fAnalysisType == "AOD") {
 362
 363     doAODEvent();
 364
 365   }//AOD--analysis-----
 366
 367   else if(fAnalysisType == "MCAOD") {
 368
 369     doMCAODEvent();
 370
 371   }
 372
 373   else return;
 374
 375
 376
 377
 378   fEventCounter->Fill(8);
 379
 380   PostData(1, fListOfHistosQA);
 381   PostData(2, fListOfHistos);
 382
 383 }
 384
 385 //--------------------------------------------------------------------------------------
 386 void AliEbyEHigherMomentsEffContTask::doAODEvent(){
 387
 388   //-------------------
 389   Double_t nPlusCharge = 0.;
 390   Double_t nMinusCharge = 0.;
 391   Double_t nPartile = 0.;
 392   Double_t nAntiParticle = 0.;
 393   Int_t gPid = 0.;
 394
 395
 396   //connect to the analysis mannager-----
 397
 398   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
 399   if (!manager) {
 400     cout<<"ERROR: Analysis manager not found."<<endl;
 401     return;
 402   }
 403   //coneect to the inputHandler------------
 404   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
 405   if (!inputHandler) {
 406     cout<<"ERROR: Input handler not found."<<endl;
 407     return;
 408   }
 409
 410   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 411   if (!fAOD)
 412     {
 413       cout<< "ERROR 01: AOD not found " <<endl;
 414       return;
 415     }
 416
 417   fPIDResponse =inputHandler->GetPIDResponse();
 418
 419
 420   if (!fPIDResponse){
 421     AliFatal("This Task needs the PID response attached to the inputHandler");
 422     return;
 423   }
 424
 425
 426   AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
 427   if(!aodHeader) AliFatal("Not a standard AOD");
 428
 429   Int_t cent = -1;
 430   cent =  aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data());
 431
 432   if (cent == 0)
 433     fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data());
 434   else if (cent == 10 || cent == -1.)
 435     fCentrality = -1;
 436   else if (cent > 0 && cent < 9)
 437     fCentrality = cent + 1;
 438
 439   if(fCentrality < 0 || fCentrality >= 10) return;
 440
 441   fEventCounter->Fill(2);
 442
 443   if(!ProperVertex()) return;
 444
 445   Int_t nTracks = fAOD->GetNumberOfTracks();
 446
 447   TExMap *trackMap = new TExMap();//Mapping matrix----
 448
 449   for(Int_t i = 0; i < nTracks; i++) {
 450
 451     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
 452
 453     if(!aodTrack) {
 454       AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
 455       continue;
 456     }
 457
 458     Double_t tpcSignalAll = aodTrack->GetTPCsignal();
 459     fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
 460
 461     Int_t gID = aodTrack->GetID();
 462
 463     if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
 464   }//1st track loop----
 465
 466  AliAODTrack* newAodTrack;
 467
 468   for( Int_t j = 0; j < nTracks; j++ )
 469     {
 470
 471       AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
 472
 473       if(!aodTrack1) {
 474         AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
 475         continue;
 476       }
 477
 478
 479       if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
 480
 481       Int_t gID = aodTrack1->GetID();
 482
 483       //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
 484       newAodTrack = gID >= 0 ? aodTrack1 : dynamic_cast<AliAODTrack*>(fAOD->GetTrack(trackMap->GetValue(-1-gID))); //Take those global track who corresponds to TPC only track
 485       if(!newAodTrack) AliFatal("Not a standard AOD");
 486
 487       Float_t dxy = 0., dz = 0.;
 488
 489       dxy = aodTrack1->DCA();
 490       dz  = aodTrack1->ZAtDCA();
 491
 492       Double_t pt = aodTrack1->Pt();
 493       Double_t eta = aodTrack1->Eta();
 494       Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
 495       Double_t chi2ndf = aodTrack1->Chi2perNDF();
 496
 497       /*
 498       if( fabs(dxy) > fDCAxy ) continue;
 499       if( fabs(dz) > fDCAz ) continue;
 500       //Extra cut on DCA---( Similar to BF Task.. )
 501       if( fDCAxy !=0 && fDCAz !=0 ){
 502         if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
 503       }
 504       */
 505       if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
 506       if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
 507       if( nclus < fTPCNClus ) continue;
 508       if( chi2ndf > fChi2perNDF ) continue;
 509
 510
 511       fHistQA[6]->Fill(dxy);
 512       fHistQA[7]->Fill(dz);
 513       fHistQA[8]->Fill(pt);
 514       fHistQA[9]->Fill(eta);
 515       fHistQA[10]->Fill(aodTrack1->Phi());
 516       fHistQA[11]->Fill(nclus);
 517       fHistQA[12]->Fill(chi2ndf);
 518       fHistDCA->Fill(dxy,dz);
 519
 520       Short_t gCharge = aodTrack1->Charge();
 521
 522       if(gCharge > 0) nPlusCharge += 1.;
 523       if(gCharge < 0) nMinusCharge += 1.;
 524
 525       if( fUsePid ) {
 526
 527         gPid = (Int_t)fParticleSpecies;
 528
 529         Double_t rap =  newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
 530         Double_t tpcSignal = newAodTrack->GetTPCsignal();
 531         //Double_t rap =  aodTrack1->Y(AliPID::ParticleMass(fParticleSpecies));
 532         //Double_t tpcSignal = aodTrack1->GetTPCsignal();
 533
 534         if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
 535
 536         fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
 537
 538         Float_t nsigmaTPCPID = -999.;
 539         //Float_t nsigmaTOFPID = -999.;
 540         //Float_t nsigmaTPCTOFPID = -999.;
 541
 542         nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
 543         //nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
 544
 545         if ( nsigmaTPCPID < fNSigmaCut  ){
 546
 547           if (gCharge > 0) nPartile +=1.;
 548           if( gCharge < 0 ) nAntiParticle +=1.;
 549
 550         }
 551       }//fUsepid-----
 552
 553     }//--------- Track Loop to select with filterbit
 554
 555
 556
 557   Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), nPlusCharge, nMinusCharge};
 558   Double_t fContainerPid[3] = { static_cast<Double_t>(fCentrality), nPartile, nAntiParticle};
 559
 560   if( fUsePid ){
 561
 562     fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
 563
 564     // cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl;
 565   }
 566   else{
 567
 568     fTHnCentNplusNminusCh->Fill(fContainerCh);
 569   }
 570
 571   //cout << "nCentrality "<< fCentrality <<", nPlusCharge="<< nPlusCharge << ", nMinusCharge=" << nMinusCharge << endl;
 572
 573   fEventCounter->Fill(7);
 574
 575   return;
 576
 577 }
 578 //--------------------------------------------------------------------------------------
 579
 580
 581 //--------------------------------------------------------------------------------------
 582 void AliEbyEHigherMomentsEffContTask::doMCAODEvent(){
 583
 584
 585   //---------
 586   Double_t nPlusCharge = 0.;
 587   Double_t nMinusCharge = 0.;
 588   Double_t nPartile = 0.;
 589   Double_t nAntiParticle = 0.;
 590   Int_t gPid = 0;
 591
 592   //Connect to Anlaysis manager------
 593   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
 594   if (!manager) {
 595     cout<<"ERROR: Analysis manager not found."<<endl;
 596     return;
 597   }
 598
 599   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
 600   if (!inputHandler) {
 601     cout<<"ERROR: Input handler not found."<<endl;
 602     return;
 603   }
 604
 605   //AOD event------
 606   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 607   if (!fAOD)
 608     {
 609       cout<< "ERROR 01: AOD not found " <<endl;
 610       return;
 611     }
 612
 613   fPIDResponse =inputHandler->GetPIDResponse();
 614
 615
 616   if (!fPIDResponse){
 617     AliFatal("This Task needs the PID response attached to the inputHandler");
 618     return;
 619   }
 620   // -- Get MC header
 621   // ------------------------------------------------------------------
 622
 623   fArrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 624   if (!fArrayMC) {
 625     AliFatal("Error: MC particles branch not found!\n");
 626   }
 627
 628   AliAODMCHeader *mcHdr=NULL;
 629   mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 630   if(!mcHdr) {
 631     printf("MC header branch not found!\n");
 632     return;
 633   }
 634
 635   AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
 636   if(!aodHeader) AliFatal("Not a standard AOD");
 637
 638   Int_t cent = -1;
 639   cent =  aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data());
 640
 641   if (cent == 0)
 642     fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data());
 643   else if (cent == 10 || cent == -1.)
 644     fCentrality = -1;
 645   else if (cent > 0 && cent < 9)
 646     fCentrality = cent + 1;
 647
 648   if( fCentrality < 0 || fCentrality >= 10) return;
 649
 650   fEventCounter->Fill(2);
 651
 652
 653
 654   if(!ProperVertex()) return;
 655
 656   Int_t nTracks = fAOD->GetNumberOfTracks();
 657
 658   //Creat Lables ---
 659   fLabel = new Int_t*[2];
 660   fLabel[0] = new Int_t[nTracks]; //All charged hadrons----
 661   fLabel[1] = new Int_t[nTracks]; //For Pid-----------
 662   //Initialize labels----
 663
 664   if(!fLabel[0]){
 665     AliError("Can't Get fLabel[0] ");
 666     return;
 667   }
 668   if(!fLabel[1]){
 669     AliError("Can't Get fLabel[1] ");
 670     return;
 671   }
 672
 673   for(Int_t i=0; i < 2; i++){
 674     for(Int_t j=0; j < nTracks; j++){
 675       fLabel[i][j] = 0;
 676     }
 677   }
 678
 679
 680   TExMap *trackMap = new TExMap();//Mapping matrix----
 681
 682   for(Int_t i = 0; i < nTracks; i++) {
 683
 684     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
 685
 686     if(!aodTrack) {
 687       AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
 688       continue;
 689     }
 690
 691     Double_t tpcSignalAll = aodTrack->GetTPCsignal();
 692     fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
 693
 694     Int_t gID = aodTrack->GetID();
 695
 696     if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
 697
 698   }//1st track loop----
 699
 700   AliAODTrack* newAodTrack;
 701
 702   for( Int_t j = 0; j < nTracks; j++ ){
 703
 704     AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
 705
 706     if(!aodTrack1) {
 707       AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
 708       continue;
 709     }
 710
 711
 712     if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
 713
 714
 715     Int_t gID = aodTrack1->GetID();
 716
 717     //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
 718
 719     newAodTrack = gID >= 0 ? aodTrack1 : dynamic_cast<AliAODTrack*>(fAOD->GetTrack(trackMap->GetValue(-1-gID))); //Take those global track who corresponds to TPC only track
 720     if(!newAodTrack) AliFatal("Not a standard AOD");
 721
 722     Float_t dxy = 0., dz = 0.;
 723
 724     dxy = aodTrack1->DCA();
 725     dz  = aodTrack1->ZAtDCA();
 726
 727     //cout << dxy << endl;
 728     Double_t pt = aodTrack1->Pt();
 729     Double_t eta = aodTrack1->Eta();
 730     Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
 731     Double_t chi2ndf = aodTrack1->Chi2perNDF();
 732
 733     /*
 734     if( fabs(dxy) > fDCAxy ) continue;
 735     if( fabs(dz) > fDCAz ) continue;
 736     //Extra cut on DCA---( Similar to BF Task.. )
 737     if( fDCAxy !=0 && fDCAz !=0 ){
 738       if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
 739     }
 740     */
 741
 742
 743     if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
 744     if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
 745     // if( nclus < fTPCNClus ) continue;
 746     //if( chi2ndf > fChi2perNDF ) continue;
 747
 748
 749     fHistQA[6]->Fill(dxy);
 750     fHistQA[7]->Fill(dz);
 751     fHistQA[8]->Fill(pt);
 752     fHistQA[9]->Fill(eta);
 753     fHistQA[10]->Fill(aodTrack1->Phi());
 754     fHistQA[11]->Fill(nclus);
 755     fHistQA[12]->Fill(chi2ndf);
 756     fHistDCA->Fill(dxy,dz);
 757
 758     Short_t gCharge = aodTrack1->Charge();
 759
 760     if( gCharge == 0 ) continue;
 761
 762     Int_t label  = TMath::Abs(aodTrack1->GetLabel());
 763     //fill the labels--------
 764     fLabel[0][j] = label;
 765
 766     if(gCharge > 0) nPlusCharge += 1.;
 767     if(gCharge < 0) nMinusCharge += 1.;
 768
 769     if( fUsePid ) {
 770
 771       gPid = (Int_t)fParticleSpecies;;
 772       Double_t rap =  newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
 773       Double_t tpcSignal = newAodTrack->GetTPCsignal();
 774
 775       if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
 776
 777       fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
 778
 779       Float_t nsigmaTPCPID = -999.;
 780       //Float_t nsigmaTOFPID = -999.;
 781       //Float_t nsigmaTPCTOFPID = -999.;
 782
 783       nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
 784       //nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
 785
 786       if ( nsigmaTPCPID < fNSigmaCut  ){
 787
 788         //Fill the labels----
 789         fLabel[1][j] = label;
 790
 791         if (gCharge > 0) nPartile +=1.;
 792         if( gCharge < 0 ) nAntiParticle +=1.;
 793
 794       }
 795     }//fUsepid-----
 796     //Check Contamination-------------------
 797
 798     if( fCheckCont ){
 799       CheckContTrackAOD(label, gCharge, j);
 800     }
 801
 802   }//--------- Track Loop to select with filterbit
 803
 804   if(fEff) FillEffSparse();
 805
 806   //cout << "Cent=" << fCentrality << " MC-PlusChrg=" << nPlusCharge << " MC-MinusChrg=" << nMinusCharge << endl;
 807
 808   //cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl;
 809
 810   Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), nPlusCharge, nMinusCharge};
 811   Double_t fContainerPid[3] = { static_cast<Double_t>(fCentrality), nPartile, nAntiParticle};
 812
 813   if( fUsePid ){
 814
 815     fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
 816
 817   }
 818   else{
 819
 820     fTHnCentNplusNminusCh->Fill(fContainerCh);
 821   }
 822
 823   fEventCounter->Fill(7);
 824
 825   return;
 826
 827 }
 828 //---------------------------------------------------------------------------------------
 829
 830 //---------------------------------------------------------------------------------------
 831 Bool_t AliEbyEHigherMomentsEffContTask::ProperVertex(){
 832
 833   Bool_t IsVtx = kFALSE;
 834
 835   const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
 836
 837   if(vertex) {
 838     Double32_t fCov[6];
 839     vertex->GetCovarianceMatrix(fCov);
 840     if(vertex->GetNContributors() > 0) {
 841       if(fCov[5] != 0) {
 842
 843         Double_t lvx = vertex->GetX();
 844         Double_t lvy = vertex->GetY();
 845         Double_t lvz = vertex->GetZ();
 846
 847         fEventCounter->Fill(5);
 848
 849         fHistQA[3]->Fill(lvx);fHistQA[4]->Fill(lvy);fHistQA[5]->Fill(lvz);
 850
 851         if(TMath::Abs(lvx) < fVxMax) {
 852           if(TMath::Abs(lvy) < fVyMax) {
 853             if(TMath::Abs(lvz) < fVzMax) {
 854
 855               fEventCounter->Fill(6);
 856               fHistQA[0]->Fill(lvx);fHistQA[1]->Fill(lvy);fHistQA[2]->Fill(lvz);
 857
 858               IsVtx = kTRUE;
 859
 860             }//Z-Vertex cut---
 861           }//Y-vertex cut--
 862         }//X-vertex cut---
 863       }//Covariance------
 864     }//Contributors check---
 865   }//If vertex-----------
 866
 867   return IsVtx;
 868 }
 869 //---------------------------------------------------------------------------------------
 870
 871 //---------------------------------------------------------------------------------------
 872 void AliEbyEHigherMomentsEffContTask::FillEffSparse(){
 873
 874   //For Efficiency-------------------------------------
 875
 876   Int_t nTracks = fAOD->GetNumberOfTracks();
 877   Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies);
 878
 879   Int_t nMCTrack = fArrayMC->GetEntriesFast();
 880
 881
 882   for (Int_t iMC = 0; iMC < nMCTrack; iMC++) {
 883
 884     AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
 885
 886     if(!partMC){
 887       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
 888       continue;
 889     }
 890
 891     /* TDatabasePDG *db = TDatabasePDG::Instance();
 892     TParticlePDG *part = 0x0;
 893
 894     if (!(part = db->GetParticle(partMC->PdgCode()))) continue;
 895     if ((part = db->GetParticle(partMC->PdgCode()))->Charge() == 0.0) continue;
 896     */
 897     if(!(partMC->PdgCode())) continue;
 898     if(partMC->Charge() == 0) continue;
 899     //-(1) Check for primary----
 900     if(!partMC->IsPhysicalPrimary()) continue;
 901
 902     //-(2) Basic track cuts--------
 903
 904     if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
 905     if (partMC->Pt() < fPtLowerLimit ||  partMC->Pt() > fPtHigherLimit) continue;
 906
 907     Double_t rap =  partMC->Y();
 908
 909     if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
 910
 911     if(TMath::Abs(partMC->GetPdgCode()) != gPdgCode) continue;
 912
 913     Short_t gCharge = partMC->Charge();
 914
 915     Float_t sign      = (gCharge < 0 ) ? -1. : 1.;
 916
 917     Float_t findable  = 1.;// Float_t(fHelper->IsParticleFindable(idxMC)); XXX
 918
 919     // (3) Initialize rec. variables-----
 920     Float_t recStatus = 0.;
 921     Float_t recPid    = 0.;
 922     // (4) Initialize track parameters------
 923     Float_t etaRec = 0.;
 924     Float_t phiRec = 0.;
 925     Float_t ptRec  = 0.;
 926
 927
 928
 929     // (5) Loop over all Labels--------
 930
 931
 932     for( Int_t iRec = 0; iRec < nTracks; iRec++ ){
 933
 934       if( iMC != fLabel[0][iRec] ) continue;
 935       else{
 936         recStatus = 1.;
 937         if( iMC == fLabel[1][iRec] )
 938           recPid = 1.;
 939
 940           AliAODTrack* aodTrack = NULL;
 941
 942           if(fAOD)
 943
 944             aodTrack = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iRec));
 945             if(!aodTrack) AliFatal("Not a standard AOD");
 946
 947           if(aodTrack){
 948
 949             etaRec = aodTrack->Eta();
 950             phiRec = aodTrack->Phi();
 951             ptRec  = aodTrack->Pt();
 952
 953           }//aodTrack
 954
 955           break;
 956       }
 957
 958
 959
 960     }//all charged track rec. check---
 961
 962
 963     //Fill ThnSparse----
 964
 965     Double_t EffContainer[15] = {static_cast<Double_t>(fCentrality), partMC->Eta(), partMC->Y(), partMC->Phi(), partMC->Pt(), static_cast<Double_t>(sign),static_cast<Double_t>(findable), static_cast<Double_t>(recStatus), static_cast<Double_t>(recPid), static_cast<Double_t>(etaRec), static_cast<Double_t>(phiRec), static_cast<Double_t>(ptRec), partMC->Eta()-etaRec, partMC->Phi()-phiRec, partMC->Pt()-ptRec};
 966
 967     fTHnEff->Fill(EffContainer);
 968
 969   }//iMC Track loop-----
 970
 971   return;
 972 }
 973 //---------------------------------------------------------------------------------------
 974
 975 //---------------------------------------------------------------------------------------
 976 void AliEbyEHigherMomentsEffContTask::CheckContTrackAOD(Int_t label, Float_t sign, Int_t idxTrack){
 977
 978   AliAODMCParticle* particle = (AliAODMCParticle*)fArrayMC->At(label);
 979
 980   if (!particle)
 981     return;
 982
 983   Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies);
 984   Bool_t isSecondaryFromWeakDecay = kFALSE;
 985   Bool_t isSecondaryFromMaterial  = kFALSE;
 986
 987   // -- Check if correctly identified
 988   if (particle->GetPdgCode() == (sign*gPdgCode)) {
 989
 990     // Check if is physical primary -> all ok
 991     if(particle->IsPhysicalPrimary())
 992       return;
 993
 994     // -- Check if secondaries from material or weak decay
 995     isSecondaryFromWeakDecay = particle->IsSecondaryFromWeakDecay();
 996     isSecondaryFromMaterial  = particle->IsSecondaryFromMaterial();
 997
 998   }
 999
1000   // -- Get MC pdg
1001   Float_t contSign = 0.;
1002
1003   // TDatabasePDG *db = TDatabasePDG::Instance();
1004   //TParticlePDG *part = 0x0;
1005   /*
1006      if(( part = db->GetParticle(particle->PdgCode()))->Charge() == 0.) contSign =  0.;
1007      else if((part = db->GetParticle(particle->PdgCode()))->Charge() <  0.) contSign = -1.;
1008      else if((part = db->GetParticle(particle->PdgCode()))->Charge() >  0.) contSign =  1.;
1009   */
1010
1011   if(particle->Charge() == 0.) contSign =  0.;
1012   if(particle->Charge() < 0. ) contSign = -1;
1013   if(particle->Charge() > 0. ) contSign = 1.;
1014
1015   // -- Get contaminating particle
1016   Float_t contPart = 0;
1017   if        (isSecondaryFromWeakDecay)                   contPart = 7; // probeParticle from WeakDecay
1018   else if   (isSecondaryFromMaterial)                    contPart = 8; // probeParticle from Material
1019   else {
1020     if      (TMath::Abs(particle->GetPdgCode()) ==  211) contPart = 1; // pion
1021     else if (TMath::Abs(particle->GetPdgCode()) ==  321) contPart = 2; // kaon
1022     else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton
1023     else if (TMath::Abs(particle->GetPdgCode()) ==   11) contPart = 4; // electron
1024     else if (TMath::Abs(particle->GetPdgCode()) ==   13) contPart = 5; // muon
1025     else                                                 contPart = 6; // other
1026   }
1027
1028   // -- Get Reconstructed values
1029   Float_t etaRec = 0.;
1030   Float_t phiRec = 0.;
1031   Float_t ptRec  = 0.;
1032
1033   AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(idxTrack));
1034
1035   if (aodTrack) {
1036     // if no track present (which should not happen)
1037     // -> pt = 0. , which is not in the looked at range
1038
1039     // -- Get Reconstructed eta/phi/pt
1040     etaRec = aodTrack->Eta();
1041     phiRec = aodTrack->Phi();
1042     ptRec  = aodTrack->Pt();
1043   }
1044
1045   // -- Fill THnSparse
1046   Double_t ContContainer[14] = { static_cast<Double_t>(fCentrality), particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), static_cast<Double_t>(sign), static_cast<Double_t>(contPart), static_cast<Double_t>(contSign), static_cast<Double_t>(etaRec), static_cast<Double_t>(phiRec), static_cast<Double_t>(ptRec), particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec };
1047
1048   fTHnCont->Fill(ContContainer);
1049
1050   return;
1051
1052 }
1053 //------------------------------------------------------------------------
1054 void AliEbyEHigherMomentsEffContTask::Terminate( Option_t * ){
1055
1056   Info("AliEbyEHigherMomentsEffContTask"," Task Successfully finished");
1057
1058 }
1059
1060 //------------------------------------------------------------------------