[u/mrichter/AliRoot.git] / PWGCF / EBYE / Fluctuations / AliEbyEHigherMomentsTask.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: Satyajit Jena.                                                 *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


//=========================================================================//
//                                                                         //
//           Analysis Task for Net-Charge Higher Moment Analysis           //
//              Author: Satyajit Jena || Nirbhay K. Behera                 //
//                      sjena@cern.ch || nbehera@cern.ch                   //
//           Last Modified: 30/01/2014: only for net-charge part           //
//=========================================================================//
#include "TChain.h"
#include "TList.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "THnSparse.h"
#include "TCanvas.h"

#include "AliAnalysisTaskSE.h"
#include "AliAnalysisManager.h"
#include "AliAODEvent.h"
#include "AliVTrack.h"
#include "AliAODTrack.h"
#include "AliAODInputHandler.h"
#include "AliAODEvent.h"
#include "AliStack.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliGenHijingEventHeader.h"
#include "AliGenEventHeader.h"
#include "AliPID.h"
#include "AliAODPid.h"                
#include "AliPIDResponse.h"
#include "AliAODpidUtil.h"        
#include "AliPIDCombined.h"
#include "AliHelperPID.h"

#include "AliEbyEHigherMomentsTask.h"

using std::cout;
using std::endl;

ClassImp(AliEbyEHigherMomentsTask)


//-----------------------------------------------------------------------
AliEbyEHigherMomentsTask::AliEbyEHigherMomentsTask( const char *name )
: AliAnalysisTaskSE( name ),
  fListOfHistos(0),
  fArrayMC(0),
  fAnalysisType(0),
  fCentralityEstimator("V0M"),
  fCentrality(0),
  fVxMax(3.),
  fVyMax(3.),
  fVzMax(10.),
  fPtLowerLimit(0.3),
  fPtHigherLimit(1.5),
  fNptBins(0),
  fBin(0),
  fEtaLowerLimit(-0.8),
  fEtaHigherLimit(0.8),
  fAODtrackCutBit(128),
  fEventCounter(0),
  fHistVxVy(0),
  fTHnCentNplusNminusCh(0),
  fTHnCentNplusNminusChTruth(0),
  fPtBinNplusNminusCh(0),
  fPtBinNplusNminusChTruth(0)
{ 
  
  for ( Int_t i = 0; i < 4; i++) { 
    fHistQA[i] = NULL;
  }
 
  DefineOutput(1, TList::Class()); // Outpput....
  //DefineOutput(2, TList::Class()); 
}

AliEbyEHigherMomentsTask::~AliEbyEHigherMomentsTask() {

  if(fListOfHistos) delete fListOfHistos;
}

//---------------------------------------------------------------------------------
void AliEbyEHigherMomentsTask::UserCreateOutputObjects() {
  
  fListOfHistos = new TList();
  fListOfHistos->SetOwner(kTRUE);
  
  fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
  fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
  fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-80% centrality");
  fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
  fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
  fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
  fListOfHistos->Add(fEventCounter);
  
  //For QA-Histograms
  fHistQA[0] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);  
  fHistQA[1] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
  fHistQA[2] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
  fHistQA[3] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
 
  for(Int_t i = 0; i < 4; i++)
    {
      fListOfHistos->Add(fHistQA[i]);
    }
  
  fHistVxVy = new TH2D("fHistVxVy","Vertex-x Vs Vertex-y", 200, -1., 1., 200, -1., 1.);
  fListOfHistos->Add(fHistVxVy);
 
 
  const Int_t nDim = 3;
  
  Int_t fBinsCh[nDim] = {100, 1500, 1500};
  Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 };
  Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5};
  
  const Int_t dim = fNptBins*2 + 1;
  //const Int_t dim = ;
  Int_t bin[dim];
  Double_t min[dim];
  Double_t max[dim];
  bin[0] = 100; min[0] = -0.5; max[0] = 99.5;

  for(Int_t i = 1; i < dim; i++){
    
    bin[i] = 900;
    min[i] = -0.5;
    max[i] = 899.5;
 
  } 
  
  fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh); 
  fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality");
  fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus");
  fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus");
  fListOfHistos->Add(fTHnCentNplusNminusCh);

  fPtBinNplusNminusCh = new THnSparseI("fPtBinNplusNminusCh","cent-nplus-nminus", dim, bin, min, max);
  
  fListOfHistos->Add(fPtBinNplusNminusCh);
  
  if( fAnalysisType == "MCAOD" ){   
    
    fTHnCentNplusNminusChTruth = new THnSparseD("fTHnCentNplusNminusChTruth","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh); 
    fTHnCentNplusNminusChTruth->GetAxis(0)->SetTitle("Centrality");
    fTHnCentNplusNminusChTruth->GetAxis(1)->SetTitle("Nplus");
    fTHnCentNplusNminusChTruth->GetAxis(2)->SetTitle("Nminus");
    fListOfHistos->Add(fTHnCentNplusNminusChTruth);
    
    fPtBinNplusNminusChTruth = new THnSparseI("fPtBinNplusNminusChTruth","cent-nplus-nminus", dim, bin, min, max);
    fListOfHistos->Add(fPtBinNplusNminusChTruth);
    
  }//MCAOD---
  
 
  //PostData(1, fListOfHistosQA);
  PostData(1, fListOfHistos);   
  
  
}


