PWGCF/EBYE/Fluctuations/AliHigherMomentsToyModel.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Nirbhay K. Behera                                              *
   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 //           Toy Model 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 "AliHigherMomentsToyModel.h"
  66
  67 using std::cout;
  68 using std::endl;
  69
  70
  71 ClassImp(AliHigherMomentsToyModel)
  72
  73
  74 //-----------------------------------------------------------------------
  75 AliHigherMomentsToyModel::AliHigherMomentsToyModel( 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   fUsePid(kFALSE),
 101   fEventCounter(0),
 102   fHistDCA(0),
 103   fTPCSig(0),
 104   fTPCSigA(0),
 105   fTHnCentNplusNminusCh(0),
 106   fTHnCentNplusNminusChTruth(0),
 107   fTHnCentNplusNminus(0)
 108 {
 109
 110   for ( Int_t i = 0; i < 13; i++) {
 111     fHistQA[i] = NULL;
 112   }
 113
 114   for ( Int_t i = 0; i < 5; i++ ){
 115     fTHnCentNplusNminusPid[i] = NULL;
 116     fTHnCentNplusNminusPidTruth[i] = NULL;
 117   }
 118
 119   DefineOutput(1, TList::Class()); // Outpput....
 120   DefineOutput(2, TList::Class());
 121
 122 }
 123
 124 AliHigherMomentsToyModel::~AliHigherMomentsToyModel() {
 125
 126   if(fListOfHistosQA) delete fListOfHistosQA;
 127   if(fListOfHistos) delete fListOfHistos;
 128
 129 }
 130
 131 //---------------------------------------------------------------------------------
 132 void AliHigherMomentsToyModel::UserCreateOutputObjects() {
 133   //For QA-Histograms
 134   fListOfHistosQA = new TList();
 135   fListOfHistosQA->SetOwner(kTRUE);
 136
 137   fListOfHistos = new TList();
 138   fListOfHistos->SetOwner(kTRUE);
 139
 140   fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
 141   fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
 142   fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-90% centrality");
 143   fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
 144   fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
 145   fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
 146   fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
 147   fListOfHistosQA->Add(fEventCounter);
 148
 149   fHistQA[0] = new TH1D("fHistQAvx", "Histo Vx After Cut", 400, -4., 4.);
 150   fHistQA[1] = new TH1D("fHistQAvy", "Histo Vy After Cut", 400, -4., 4.);
 151   fHistQA[2] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
 152   fHistQA[3] = new TH1D("fHistQAvxA", "Histo Vx All", 500, -5., 5.);
 153   fHistQA[4] = new TH1D("fHistQAvyA", "Histo Vy All", 500, -5., 5.);
 154   fHistQA[5] = new TH1D("fHistQAvzA", "Histo Vz All", 500, -25., 25.);
 155   fHistQA[6] = new TH1D("fHistQADcaXyC", "Histo DCAxy after cut", 500, -5., 5.);
 156   fHistQA[7] = new TH1D("fHistQADcaZC", "Histo DCAz after cut", 500, -5., 5.);
 157   fHistQA[8] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
 158   fHistQA[9] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
 159   fHistQA[10] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
 160   fHistQA[11] = new TH1D("fHistQANCls","Number of TPC cluster",200,0,200);
 161   fHistQA[12] = new TH1D("fHistQAChi2","Chi2 per NDF",100,0,10);
 162
 163   for(Int_t i = 0; i < 13; i++)
 164     {
 165       fListOfHistosQA->Add(fHistQA[i]);
 166     }
 167
 168   fHistDCA = new TH2D("fHistDCA","DCAxy Vs DCAz", 500, -5., 5., 500, -5., 5.);
 169   fTPCSig = new TH2D("fTPCSig","TPC signal",200, 0.0, 10. ,1000,0.,1000);
 170   fTPCSig->SetMarkerColor(kRed);
 171   fTPCSigA = new TH2D("fTPCSigA","TPC signal all ",200, 0.0, 10. ,1000,0.,1000);
 172   fListOfHistosQA->Add(fHistDCA);
 173   fListOfHistosQA->Add(fTPCSig);
 174   fListOfHistosQA->Add(fTPCSigA);
 175
 176
 177   const Int_t nDim = 3;
 178   const Int_t nPid = 5;
 179
 180   TString Species[nPid] = {"Electron","Muon","Pion","Kaon","Proton"};
 181
 182   Int_t fBins[nPid][nDim];
 183   Double_t fMin[nPid][nDim];
 184   Double_t fMax[nPid][nDim];
 185
 186   for( Int_t i = 0; i < nPid; i++ ){
 187     for( Int_t j = 0; j < nDim; j++ ){
 188       fBins[i][j] = 0;
 189       fMin[i][j] = 0.;
 190       fMax[i][j] = 0.;
 191     }
 192   }
 193
 194   Int_t fBinsCh[nDim] = {100, 1500, 1500};
 195   Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 };
 196   Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5};
 197
 198   fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
 199   fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality");
 200   fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus");
 201   fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus");
 202   fListOfHistos->Add(fTHnCentNplusNminusCh);
 203
 204   if( fAnalysisType == "MCAOD" ){
 205
 206     fTHnCentNplusNminusChTruth = new THnSparseD("fTHnCentNplusNminusChTruth","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
 207     fTHnCentNplusNminusChTruth->GetAxis(0)->SetTitle("Centrality");
 208     fTHnCentNplusNminusChTruth->GetAxis(1)->SetTitle("Nplus");
 209     fTHnCentNplusNminusChTruth->GetAxis(2)->SetTitle("Nminus");
 210     fListOfHistos->Add(fTHnCentNplusNminusChTruth);
 211   }//MCAOD---
 212
 213
 214   TString hname1, hname11;
 215   TString htitle1, axisTitle1, axisTitle2;
 216
 217
 218
 219   if( fUsePid ){
 220
 221     Int_t gPid = (Int_t)fParticleSpecies;
 222
 223     if( gPid > 1 && gPid < 5 ){
 224       //Pion----
 225       fBins[2][0] = 100; fBins[2][1] = 1000; fBins[2][2] = 1000;
 226       fMin[2][0] = -0.5; fMin[2][1] = -0.5; fMin[2][2] = -0.5;
 227       fMax[2][0] = 99.5; fMax[2][1] = 999.5; fMax[2][2] = 999.5;
 228       //Kaon------
 229       fBins[3][0] = 100; fBins[3][1] = 600; fBins[3][2] = 600;
 230       fMin[3][0] = -0.5; fMin[3][1] = -0.5; fMin[3][2] = -0.5;
 231       fMax[3][0] = 99.5; fMax[3][1] = 599.5; fMax[3][2] = 599.5;
 232       //Proton-----
 233       fBins[4][0] = 100; fBins[4][1] = 400; fBins[4][2] = 400;
 234       fMin[4][0] = -0.5; fMin[4][1] = -0.5; fMin[4][2] = -0.5;
 235       fMax[4][0] = 99.5; fMax[4][1] = 399.5; fMax[4][2] = 399.5;
 236
 237       hname1 = "fCentNplusNminusPid"; hname1 +=gPid;
 238       htitle1 = Species[gPid]; htitle1 +=" And Neg-"; htitle1 +=Species[gPid];
 239       axisTitle1 = Species[gPid];
 240       axisTitle2 ="Neg-"; axisTitle2 += Species[gPid];
 241
 242       fTHnCentNplusNminusPid[gPid] = new THnSparseD(hname1.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
 243
 244       fTHnCentNplusNminusPid[gPid]->GetAxis(0)->SetTitle("Centrality");
 245       fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
 246       fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle2.Data());
 247
 248       fListOfHistos->Add(fTHnCentNplusNminusPid[gPid]);
 249
 250       if( fAnalysisType == "MCAOD" ){
 251
 252         hname11 = "fCentNplusNminusPidTruth"; hname11 +=gPid;
 253         fTHnCentNplusNminusPidTruth[gPid] = new THnSparseD(hname11.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
 254
 255         fTHnCentNplusNminusPidTruth[gPid]->GetAxis(0)->SetTitle("Centrality");
 256         fTHnCentNplusNminusPidTruth[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
 257         fTHnCentNplusNminusPidTruth[gPid]->GetAxis(1)->SetTitle(axisTitle2.Data());
 258
 259         fListOfHistos->Add(fTHnCentNplusNminusPidTruth[gPid]);
 260
 261       }//MCAOD-----
 262     }//Pion, Koan and Proton-------
 263
 264     else{
 265
 266       Int_t fBinsX[nDim] = {100, 1500, 1500};
 267       Double_t fMinX[nDim] = { -0.5, -0.5, -0.5 };
 268       Double_t fMaxX[nDim] = { 99.5, 1499.5, 1499.5};
 269
 270       fTHnCentNplusNminus = new THnSparseD("fTHnCentNplusNminus","Cent-NplusChrg-NminusChrg", nDim, fBinsX, fMinX, fMaxX);
 271       fTHnCentNplusNminus->GetAxis(0)->SetTitle("Centrality");
 272       fTHnCentNplusNminus->GetAxis(1)->SetTitle("UnwantedPlus");
 273       fTHnCentNplusNminus->GetAxis(2)->SetTitle("UnwantedMinus");
 274       fListOfHistos->Add(fTHnCentNplusNminus);
 275
 276     }//Unwanted particle -------
 277
 278   }//fUsePid-------
 279
 280
 281
 282   PostData(1, fListOfHistosQA);
 283   PostData(2, fListOfHistos);
 284
 285 }
 286
 287
 288 //----------------------------------------------------------------------------------
 289 void AliHigherMomentsToyModel::UserExec( Option_t * ){
 290
 291   fEventCounter->Fill(1);
 292
 293
 294   if(fAnalysisType == "AOD") {
 295
 296     doAODEvent();
 297
 298   }//AOD--analysis-----
 299
 300   else if(fAnalysisType == "MCAOD") {
 301
 302     doMCAODEvent();
 303
 304   }
 305
 306   else return;
 307
 308
 309
 310
 311   fEventCounter->Fill(8);
 312
 313   PostData(1, fListOfHistosQA);
 314   PostData(2, fListOfHistos);
 315
 316 }
 317
 318 //--------------------------------------------------------------------------------------
 319 void AliHigherMomentsToyModel::doAODEvent(){
 320
 321   //-------------------
 322   Double_t nPlusCharge = 0.;
 323   Double_t nMinusCharge = 0.;
 324   Double_t nPartile = 0.;
 325   Double_t nAntiParticle = 0.;
 326   Int_t gPid = 0.;
 327
 328
 329   //connect to the analysis mannager-----
 330
 331   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
 332   if (!manager) {
 333     cout<<"ERROR: Analysis manager not found."<<endl;
 334     return;
 335   }
 336   //coneect to the inputHandler------------
 337   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
 338   if (!inputHandler) {
 339     cout<<"ERROR: Input handler not found."<<endl;
 340     return;
 341   }
 342
 343   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 344   if (!fAOD)
 345     {
 346       cout<< "ERROR 01: AOD not found " <<endl;
 347       return;
 348     }
 349
 350   fPIDResponse =inputHandler->GetPIDResponse();
 351
 352
 353   if (!fPIDResponse){
 354     AliFatal("This Task needs the PID response attached to the inputHandler");
 355     return;
 356   }
 357
 358
 359   AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
 360   if(!aodHeader) AliFatal("Not a standard AOD");
 361
 362   fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 363   /* Int_t cent = -1;
 364   cent =  aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data());
 365
 366   if (cent == 0)
 367     fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data());
 368   else if (cent == 10 || cent == -1.)
 369     fCentrality = -1;
 370   else if (cent > 0 && cent < 9)
 371     fCentrality = cent + 1;
 372   */
 373   if(fCentrality < 0 || fCentrality >= 91) return;
 374
 375   fEventCounter->Fill(2);
 376
 377   if(!ProperVertex()) return;
 378
 379   Int_t nTracks = fAOD->GetNumberOfTracks();
 380
 381   TExMap *trackMap = new TExMap();//Mapping matrix----
 382
 383   for(Int_t i = 0; i < nTracks; i++) {
 384
 385     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
 386
 387     if(!aodTrack) {
 388       AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
 389       continue;
 390     }
 391
 392     Double_t tpcSignalAll = aodTrack->GetTPCsignal();
 393     fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
 394
 395     Int_t gID = aodTrack->GetID();
 396
 397     if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
 398   }//1st track loop----
 399
 400  AliAODTrack* newAodTrack;
 401
 402   for( Int_t j = 0; j < nTracks; j++ )
 403     {
 404
 405       AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
 406
 407       if(!aodTrack1) {
 408         AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
 409         continue;
 410       }
 411
 412
 413       if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
 414
 415       Int_t gID = aodTrack1->GetID();
 416
 417       //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
 418       newAodTrack = gID >= 0 ? aodTrack1 : dynamic_cast<AliAODTrack*>(fAOD->GetTrack(trackMap->GetValue(-1-gID))); //Take those global track who corresponds to TPC only track
 419       if(!newAodTrack) AliFatal("Not a standard AOD");
 420
 421       Float_t dxy = 0., dz = 0.;
 422
 423       dxy = aodTrack1->DCA();
 424       dz  = aodTrack1->ZAtDCA();
 425
 426       Double_t pt = aodTrack1->Pt();
 427       Double_t eta = aodTrack1->Eta();
 428       Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
 429       Double_t chi2ndf = aodTrack1->Chi2perNDF();
 430
 431       /*
 432       if( fabs(dxy) > fDCAxy ) continue;
 433       if( fabs(dz) > fDCAz ) continue;
 434       //Extra cut on DCA---( Similar to BF Task.. )
 435       if( fDCAxy !=0 && fDCAz !=0 ){
 436         if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
 437       }
 438       */
 439       if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
 440       if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
 441       if( nclus < fTPCNClus ) continue;
 442       if( chi2ndf > fChi2perNDF ) continue;
 443
 444
 445       fHistQA[6]->Fill(dxy);
 446       fHistQA[7]->Fill(dz);
 447       fHistQA[8]->Fill(pt);
 448       fHistQA[9]->Fill(eta);
 449       fHistQA[10]->Fill(aodTrack1->Phi());
 450       fHistQA[11]->Fill(nclus);
 451       fHistQA[12]->Fill(chi2ndf);
 452       fHistDCA->Fill(dxy,dz);
 453
 454       Short_t gCharge = aodTrack1->Charge();
 455
 456       if(gCharge > 0) nPlusCharge += 1.;
 457       if(gCharge < 0) nMinusCharge += 1.;
 458
 459       if( fUsePid ) {
 460
 461         gPid = (Int_t)fParticleSpecies;
 462
 463         Double_t rap =  newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
 464         Double_t tpcSignal = newAodTrack->GetTPCsignal();
 465         //Double_t rap =  aodTrack1->Y(AliPID::ParticleMass(fParticleSpecies));
 466         //Double_t tpcSignal = aodTrack1->GetTPCsignal();
 467
 468         if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
 469
 470         fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
 471
 472         Float_t nsigmaTPCPID = -999.;
 473         //Float_t nsigmaTOFPID = -999.;
 474         //Float_t nsigmaTPCTOFPID = -999.;
 475
 476         nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
 477         //nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
 478
 479         if ( nsigmaTPCPID < fNSigmaCut  ){
 480
 481           if (gCharge > 0) nPartile +=1.;
 482           if( gCharge < 0 ) nAntiParticle +=1.;
 483
 484         }
 485       }//fUsepid-----
 486
 487     }//--------- Track Loop to select with filterbit
 488
 489
 490
 491   Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), nPlusCharge, nMinusCharge};
 492   Double_t fContainerPid[3] = { static_cast<Double_t>(fCentrality), nPartile, nAntiParticle};
 493
 494
 495   fTHnCentNplusNminusCh->Fill(fContainerCh);
 496
 497   if( fUsePid ){
 498     gPid = (Int_t)fParticleSpecies;
 499     fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
 500
 501     // cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl;
 502   }
 503
 504
 505   //cout << "nCentrality "<< fCentrality <<", nPlusCharge="<< nPlusCharge << ", nMinusCharge=" << nMinusCharge << endl;
 506
 507   fEventCounter->Fill(7);
 508
 509   return;
 510
 511 }
 512 //--------------------------------------------------------------------------------------
 513
 514
 515 //--------------------------------------------------------------------------------------
 516 void AliHigherMomentsToyModel::doMCAODEvent(){
 517
 518
 519   //---------
 520   Double_t nPlusCharge = 0.;
 521   Double_t nMinusCharge = 0.;
 522
 523   Double_t nPlusChargeTruth = 0.;
 524   Double_t nMinusChargeTruth = 0.;
 525
 526   Double_t nPartile = 0.;
 527   Double_t nAntiParticle = 0.;
 528   Double_t nPartileTruth = 0.;
 529   Double_t nAntiParticleTruth = 0.;
 530
 531   Int_t gPid = 0;
 532   Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies);
 533
 534   //Connect to Anlaysis manager------
 535   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
 536   if (!manager) {
 537     cout<<"ERROR: Analysis manager not found."<<endl;
 538     return;
 539   }
 540
 541   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
 542   if (!inputHandler) {
 543     cout<<"ERROR: Input handler not found."<<endl;
 544     return;
 545   }
 546
 547   //AOD event------
 548   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 549   if (!fAOD)
 550     {
 551       cout<< "ERROR 01: AOD not found " <<endl;
 552       return;
 553     }
 554
 555   fPIDResponse =inputHandler->GetPIDResponse();
 556
 557
 558   if (!fPIDResponse){
 559     AliFatal("This Task needs the PID response attached to the inputHandler");
 560     return;
 561   }
 562   // -- Get MC header
 563   // ------------------------------------------------------------------
 564
 565   fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 566   if (!fArrayMC) {
 567     AliFatal("Error: MC particles branch not found!\n");
 568     return;
 569   }
 570
 571   AliAODMCHeader *mcHdr=NULL;
 572   mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 573   if(!mcHdr) {
 574     printf("MC header branch not found!\n");
 575     return;
 576   }
 577
 578   AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
 579   if(!aodHeader) AliFatal("Not a standard AOD");
 580
 581   fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 582
 583   /*
 584     Int_t cent = -1;
 585     cent =  aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data());
 586
 587     if (cent == 0)
 588     fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data());
 589     else if (cent == 10 || cent == -1.)
 590     fCentrality = -1;
 591     else if (cent > 0 && cent < 9)
 592     fCentrality = cent + 1;
 593   */
 594   if( fCentrality < 0 || fCentrality >= 91) return;
 595
 596   fEventCounter->Fill(2);
 597
 598
 599
 600   if(!ProperVertex()) return;
 601
 602   Int_t nTracks = fAOD->GetNumberOfTracks();
 603
 604
 605   TExMap *trackMap = new TExMap();//Mapping matrix----
 606
 607   for(Int_t i = 0; i < nTracks; i++) {
 608
 609     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
 610
 611     if(!aodTrack) {
 612       AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
 613       continue;
 614     }
 615
 616     Double_t tpcSignalAll = aodTrack->GetTPCsignal();
 617     fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
 618
 619     Int_t gID = aodTrack->GetID();
 620
 621     if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
 622
 623   }//1st track loop----
 624
 625   AliAODTrack* newAodTrack;
 626
 627   for( Int_t j = 0; j < nTracks; j++ ){
 628
 629     AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
 630
 631     if(!aodTrack1) {
 632       AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
 633       continue;
 634     }
 635
 636
 637     if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
 638
 639
 640     Int_t gID = aodTrack1->GetID();
 641
 642     //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
 643
 644     newAodTrack = gID >= 0 ? aodTrack1 : dynamic_cast<AliAODTrack*>(fAOD->GetTrack(trackMap->GetValue(-1-gID))); //Take those global track who corresponds to TPC only track
 645     if(!newAodTrack) AliFatal("Not a standard AOD");
 646
 647
 648     //cout << dxy << endl;
 649     Double_t pt = aodTrack1->Pt();
 650     Double_t eta = aodTrack1->Eta();
 651     Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
 652     Double_t chi2ndf = aodTrack1->Chi2perNDF();
 653
 654     Float_t dxy = 0., dz = 0.;
 655
 656     if( fAODtrackCutBit == 128 ){
 657
 658       dxy = aodTrack1->DCA();
 659       dz  = aodTrack1->ZAtDCA();
 660
 661       if( fabs(dxy) > fDCAxy ) continue;
 662       if( fabs(dz) > fDCAz ) continue;
 663       //Extra cut on DCA---( Similar to BF Task.. )
 664       if( fDCAxy !=0 && fDCAz !=0 ){
 665         if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
 666       }
 667
 668       fHistQA[6]->Fill(dxy);
 669       fHistQA[7]->Fill(dz);
 670       fHistDCA->Fill(dxy,dz);
 671
 672     }
 673
 674
 675     if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
 676     if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
 677     if( nclus < fTPCNClus ) continue;
 678     if( chi2ndf > fChi2perNDF ) continue;
 679
 680
 681
 682     fHistQA[8]->Fill(pt);
 683     fHistQA[9]->Fill(eta);
 684     fHistQA[10]->Fill(aodTrack1->Phi());
 685     fHistQA[11]->Fill(nclus);
 686     fHistQA[12]->Fill(chi2ndf);
 687
 688
 689     Short_t gCharge = aodTrack1->Charge();
 690
 691     if( gCharge == 0 ) continue;
 692
 693
 694     if(gCharge > 0) nPlusCharge += 1.;
 695     if(gCharge < 0) nMinusCharge += 1.;
 696
 697     if( fUsePid ) {
 698
 699       gPid = (Int_t)fParticleSpecies;;
 700       Double_t rap =  newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
 701       Double_t tpcSignal = newAodTrack->GetTPCsignal();
 702
 703       if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
 704
 705       fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
 706
 707       Float_t nsigmaTPCPID = -999.;
 708       //Float_t nsigmaTOFPID = -999.;
 709       //Float_t nsigmaTPCTOFPID = -999.;
 710
 711       nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
 712       //nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
 713
 714       if ( nsigmaTPCPID < fNSigmaCut  ){
 715
 716
 717         if (gCharge > 0) nPartile +=1.;
 718         if( gCharge < 0 ) nAntiParticle +=1.;
 719
 720       }//nSigmaCut-----
 721     }//fUsepid-----
 722
 723   }//--------- Track Loop to select with filterbit
 724
 725   //cout << "Cent=" << fCentrality << " MC-PlusChrg=" << nPlusCharge << " MC-MinusChrg=" << nMinusCharge << endl;
 726
 727   //cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl;
 728
 729   Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), nPlusCharge, nMinusCharge};
 730   Double_t fContainerPid[3] = { static_cast<Double_t>(fCentrality), nPartile, nAntiParticle};
 731
 732   fTHnCentNplusNminusCh->Fill(fContainerCh);
 733
 734   if( fUsePid ){
 735
 736     gPid = (Int_t)fParticleSpecies;
 737     fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
 738
 739   }
 740
 741   //===========================================
 742   //--------MC Truth---------------------------
 743   //===========================================
 744
 745   Int_t nMCTrack = fArrayMC->GetEntriesFast();
 746
 747   for (Int_t iMC = 0; iMC < nMCTrack; iMC++) {
 748
 749     AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
 750
 751     if(!partMC){
 752       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
 753       continue;
 754     }
 755
 756     if(partMC->Charge() == 0) continue;
 757     if(!partMC->IsPhysicalPrimary()) continue;
 758
 759     if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
 760     if (partMC->Pt() < fPtLowerLimit ||  partMC->Pt() > fPtHigherLimit) continue;
 761
 762     Short_t gCharge = partMC->Charge();
 763
 764     if(gCharge > 0) nPlusChargeTruth += 1.;
 765     if(gCharge < 0) nMinusChargeTruth += 1.;
 766
 767     if(fUsePid){
 768
 769       Double_t rap =  partMC->Y();
 770
 771       if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
 772
 773       if(TMath::Abs(partMC->GetPdgCode()) != gPdgCode) continue;
 774
 775       if(gCharge > 0) nPartileTruth += 1.;
 776       if(gCharge < 0) nAntiParticleTruth += 1.;
 777
 778     }//if(fUsePid) ----
 779
 780   }//MC-Truth Track loop--
 781
 782   Double_t fContainerChTruth[3] = { static_cast<Double_t>(fCentrality), nPlusChargeTruth, nMinusChargeTruth };
 783   Double_t fContainerPidTruth[3] = { static_cast<Double_t>(fCentrality), nPartileTruth, nAntiParticleTruth };
 784
 785   //cout << "Cent=" << fCentrality << " MC-PlusChrgT=" << nPlusChargeTruth << " MC-MinusChrgT=" << nMinusChargeTruth << endl;
 786
 787   //cout << "nCentrality "<< fCentrality <<", nParticleT="<< nPartileTruth << ", nMinusParticleT=" << nAntiParticleTruth << endl;
 788
 789   //cout << "----------------------------" << endl;
 790   fTHnCentNplusNminusChTruth->Fill(fContainerChTruth);
 791
 792   if( fUsePid ){
 793     gPid = (Int_t)fParticleSpecies;
 794     fTHnCentNplusNminusPidTruth[gPid]->Fill(fContainerPidTruth);
 795
 796   }
 797
 798
 799   fEventCounter->Fill(7);
 800
 801   return;
 802
 803 }
 804 //---------------------------------------------------------------------------------------
 805
 806 //---------------------------------------------------------------------------------------
 807 Bool_t AliHigherMomentsToyModel::ProperVertex(){
 808
 809   Bool_t IsVtx = kFALSE;
 810
 811   const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
 812
 813   if(vertex) {
 814     Double32_t fCov[6];
 815     vertex->GetCovarianceMatrix(fCov);
 816     if(vertex->GetNContributors() > 0) {
 817       if(fCov[5] != 0) {
 818
 819         Double_t lvx = vertex->GetX();
 820         Double_t lvy = vertex->GetY();
 821         Double_t lvz = vertex->GetZ();
 822
 823         fEventCounter->Fill(5);
 824
 825         fHistQA[3]->Fill(lvx);fHistQA[4]->Fill(lvy);fHistQA[5]->Fill(lvz);
 826
 827         if(TMath::Abs(lvx) < fVxMax) {
 828           if(TMath::Abs(lvy) < fVyMax) {
 829             if(TMath::Abs(lvz) < fVzMax) {
 830
 831               fEventCounter->Fill(6);
 832               fHistQA[0]->Fill(lvx);fHistQA[1]->Fill(lvy);fHistQA[2]->Fill(lvz);
 833
 834               IsVtx = kTRUE;
 835
 836             }//Z-Vertex cut---
 837           }//Y-vertex cut--
 838         }//X-vertex cut---
 839       }//Covariance------
 840     }//Contributors check---
 841   }//If vertex-----------
 842
 843   return IsVtx;
 844 }
 845 //---------------------------------------------------------------------------------------
 846
 847
 848
 849 //------------------------------------------------------------------------
 850 void AliHigherMomentsToyModel::Terminate( Option_t * ){
 851
 852   Info("AliHigherMomentsToyModel"," Task Successfully finished");
 853
 854 }
 855
 856 //------------------------------------------------------------------------