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 //           Last Modified: 30/01/2014: only for net-charge part           //
  23 //=========================================================================//
  24 #include "TChain.h"
  25 #include "TList.h"
  26 #include "TFile.h"
  27 #include "TTree.h"
  28 #include "TH1D.h"
  29 #include "TH2D.h"
  30 #include "TH3D.h"
  31 #include "THnSparse.h"
  32 #include "TCanvas.h"
  33
  34 #include "AliAnalysisTaskSE.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliAODEvent.h"
  37 #include "AliVTrack.h"
  38 #include "AliAODTrack.h"
  39 #include "AliAODInputHandler.h"
  40 #include "AliAODEvent.h"
  41 #include "AliStack.h"
  42 #include "AliAODMCHeader.h"
  43 #include "AliAODMCParticle.h"
  44 #include "AliMCEventHandler.h"
  45 #include "AliMCEvent.h"
  46 #include "AliStack.h"
  47 #include "AliGenHijingEventHeader.h"
  48 #include "AliGenEventHeader.h"
  49 #include "AliPID.h"
  50 #include "AliAODPid.h"
  51 #include "AliPIDResponse.h"
  52 #include "AliAODpidUtil.h"
  53 #include "AliPIDCombined.h"
  54 #include "AliHelperPID.h"
  55
  56 #include "AliEbyEHigherMomentsTask.h"
  57
  58 using std::cout;
  59 using std::endl;
  60
  61 ClassImp(AliEbyEHigherMomentsTask)
  62
  63
  64 //-----------------------------------------------------------------------
  65 AliEbyEHigherMomentsTask::AliEbyEHigherMomentsTask( const char *name )
  66 : AliAnalysisTaskSE( name ),
  67   fListOfHistos(0),
  68   fArrayMC(0),
  69   fAnalysisType(0),
  70   fCentralityEstimator("V0M"),
  71   fCentrality(0),
  72   fVxMax(3.),
  73   fVyMax(3.),
  74   fVzMax(10.),
  75   fPtLowerLimit(0.3),
  76   fPtHigherLimit(1.5),
  77   fNptBins(0),
  78   fBin(0),
  79   fEtaLowerLimit(-0.8),
  80   fEtaHigherLimit(0.8),
  81   fAODtrackCutBit(128),
  82   fEventCounter(0),
  83   fHistVxVy(0),
  84   fTHnCentNplusNminusCh(0),
  85   fTHnCentNplusNminusChTruth(0),
  86   fPtBinNplusNminusCh(0),
  87   fPtBinNplusNminusChTruth(0)
  88 {
  89
  90   for ( Int_t i = 0; i < 4; i++) {
  91     fHistQA[i] = NULL;
  92   }
  93
  94   DefineOutput(1, TList::Class()); // Outpput....
  95   //DefineOutput(2, TList::Class());
  96 }
  97
  98 AliEbyEHigherMomentsTask::~AliEbyEHigherMomentsTask() {
  99
 100   if(fListOfHistos) delete fListOfHistos;
 101 }
 102
 103 //---------------------------------------------------------------------------------
 104 void AliEbyEHigherMomentsTask::UserCreateOutputObjects() {
 105
 106   fListOfHistos = new TList();
 107   fListOfHistos->SetOwner(kTRUE);
 108
 109   fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
 110   fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
 111   fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-80% centrality");
 112   fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
 113   fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
 114   fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
 115   fListOfHistos->Add(fEventCounter);
 116
 117   //For QA-Histograms
 118   fHistQA[0] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
 119   fHistQA[1] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
 120   fHistQA[2] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
 121   fHistQA[3] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
 122
 123   for(Int_t i = 0; i < 4; i++)
 124     {
 125       fListOfHistos->Add(fHistQA[i]);
 126     }
 127
 128   fHistVxVy = new TH2D("fHistVxVy","Vertex-x Vs Vertex-y", 200, -1., 1., 200, -1., 1.);
 129   fListOfHistos->Add(fHistVxVy);
 130
 131
 132   const Int_t nDim = 3;
 133
 134   Int_t fBinsCh[nDim] = {100, 1500, 1500};
 135   Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 };
 136   Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5};
 137
 138   const Int_t dim = fNptBins*2 + 1;
 139   //const Int_t dim = ;
 140   Int_t bin[dim];
 141   Double_t min[dim];
 142   Double_t max[dim];
 143   bin[0] = 100; min[0] = -0.5; max[0] = 99.5;
 144
 145   for(Int_t i = 1; i < dim; i++){
 146
 147     bin[i] = 900;
 148     min[i] = -0.5;
 149     max[i] = 899.5;
 150
 151   }
 152
 153   fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
 154   fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality");
 155   fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus");
 156   fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus");
 157   fListOfHistos->Add(fTHnCentNplusNminusCh);
 158
 159   fPtBinNplusNminusCh = new THnSparseI("fPtBinNplusNminusCh","cent-nplus-nminus", dim, bin, min, max);
 160
 161   fListOfHistos->Add(fPtBinNplusNminusCh);
 162
 163   if( fAnalysisType == "MCAOD" ){
 164
 165     fTHnCentNplusNminusChTruth = new THnSparseD("fTHnCentNplusNminusChTruth","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
 166     fTHnCentNplusNminusChTruth->GetAxis(0)->SetTitle("Centrality");
 167     fTHnCentNplusNminusChTruth->GetAxis(1)->SetTitle("Nplus");
 168     fTHnCentNplusNminusChTruth->GetAxis(2)->SetTitle("Nminus");
 169     fListOfHistos->Add(fTHnCentNplusNminusChTruth);
 170
 171     fPtBinNplusNminusChTruth = new THnSparseI("fPtBinNplusNminusChTruth","cent-nplus-nminus", dim, bin, min, max);
 172     fListOfHistos->Add(fPtBinNplusNminusChTruth);
 173
 174   }//MCAOD---
 175
 176
 177
 178   //PostData(1, fListOfHistosQA);
 179   PostData(1, fListOfHistos);
 180
 181
 182 }
 183
 184
 185 //----------------------------------------------------------------------------------
 186 void AliEbyEHigherMomentsTask::UserExec( Option_t * ){
 187
 188   fEventCounter->Fill(1);
 189
 190   if(fAnalysisType == "AOD") {
 191
 192     doAODEvent();
 193
 194   }//AOD--analysis-----
 195
 196   else if(fAnalysisType == "MCAOD") {
 197
 198     doMCAODEvent();
 199
 200   }
 201
 202   else return;
 203
 204
 205
 206 }
 207
 208 //--------------------------------------------------------------------------------------
 209 void AliEbyEHigherMomentsTask::doAODEvent(){
 210
 211   Double_t positiveSum = 0.;
 212   Double_t negativeSum = 0.;
 213   const Int_t dim = fNptBins*2;
 214   Double_t PtCh[dim];
 215
 216   for(Int_t idx = 0; idx < dim; idx++){
 217     PtCh[idx] = 0.;
 218   }
 219
 220
 221   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
 222   if (!manager) {
 223     cout<<"ERROR: Analysis manager not found."<<endl;
 224     return;
 225   }
 226   //coneect to the inputHandler------------
 227   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
 228   if (!inputHandler) {
 229     cout<<"ERROR: Input handler not found."<<endl;
 230     return;
 231   }
 232
 233   AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 234
 235   if (!fAOD)
 236     {
 237       cout<< "ERROR: AOD not found " <<endl;
 238       return;
 239     }
 240
 241   if(!ProperVertex(fAOD)) return;
 242
 243   AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
 244   if(!aodHeader) AliFatal("Not a standard AOD");
 245
 246   fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 247
 248   if(fCentrality < 0 || fCentrality >= 81) return;
 249
 250   fEventCounter->Fill(2);
 251
 252   Int_t nTracks = fAOD->GetNumberOfTracks();
 253
 254   for(Int_t i = 0; i < nTracks; i++) {
 255
 256     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
 257
 258     if(!aodTrack) {
 259       AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
 260       continue;
 261     }
 262
 263     if(!AcceptTrack(aodTrack) ) continue;
 264
 265
 266     fBin = GetPtBin(aodTrack->Pt());
 267
 268     Short_t gCharge = aodTrack->Charge();
 269     if(gCharge > 0){
 270       positiveSum += 1.;
 271       PtCh[fBin] += 1;
 272     }
 273     if(gCharge < 0){
 274       negativeSum += 1.;
 275       PtCh[fNptBins+fBin] += 1;
 276     }
 277
 278
 279   }//--------- Track Loop to select with filterbit
 280
 281   Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), positiveSum, negativeSum};
 282
 283   fTHnCentNplusNminusCh->Fill(fContainerCh);
 284
 285   //cout << fCentrality << " " << positiveSum <<" " << negativeSum << endl;
 286
 287   Double_t ptContainer[dim+1];
 288
 289   ptContainer[0] = fCentrality;
 290
 291   for(Int_t i = 1; i <= dim; i++){
 292     ptContainer[i] = PtCh[i-1];
 293     //cout << PtCh[i-1] <<" ,";
 294   }
 295   //cout << endl;
 296
 297   fPtBinNplusNminusCh->Fill(ptContainer);
 298
 299
 300   fEventCounter->Fill(7);
 301   PostData(1, fListOfHistos);
 302   return;
 303
 304 }
 305 //--------------------------------------------------------------------------------------
 306 //--------------------------------------------------------------------------------------
 307 void AliEbyEHigherMomentsTask::doMCAODEvent(){
 308
 309
 310   //---------
 311   Double_t positiveSumMCRec = 0.;
 312   Double_t negativeSumMCRec = 0.;
 313
 314   Double_t positiveSumMCTruth = 0.;
 315   Double_t negativeSumMCTruth = 0.;
 316
 317   const Int_t dim = fNptBins*2;
 318   Double_t PtCh[dim];
 319   Double_t ptChMC[dim];
 320   for(Int_t idx = 0; idx < dim; idx++){
 321     PtCh[idx] = 0;
 322     ptChMC[idx] = 0;
 323   }
 324
 325   //Connect to Anlaysis manager------
 326   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
 327   if (!manager) {
 328     cout<<"ERROR: Analysis manager not found."<<endl;
 329     return;
 330   }
 331
 332   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
 333   if (!inputHandler) {
 334     cout<<"ERROR: Input handler not found."<<endl;
 335     return;
 336   }
 337
 338   //AOD event------
 339   AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 340
 341   if (!fAOD)
 342     {
 343       cout<< "ERROR: AOD not found " <<endl;
 344       return;
 345     }
 346
 347   // -- Get MC header
 348   // ------------------------------------------------------------------
 349
 350   fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 351   if (!fArrayMC) {
 352     AliFatal("Error: MC particles branch not found!\n");
 353     return;
 354   }
 355
 356   AliAODMCHeader *mcHdr=NULL;
 357   mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 358   if(!mcHdr) {
 359     printf("MC header branch not found!\n");
 360     return;
 361   }
 362
 363   if(!ProperVertex(fAOD)) return;
 364
 365   AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
 366   if(!aodHeader) AliFatal("Not a standard AOD");
 367
 368   fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 369
 370   if( fCentrality < 0 || fCentrality >= 81) return;
 371
 372   fEventCounter->Fill(2);
 373
 374   Int_t nTracks = fAOD->GetNumberOfTracks();
 375
 376   for(Int_t i = 0; i < nTracks; i++) {
 377
 378     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
 379
 380     if(!aodTrack) {
 381       AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
 382       continue;
 383     }
 384
 385     if(!AcceptTrack(aodTrack) ) continue;
 386
 387
 388     fBin = GetPtBin(aodTrack->Pt());
 389
 390     //cout << "Pt Bin " << fBin << endl;
 391
 392     Short_t gCharge = aodTrack->Charge();
 393     if(gCharge > 0){
 394       positiveSumMCRec += 1.;
 395      PtCh[fBin] += 1;
 396     }
 397     if(gCharge < 0){
 398       negativeSumMCRec += 1.;
 399       PtCh[fNptBins+fBin] += 1;
 400     }
 401
 402   }//--------- Track Loop to select with filterbit
 403
 404   //===============================================================================
 405   //---------------------MC Truth--------------------------------------------------
 406   //===============================================================================
 407
 408   Int_t nMCTrack = fArrayMC->GetEntriesFast();
 409
 410   for (Int_t iMC = 0; iMC < nMCTrack; iMC++) {
 411
 412     AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
 413
 414     if(!partMC){
 415       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
 416       continue;
 417     }
 418
 419     if(partMC->Charge() == 0) continue;
 420     if(!partMC->IsPhysicalPrimary()) continue;
 421
 422     if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
 423     if (partMC->Pt() < fPtLowerLimit ||  partMC->Pt() > fPtHigherLimit) continue;
 424
 425     Short_t gCharge = partMC->Charge();
 426     Int_t mcbin = GetPtBin(partMC->Pt());
 427
 428     if(gCharge > 0){
 429       positiveSumMCTruth += 1.;
 430       ptChMC[mcbin] += 1.;
 431     }
 432     if(gCharge < 0){
 433       negativeSumMCTruth += 1.;
 434       ptChMC[fNptBins+mcbin] += 1.;
 435     }
 436
 437
 438   }//MC-Truth Track loop--
 439
 440   Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), positiveSumMCRec, negativeSumMCRec};//Reco. values ch. hadrons
 441   Double_t fContainerChTruth[3] = { static_cast<Double_t>(fCentrality), positiveSumMCTruth, negativeSumMCTruth};
 442
 443   fTHnCentNplusNminusCh->Fill(fContainerCh);//Fill the rec. ch. particles---
 444   fTHnCentNplusNminusChTruth->Fill(fContainerChTruth);//MC -Truth ch. particles
 445
 446   //cout << fCentrality << " " << positiveSumMCRec << " " << negativeSumMCRec <<endl;
 447   //cout <<fCentrality<<"   " << positiveSumMCTruth << " " << negativeSumMCTruth << endl;
 448   //cout <<"   " << posPidSumMCRec << " " << negPidSumMCRec << endl;
 449   //cout <<"   " << posPidSumMCTruth << " " << negPidSumMCTruth << endl;
 450   //cout <<"---------------------------------" << endl;
 451
 452   Double_t ptContainer[dim+1];
 453   Double_t ptContainerMC[dim+1];
 454   ptContainer[0] = fCentrality;
 455   ptContainerMC[0] = fCentrality;
 456
 457   for(Int_t i = 1; i <= dim; i++){
 458     ptContainer[i] = PtCh[i-1];
 459     ptContainerMC[i] = ptChMC[i-1];
 460     //cout <<" "<< PtCh[i-1]<<endl;
 461     //  cout<< " MC=" << ptChMC[i-1];
 462   }
 463
 464   //cout << endl;
 465
 466   fPtBinNplusNminusCh->Fill(ptContainer);
 467   fPtBinNplusNminusChTruth->Fill(ptContainerMC);
 468
 469   fEventCounter->Fill(7);
 470   PostData(1, fListOfHistos);
 471   return;
 472
 473 }
 474
 475 //---------------------------------------------------------------------------------------
 476 Bool_t AliEbyEHigherMomentsTask::ProperVertex(AliAODEvent *fAOD) const{
 477
 478   Bool_t IsVtx = kFALSE;
 479
 480   const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
 481
 482   if(vertex) {
 483     Double32_t fCov[6];
 484     vertex->GetCovarianceMatrix(fCov);
 485     if(vertex->GetNContributors() > 0) {
 486       if(fCov[5] != 0) {
 487
 488         Double_t lvx = vertex->GetX();
 489         Double_t lvy = vertex->GetY();
 490         Double_t lvz = vertex->GetZ();
 491
 492         fEventCounter->Fill(5);
 493
 494         if(TMath::Abs(lvx) < fVxMax) {
 495           if(TMath::Abs(lvy) < fVyMax) {
 496             if(TMath::Abs(lvz) < fVzMax) {
 497
 498               fEventCounter->Fill(6);
 499               fHistQA[0]->Fill(lvz);
 500               fHistVxVy->Fill(lvx,lvy);
 501               IsVtx = kTRUE;
 502
 503             }//Z-Vertex cut---
 504           }//Y-vertex cut--
 505         }//X-vertex cut---
 506       }//Covariance------
 507     }//Contributors check---
 508   }//If vertex-----------
 509
 510   AliCentrality *centrality = fAOD->GetCentrality();
 511   if (centrality->GetQuality() != 0) IsVtx = kFALSE;
 512
 513   return IsVtx;
 514 }
 515 //---------------------------------------------------------------------------------------
 516
 517 //---------------------------------------------------------------------------------------
 518 Bool_t AliEbyEHigherMomentsTask::AcceptTrack(AliAODTrack* track) const{
 519
 520   if(!track) return kFALSE;
 521   if( track->Charge() == 0 ) return kFALSE;
 522
 523   Double_t pt = track->Pt();
 524   Double_t eta = track->Eta();
 525   if(!track->TestFilterBit(fAODtrackCutBit)) return kFALSE;
 526   if( pt < fPtLowerLimit || pt > fPtHigherLimit ) return kFALSE;
 527   if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) return kFALSE;
 528
 529   fHistQA[1]->Fill(pt);
 530   fHistQA[2]->Fill(eta);
 531   fHistQA[3]->Fill(track->Phi());
 532
 533   return kTRUE;
 534 }
 535 //------------------------------------------------------------------------
 536
 537 //------------------------------------------------------------------------
 538 Int_t AliEbyEHigherMomentsTask::GetPtBin(Double_t pt){
 539
 540   Int_t bin = 0;
 541
 542   Double_t BinSize = (fPtHigherLimit - fPtLowerLimit)/fNptBins;
 543
 544   for(Int_t iBin = 0; iBin < fNptBins; iBin++){
 545
 546     Double_t xLow = fPtLowerLimit + iBin*BinSize;
 547     Double_t xHigh = fPtLowerLimit + (iBin+1)*BinSize;
 548
 549     if( pt >= xLow && pt < xHigh){
 550       bin = iBin;
 551       //cout << "Pt Bin #" << bin <<"pt=" << pt << endl;
 552       break;
 553     }
 554
 555   }//for
 556
 557
 558   return bin;
 559 }
 560 //------------------------------------------------------------------------
 561 //------------------------------------------------------------------------
 562 void AliEbyEHigherMomentsTask::Terminate( Option_t * ){
 563
 564   Info("AliEbyEHigherMomentTask"," Task Successfully finished");
 565
 566 }
 567
 568 //------------------------------------------------------------------------