//----------------------------------------------------------------------------------
void AliEbyEHigherMomentsTask::UserExec( Option_t * ){
  
  fEventCounter->Fill(1);
  
  if(fAnalysisType == "AOD") {
    
    doAODEvent();
    
  }//AOD--analysis-----
  
  else if(fAnalysisType == "MCAOD") {
    
    doMCAODEvent();
    
  }
  
  else return;
 
  
}

//--------------------------------------------------------------------------------------
void AliEbyEHigherMomentsTask::doAODEvent(){
  
  Double_t positiveSum = 0.;
  Double_t negativeSum = 0.;
  const Int_t dim = fNptBins*2;
  Double_t PtCh[dim];

  for(Int_t idx = 0; idx < dim; idx++){
    PtCh[idx] = 0.;
  }


  AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
  if (!manager) {
    cout<<"ERROR: Analysis manager not found."<<endl;
    return;
  }
  //coneect to the inputHandler------------
  AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
  if (!inputHandler) {
    cout<<"ERROR: Input handler not found."<<endl;
    return;
  }
  
  AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());

  if (!fAOD)
    {
      cout<< "ERROR: AOD not found " <<endl;
      return;
    }
  
  if(!ProperVertex(fAOD)) return;   
  
  AliAODHeader *aodHeader = fAOD->GetHeader();
  
  fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
  
  if(fCentrality < 0 || fCentrality >= 81) return;
  
  fEventCounter->Fill(2);
  
  Int_t nTracks = fAOD->GetNumberOfTracks();
  
  for(Int_t i = 0; i < nTracks; i++) {
    
    AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i)); 
    
    if(!aodTrack) {
      AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
      continue;
    }
    
    if(!AcceptTrack(aodTrack) ) continue;

    
    fBin = GetPtBin(aodTrack->Pt());
    
    Short_t gCharge = aodTrack->Charge();
    if(gCharge > 0){
      positiveSum += 1.;
      PtCh[fBin] += 1;
    }
    if(gCharge < 0){
      negativeSum += 1.;
      PtCh[fNptBins+fBin] += 1;
    }
   
   
  }//--------- Track Loop to select with filterbit
  
  Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), positiveSum, negativeSum};
  
  fTHnCentNplusNminusCh->Fill(fContainerCh);

  //cout << fCentrality << " " << positiveSum <<" " << negativeSum << endl;
  
  Double_t ptContainer[dim+1];

  ptContainer[0] = fCentrality;

  for(Int_t i = 1; i <= dim; i++){
    ptContainer[i] = PtCh[i-1];
    //cout << PtCh[i-1] <<" ,";
  }
  //cout << endl;

  fPtBinNplusNminusCh->Fill(ptContainer);
  
  
  fEventCounter->Fill(7);
  PostData(1, fListOfHistos);
  return;
  
}
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
void AliEbyEHigherMomentsTask::doMCAODEvent(){
  
  
  //---------
  Double_t positiveSumMCRec = 0.;
  Double_t negativeSumMCRec = 0.;
 
  Double_t positiveSumMCTruth = 0.;
  Double_t negativeSumMCTruth = 0.;
 
  const Int_t dim = fNptBins*2;
  Double_t PtCh[dim];
  Double_t ptChMC[dim];
  for(Int_t idx = 0; idx < dim; idx++){
    PtCh[idx] = 0;
    ptChMC[idx] = 0;
  }
  
  //Connect to Anlaysis manager------
  AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
  if (!manager) {
    cout<<"ERROR: Analysis manager not found."<<endl;
    return;
  }
  
  AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
  if (!inputHandler) {
    cout<<"ERROR: Input handler not found."<<endl;
    return;
  }
  
  //AOD event------
  AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
  
  if (!fAOD)
    {
      cout<< "ERROR: AOD not found " <<endl;
      return;
    }

  // -- Get MC header
  // ------------------------------------------------------------------
  
  fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
  if (!fArrayMC) {
    AliFatal("Error: MC particles branch not found!\n");
    return;
  }
  
  AliAODMCHeader *mcHdr=NULL;
  mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());  
  if(!mcHdr) {
    printf("MC header branch not found!\n");
    return;
  }
  
  if(!ProperVertex(fAOD)) return;

  AliAODHeader *aodHeader = fAOD->GetHeader();
  
  fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());

  if( fCentrality < 0 || fCentrality >= 81) return;
  
  fEventCounter->Fill(2);
 
  Int_t nTracks = fAOD->GetNumberOfTracks();
  
  for(Int_t i = 0; i < nTracks; i++) {
    
    AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i)); 
    
    if(!aodTrack) {
      AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
      continue;
    }
    
    if(!AcceptTrack(aodTrack) ) continue;

    
    fBin = GetPtBin(aodTrack->Pt());
    
    //cout << "Pt Bin " << fBin << endl;

    Short_t gCharge = aodTrack->Charge();
    if(gCharge > 0){
      positiveSumMCRec += 1.;
     PtCh[fBin] += 1;
    }
    if(gCharge < 0){
      negativeSumMCRec += 1.;    
      PtCh[fNptBins+fBin] += 1;
    }
   
  }//--------- Track Loop to select with filterbit
  
  //===============================================================================
  //---------------------MC Truth--------------------------------------------------
  //===============================================================================
  
  Int_t nMCTrack = fArrayMC->GetEntriesFast();
  
  for (Int_t iMC = 0; iMC < nMCTrack; iMC++) {
    
    AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
    
    if(!partMC){
      AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
      continue;
    }
    
    if(partMC->Charge() == 0) continue;
    if(!partMC->IsPhysicalPrimary()) continue;
    
    if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
    if (partMC->Pt() < fPtLowerLimit ||  partMC->Pt() > fPtHigherLimit) continue;
    
    Short_t gCharge = partMC->Charge();
    Int_t mcbin = GetPtBin(partMC->Pt());
    
    if(gCharge > 0){
      positiveSumMCTruth += 1.;
      ptChMC[mcbin] += 1.;
    }
    if(gCharge < 0){
      negativeSumMCTruth += 1.;
      ptChMC[fNptBins+mcbin] += 1.;
    }
   
   
  }//MC-Truth Track loop--

  Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), positiveSumMCRec, negativeSumMCRec};//Reco. values ch. hadrons
  Double_t fContainerChTruth[3] = { static_cast<Double_t>(fCentrality), positiveSumMCTruth, negativeSumMCTruth};  
 
  fTHnCentNplusNminusCh->Fill(fContainerCh);//Fill the rec. ch. particles---
  fTHnCentNplusNminusChTruth->Fill(fContainerChTruth);//MC -Truth ch. particles
  
  //cout << fCentrality << " " << positiveSumMCRec << " " << negativeSumMCRec <<endl;
  //cout <<fCentrality<<"   " << positiveSumMCTruth << " " << negativeSumMCTruth << endl;
  //cout <<"   " << posPidSumMCRec << " " << negPidSumMCRec << endl;
  //cout <<"   " << posPidSumMCTruth << " " << negPidSumMCTruth << endl;
  //cout <<"---------------------------------" << endl;
  
  Double_t ptContainer[dim+1];
  Double_t ptContainerMC[dim+1];
  ptContainer[0] = fCentrality;
  ptContainerMC[0] = fCentrality;

  for(Int_t i = 1; i <= dim; i++){
    ptContainer[i] = PtCh[i-1];
    ptContainerMC[i] = ptChMC[i-1];
    //cout <<" "<< PtCh[i-1]<<endl;
    //  cout<< " MC=" << ptChMC[i-1];
  }

  //cout << endl;

  fPtBinNplusNminusCh->Fill(ptContainer);
  fPtBinNplusNminusChTruth->Fill(ptContainerMC);

  fEventCounter->Fill(7);
  PostData(1, fListOfHistos);
  return;
  
}

//---------------------------------------------------------------------------------------
Bool_t AliEbyEHigherMomentsTask::ProperVertex(AliAODEvent *fAOD) const{
  
  Bool_t IsVtx = kFALSE;
  
  const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
  
  if(vertex) {
    Double32_t fCov[6];
    vertex->GetCovarianceMatrix(fCov);
    if(vertex->GetNContributors() > 0) {
      if(fCov[5] != 0) {
	
	Double_t lvx = vertex->GetX();
	Double_t lvy = vertex->GetY();
	Double_t lvz = vertex->GetZ();
	
	fEventCounter->Fill(5);
      
	if(TMath::Abs(lvx) < fVxMax) {
	  if(TMath::Abs(lvy) < fVyMax) {
	    if(TMath::Abs(lvz) < fVzMax) {
	      
	      fEventCounter->Fill(6);
	      fHistQA[0]->Fill(lvz);
	      fHistVxVy->Fill(lvx,lvy);
	      IsVtx = kTRUE; 
	      
	    }//Z-Vertex cut---
 	  }//Y-vertex cut--
	}//X-vertex cut---
      }//Covariance------
    }//Contributors check---
  }//If vertex-----------

  AliCentrality *centrality = fAOD->GetCentrality();
  if (centrality->GetQuality() != 0) IsVtx = kFALSE;
  
  return IsVtx;
}
//---------------------------------------------------------------------------------------

//---------------------------------------------------------------------------------------
Bool_t AliEbyEHigherMomentsTask::AcceptTrack(AliAODTrack* track) const{
  
  if(!track) return kFALSE;
  if( track->Charge() == 0 ) return kFALSE;

  Double_t pt = track->Pt();
  Double_t eta = track->Eta();
  if(!track->TestFilterBit(fAODtrackCutBit)) return kFALSE;
  if( pt < fPtLowerLimit || pt > fPtHigherLimit ) return kFALSE;
  if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) return kFALSE;
  
  fHistQA[1]->Fill(pt);
  fHistQA[2]->Fill(eta);
  fHistQA[3]->Fill(track->Phi());
 
  return kTRUE;
}
//------------------------------------------------------------------------

//------------------------------------------------------------------------
Int_t AliEbyEHigherMomentsTask::GetPtBin(Double_t pt){
  
  Int_t bin = 0;

  Double_t BinSize = (fPtHigherLimit - fPtLowerLimit)/fNptBins;
  
  for(Int_t iBin = 0; iBin < fNptBins; iBin++){
    
    Double_t xLow = fPtLowerLimit + iBin*BinSize;
    Double_t xHigh = fPtLowerLimit + (iBin+1)*BinSize;
    
    if( pt >= xLow && pt < xHigh){
      bin = iBin;
      //cout << "Pt Bin #" << bin <<"pt=" << pt << endl;
      break;
    }
    
  }//for 
  
  
  return bin;
}
//------------------------------------------------------------------------
//------------------------------------------------------------------------
void AliEbyEHigherMomentsTask::Terminate( Option_t * ){
  
  Info("AliEbyEHigherMomentTask"," Task Successfully finished");
  
}

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