--- /dev/null
+/**************************************************************************
+ * 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 //
+// //
+//=========================================================================//
+
+#include "TChain.h"
+#include "TList.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1D.h"
+#include "TH2F.h"
+#include "THnSparse.h"
+#include "TCanvas.h"
+#include "TParticle.h"
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+
+#include "AliESDVertex.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliAODInputHandler.h"
+#include "AliESD.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliStack.h"
+#include "AliESDtrackCuts.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliMCParticle.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 "AliEbyEHigherMomentsEffContTask.h"
+
+using std::cout;
+using std::endl;
+
+
+ClassImp(AliEbyEHigherMomentsEffContTask)
+
+
+//-----------------------------------------------------------------------
+AliEbyEHigherMomentsEffContTask::AliEbyEHigherMomentsEffContTask( const char *name )
+: AliAnalysisTaskSE( name ),
+ fListOfHistosQA(0),
+ fListOfHistos(0),
+ fAOD(0),
+ fArrayMC(0),
+ fPIDResponse(0),
+ fParticleSpecies(AliPID::kProton),
+ fAnalysisType("AOD"),
+ fCentralityEstimator("V0M"),
+ fCentrality(0),
+ fVxMax(3.),
+ fVyMax(3.),
+ fVzMax(10.),
+ fDCAxy(2.4),
+ fDCAz(3.),
+ fPtLowerLimit(0.3),
+ fPtHigherLimit(1.5),
+ fEtaLowerLimit(-0.8),
+ fEtaHigherLimit(0.8),
+ fRapidityCut(0.5),
+ fNSigmaCut(3.),
+ fTPCNClus(80),
+ fChi2perNDF(4.),
+ fAODtrackCutBit(1024),
+ fLabel(NULL),
+ fUsePid(kFALSE),
+ fCheckCont(kFALSE),
+ fEff(kFALSE),
+ fEventCounter(0),
+ fHistDCA(0),
+ fTPCSig(0),
+ fTPCSigA(0),
+ fTHnCentNplusNminusCh(0),
+ fTHnCentNplusNminus(0),
+ fTHnEff(0),
+ fTHnCont(0)
+{
+
+ for ( Int_t i = 0; i < 13; i++) {
+ fHistQA[i] = NULL;
+ }
+
+ for ( Int_t i = 0; i < 5; i++ ){
+ fTHnCentNplusNminusPid[i] = NULL;
+ }
+
+ DefineOutput(1, TList::Class()); // Outpput....
+ DefineOutput(2, TList::Class());
+
+}
+
+AliEbyEHigherMomentsEffContTask::~AliEbyEHigherMomentsEffContTask() {
+
+ if(fListOfHistosQA) delete fListOfHistosQA;
+ if(fListOfHistos) delete fListOfHistos;
+ if(fLabel[0] ) delete [] (fLabel[0]);
+ if(fLabel[1] ) delete [] (fLabel[1]);
+
+}
+
+//---------------------------------------------------------------------------------
+void AliEbyEHigherMomentsEffContTask::UserCreateOutputObjects() {
+ //For QA-Histograms
+ fListOfHistosQA = new TList();
+ fListOfHistosQA->SetOwner(kTRUE);
+
+ 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-90% centrality");
+ fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
+ fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
+ fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
+ fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
+ fListOfHistosQA->Add(fEventCounter);
+
+ fHistQA[0] = new TH1D("fHistQAvx", "Histo Vx After Cut", 400, -4., 4.);
+ fHistQA[1] = new TH1D("fHistQAvy", "Histo Vy After Cut", 400, -4., 4.);
+ fHistQA[2] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
+ fHistQA[3] = new TH1D("fHistQAvxA", "Histo Vx All", 500, -5., 5.);
+ fHistQA[4] = new TH1D("fHistQAvyA", "Histo Vy All", 500, -5., 5.);
+ fHistQA[5] = new TH1D("fHistQAvzA", "Histo Vz All", 500, -25., 25.);
+ fHistQA[6] = new TH1D("fHistQADcaXyC", "Histo DCAxy after cut", 500, -5., 5.);
+ fHistQA[7] = new TH1D("fHistQADcaZC", "Histo DCAz after cut", 500, -5., 5.);
+ fHistQA[8] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
+ fHistQA[9] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
+ fHistQA[10] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
+ fHistQA[11] = new TH1D("fHistQANCls","Number of TPC cluster",200,0,200);
+ fHistQA[12] = new TH1D("fHistQAChi2","Chi2 per NDF",100,0,10);
+
+ for(Int_t i = 0; i < 13; i++)
+ {
+ fListOfHistosQA->Add(fHistQA[i]);
+ }
+
+ fHistDCA = new TH2D("fHistDCA","DCAxy Vs DCAz", 500, -5., 5., 500, -5., 5.);
+ fTPCSig = new TH2D("fTPCSig","TPC signal",200, 0.0, 10. ,1000,0.,1000);
+ fTPCSig->SetMarkerColor(kRed);
+ fTPCSigA = new TH2D("fTPCSigA","TPC signal all ",200, 0.0, 10. ,1000,0.,1000);
+ fListOfHistosQA->Add(fHistDCA);
+ fListOfHistosQA->Add(fTPCSig);
+ fListOfHistosQA->Add(fTPCSigA);
+
+ const Int_t nDim = 3;
+ const Int_t nPid = 5;
+
+ TString Species[nPid] = {"Electron","Muon","Pion","Kaon","Proton"};
+
+ Int_t fBins[nPid][nDim];
+ Double_t fMin[nPid][nDim];
+ Double_t fMax[nPid][nDim];
+
+ for( Int_t i = 0; i < nPid; i++ ){
+ for( Int_t j = 0; j < nDim; j++ ){
+ fBins[i][j] = 0;
+ fMin[i][j] = 0.;
+ fMax[i][j] = 0.;
+ }
+ }
+
+
+ // if(fAnalysisType == "AOD" || fAnalysisType == "MCAOD")
+
+ //{
+ TString hname1;
+ TString htitle1, axisTitle1, axisTitle2;
+
+
+
+ if( fUsePid ){
+
+ Int_t gPid = (Int_t)fParticleSpecies;
+
+ if( gPid > 1 && gPid < 5 ){
+ //Pion----
+ fBins[2][0] = 100; fBins[2][1] = 1000; fBins[2][2] = 1000;
+ fMin[2][0] = -0.5; fMin[2][1] = -0.5; fMin[2][2] = -0.5;
+ fMax[2][0] = 99.5; fMax[2][1] = 999.5; fMax[2][2] = 999.5;
+ //Kaon------
+ fBins[3][0] = 100; fBins[3][1] = 600; fBins[3][2] = 600;
+ fMin[3][0] = -0.5; fMin[3][1] = -0.5; fMin[3][2] = -0.5;
+ fMax[3][0] = 99.5; fMax[3][1] = 599.5; fMax[3][2] = 599.5;
+ //Proton-----
+ fBins[4][0] = 100; fBins[4][1] = 400; fBins[4][2] = 400;
+ fMin[4][0] = -0.5; fMin[4][1] = -0.5; fMin[4][2] = -0.5;
+ fMax[4][0] = 99.5; fMax[4][1] = 399.5; fMax[4][2] = 399.5;
+
+ hname1 = "fCentNplusNminusPid"; hname1 +=gPid;
+ htitle1 = Species[gPid]; htitle1 +=" And Neg-"; htitle1 +=Species[gPid];
+ axisTitle1 = Species[gPid];
+ axisTitle2 ="Neg-"; axisTitle2 += Species[gPid];
+
+ fTHnCentNplusNminusPid[gPid] = new THnSparseD(hname1.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
+
+ fTHnCentNplusNminusPid[gPid]->GetAxis(0)->SetTitle("Centrality");
+ fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
+ fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle2.Data());
+
+ fListOfHistos->Add(fTHnCentNplusNminusPid[gPid]);
+
+ }//Pion, Koan and Proton-------
+
+ else{
+
+ Int_t fBinsX[nDim] = {100, 1500, 1500};
+ Double_t fMinX[nDim] = { -0.5, -0.5, -0.5 };
+ Double_t fMaxX[nDim] = { 99.5, 1499.5, 1499.5};
+
+ fTHnCentNplusNminus = new THnSparseD("fTHnCentNplusNminus","Cent-NplusChrg-NminusChrg", nDim, fBinsX, fMinX, fMaxX);
+ fTHnCentNplusNminus->GetAxis(0)->SetTitle("Centrality");
+ fTHnCentNplusNminus->GetAxis(1)->SetTitle("UnwantedPlus");
+ fTHnCentNplusNminus->GetAxis(2)->SetTitle("UnwantedMinus");
+ fListOfHistos->Add(fTHnCentNplusNminus);
+
+ }//Unwanted particle -------
+
+ }//fUsePid-------
+
+ else{
+
+ 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};
+
+ 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);
+ }//All charge hadrons---------
+
+
+
+ if( fAnalysisType == "MCAOD" ){
+
+ Double_t dCent[2] = {-0.5, 9.5};
+ Int_t iCent = 10;
+
+ Double_t dEta[2] = {-0.9, 0.9};
+ Int_t iEta = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ;
+
+ Double_t dRap[2] = {-0.5, 0.5};
+ Int_t iRap = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ;
+
+ Double_t dPhi[2] = {0.0, TMath::TwoPi()};
+ Int_t iPhi = 42;
+
+ Double_t dPt[2] = {0.1, 3.0};
+ Int_t iPt = Int_t((dPt[1]-dPt[0]) / 0.075);
+
+ Double_t dSign[2] = {-1.5, 1.5};
+ Int_t iSign = 3;
+
+ if(fEff){
+
+ Int_t nBinEff[15] = { iCent, iEta, iRap, iPhi, iPt, iSign, 2, 2, 2, iEta, iPhi, iPt, iEta, 2*iPhi+1, 2*iPt+1};
+ 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] };
+ 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] };
+
+ fTHnEff = new THnSparseF("fHnEff", "cent:etaMC:yMC:phiMC:ptMC:sign:findable:recStatus:pidStatus:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt", 15, nBinEff, nMinEff, nMaxEff);
+
+ fTHnEff->Sumw2();
+
+ fTHnEff->GetAxis(0)->SetTitle("centrality");
+ fTHnEff->GetAxis(1)->SetTitle("#eta_{MC}");
+ fTHnEff->GetAxis(2)->SetTitle("#it{y}_{MC}");
+ fTHnEff->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");
+ fTHnEff->GetAxis(4)->SetTitle("p_{T,MC} (GeV/#it{c})");
+ fTHnEff->GetAxis(5)->SetTitle("sign");
+ fTHnEff->GetAxis(6)->SetTitle("findable");
+ fTHnEff->GetAxis(7)->SetTitle("recStatus");
+ fTHnEff->GetAxis(8)->SetTitle("recPid");
+ fTHnEff->GetAxis(9)->SetTitle("#eta_{Rec}");
+ fTHnEff->GetAxis(10)->SetTitle("#varphi_{Rec} (rad)");
+ fTHnEff->GetAxis(11)->SetTitle("p_{T,Rec} (GeV/#it{c})");
+ fTHnEff->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}");
+ fTHnEff->GetAxis(13)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");
+ fTHnEff->GetAxis(14)->SetTitle("p_{T,MC}-p_{T,Rec} (GeV/#it{c})");
+
+ fListOfHistos->Add(fTHnEff);
+ }
+
+ if(fCheckCont){
+
+ Int_t nBinCont[14] = { iCent, iEta, iRap, iPhi, iPt, iSign, 8, iSign, iEta, iPhi, iPt, iEta, 2*iPhi+1, 2*iPt+1 };
+ 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]};
+ 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] };
+
+ fTHnCont = new THnSparseF("fHnCont", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",14, nBinCont, nMinCont, nMaxCont);
+ fTHnCont->Sumw2();
+
+ fTHnCont->GetAxis(0)->SetTitle("centrality");
+ fTHnCont->GetAxis(1)->SetTitle("#eta_{MC}");
+ fTHnCont->GetAxis(2)->SetTitle("#it{y}_{MC}");
+ fTHnCont->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");
+ fTHnCont->GetAxis(4)->SetTitle("p_{T,MC} (GeV/#it{c})");
+ fTHnCont->GetAxis(5)->SetTitle("sign");
+ fTHnCont->GetAxis(6)->SetTitle("contPart");
+ fTHnCont->GetAxis(7)->SetTitle("contSign");
+ fTHnCont->GetAxis(8)->SetTitle("#eta_{Rec}");
+ fTHnCont->GetAxis(9)->SetTitle("#varphi_{Rec} (rad)");
+ fTHnCont->GetAxis(10)->SetTitle("p_{T,Rec} (GeV/#it{c})");
+ fTHnCont->GetAxis(11)->SetTitle("#eta_{MC}-#eta_{Rec}");
+ fTHnCont->GetAxis(12)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");
+ fTHnCont->GetAxis(13)->SetTitle("p_{T,MC}-p_{T,Rec} (GeV/#it{c})");
+
+ fListOfHistos->Add(fTHnCont);
+ }
+
+ }//MCAOD------
+
+ PostData(1, fListOfHistosQA);
+ PostData(2, fListOfHistos);
+
+}
+
+
+//----------------------------------------------------------------------------------
+void AliEbyEHigherMomentsEffContTask::UserExec( Option_t * ){
+
+ fEventCounter->Fill(1);
+
+
+ if(fAnalysisType == "AOD") {
+
+ doAODEvent();
+
+ }//AOD--analysis-----
+
+ else if(fAnalysisType == "MCAOD") {
+
+ doMCAODEvent();
+
+ }
+
+ else return;
+
+
+
+
+ fEventCounter->Fill(8);
+
+ PostData(1, fListOfHistosQA);
+ PostData(2, fListOfHistos);
+
+}
+
+//--------------------------------------------------------------------------------------
+void AliEbyEHigherMomentsEffContTask::doAODEvent(){
+
+ //-------------------
+ Double_t nPlusCharge = 0.;
+ Double_t nMinusCharge = 0.;
+ Double_t nPartile = 0.;
+ Double_t nAntiParticle = 0.;
+ Int_t gPid = 0.;
+
+
+ //connect to the analysis mannager-----
+
+ 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;
+ }
+
+ fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+ if (!fAOD)
+ {
+ cout<< "ERROR 01: AOD not found " <<endl;
+ return;
+ }
+
+ fPIDResponse =inputHandler->GetPIDResponse();
+
+
+ if (!fPIDResponse){
+ AliFatal("This Task needs the PID response attached to the inputHandler");
+ return;
+ }
+
+
+ AliAODHeader *aodHeader = fAOD->GetHeader();
+
+ Int_t cent = -1;
+ cent = aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data());
+
+ if (cent == 0)
+ fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data());
+ else if (cent == 10 || cent == -1.)
+ fCentrality = -1;
+ else if (cent > 0 && cent < 9)
+ fCentrality = cent + 1;
+
+ if(fCentrality < 0 || fCentrality >= 10) return;
+
+ fEventCounter->Fill(2);
+
+ if(!ProperVertex()) return;
+
+ Int_t nTracks = fAOD->GetNumberOfTracks();
+
+ TExMap *trackMap = new TExMap();//Mapping matrix----
+
+ 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;
+ }
+
+ Double_t tpcSignalAll = aodTrack->GetTPCsignal();
+ fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
+
+ Int_t gID = aodTrack->GetID();
+
+ if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
+ }//1st track loop----
+
+ AliAODTrack* newAodTrack;
+
+ for( Int_t j = 0; j < nTracks; j++ )
+ {
+
+ AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
+
+ if(!aodTrack1) {
+ AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
+ continue;
+ }
+
+
+ if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
+
+ Int_t gID = aodTrack1->GetID();
+
+ //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
+ newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
+
+ Float_t dxy = 0., dz = 0.;
+
+ dxy = aodTrack1->DCA();
+ dz = aodTrack1->ZAtDCA();
+
+ Double_t pt = aodTrack1->Pt();
+ Double_t eta = aodTrack1->Eta();
+ Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
+ Double_t chi2ndf = aodTrack1->Chi2perNDF();
+
+ /*
+ if( fabs(dxy) > fDCAxy ) continue;
+ if( fabs(dz) > fDCAz ) continue;
+ //Extra cut on DCA---( Similar to BF Task.. )
+ if( fDCAxy !=0 && fDCAz !=0 ){
+ if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
+ }
+ */
+ if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
+ if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
+ if( nclus < fTPCNClus ) continue;
+ if( chi2ndf > fChi2perNDF ) continue;
+
+
+ fHistQA[6]->Fill(dxy);
+ fHistQA[7]->Fill(dz);
+ fHistQA[8]->Fill(pt);
+ fHistQA[9]->Fill(eta);
+ fHistQA[10]->Fill(aodTrack1->Phi());
+ fHistQA[11]->Fill(nclus);
+ fHistQA[12]->Fill(chi2ndf);
+ fHistDCA->Fill(dxy,dz);
+
+ Short_t gCharge = aodTrack1->Charge();
+
+ if(gCharge > 0) nPlusCharge += 1.;
+ if(gCharge < 0) nMinusCharge += 1.;
+
+ if( fUsePid ) {
+
+ gPid = (Int_t)fParticleSpecies;
+
+ Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
+ Double_t tpcSignal = newAodTrack->GetTPCsignal();
+ //Double_t rap = aodTrack1->Y(AliPID::ParticleMass(fParticleSpecies));
+ //Double_t tpcSignal = aodTrack1->GetTPCsignal();
+
+ if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
+
+ fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
+
+ Float_t nsigmaTPCPID = -999.;
+ Float_t nsigmaTOFPID = -999.;
+ //Float_t nsigmaTPCTOFPID = -999.;
+
+ nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
+ nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
+
+ if ( nsigmaTPCPID < fNSigmaCut ){
+
+ if (gCharge > 0) nPartile +=1.;
+ if( gCharge < 0 ) nAntiParticle +=1.;
+
+ }
+ }//fUsepid-----
+
+ }//--------- Track Loop to select with filterbit
+
+
+
+ Double_t fContainerCh[3] = { fCentrality, nPlusCharge, nMinusCharge};
+ Double_t fContainerPid[3] = { fCentrality, nPartile, nAntiParticle};
+
+ if( fUsePid ){
+
+ fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
+
+ // cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl;
+ }
+ else{
+
+ fTHnCentNplusNminusCh->Fill(fContainerCh);
+ }
+
+ //cout << "nCentrality "<< fCentrality <<", nPlusCharge="<< nPlusCharge << ", nMinusCharge=" << nMinusCharge << endl;
+
+ fEventCounter->Fill(7);
+
+ return;
+
+}
+//--------------------------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------------------------
+void AliEbyEHigherMomentsEffContTask::doMCAODEvent(){
+
+
+ //---------
+ Double_t nPlusCharge = 0.;
+ Double_t nMinusCharge = 0.;
+ Double_t nPartile = 0.;
+ Double_t nAntiParticle = 0.;
+ Int_t gPid = 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------
+ fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+ if (!fAOD)
+ {
+ cout<< "ERROR 01: AOD not found " <<endl;
+ return;
+ }
+
+ fPIDResponse =inputHandler->GetPIDResponse();
+
+
+ if (!fPIDResponse){
+ AliFatal("This Task needs the PID response attached to the inputHandler");
+ return;
+ }
+ // -- Get MC header
+ // ------------------------------------------------------------------
+
+ fArrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+ if (!fArrayMC) {
+ AliFatal("Error: MC particles branch not found!\n");
+ }
+
+ AliAODMCHeader *mcHdr=NULL;
+ mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+ if(!mcHdr) {
+ printf("MC header branch not found!\n");
+ return;
+ }
+
+ AliAODHeader *aodHeader = fAOD->GetHeader();
+
+ Int_t cent = -1;
+ cent = aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data());
+
+ if (cent == 0)
+ fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data());
+ else if (cent == 10 || cent == -1.)
+ fCentrality = -1;
+ else if (cent > 0 && cent < 9)
+ fCentrality = cent + 1;
+
+ if( fCentrality < 0 || fCentrality >= 10) return;
+
+ fEventCounter->Fill(2);
+
+
+
+ if(!ProperVertex()) return;
+
+ Int_t nTracks = fAOD->GetNumberOfTracks();
+
+ //Creat Lables ---
+ fLabel = new Int_t*[2];
+ fLabel[0] = new Int_t[nTracks]; //All charged hadrons----
+ fLabel[1] = new Int_t[nTracks]; //For Pid-----------
+ //Initialize labels----
+
+ if(!fLabel[0]){
+ AliError("Can't Get fLabel[0] ");
+ return;
+ }
+ if(!fLabel[1]){
+ AliError("Can't Get fLabel[1] ");
+ return;
+ }
+
+ for(Int_t i=0; i < 2; i++){
+ for(Int_t j=0; j < nTracks; j++){
+ fLabel[i][j] = 0;
+ }
+ }
+
+
+ TExMap *trackMap = new TExMap();//Mapping matrix----
+
+ 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;
+ }
+
+ Double_t tpcSignalAll = aodTrack->GetTPCsignal();
+ fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
+
+ Int_t gID = aodTrack->GetID();
+
+ if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
+
+ }//1st track loop----
+
+ AliAODTrack* newAodTrack;
+
+ for( Int_t j = 0; j < nTracks; j++ ){
+
+ AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
+
+ if(!aodTrack1) {
+ AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
+ continue;
+ }
+
+
+ if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
+
+
+ Int_t gID = aodTrack1->GetID();
+
+ //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
+
+ newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
+
+ Float_t dxy = 0., dz = 0.;
+
+ dxy = aodTrack1->DCA();
+ dz = aodTrack1->ZAtDCA();
+
+ //cout << dxy << endl;
+ Double_t pt = aodTrack1->Pt();
+ Double_t eta = aodTrack1->Eta();
+ Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
+ Double_t chi2ndf = aodTrack1->Chi2perNDF();
+
+ /*
+ if( fabs(dxy) > fDCAxy ) continue;
+ if( fabs(dz) > fDCAz ) continue;
+ //Extra cut on DCA---( Similar to BF Task.. )
+ if( fDCAxy !=0 && fDCAz !=0 ){
+ if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
+ }
+ */
+
+
+ if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
+ if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
+ // if( nclus < fTPCNClus ) continue;
+ //if( chi2ndf > fChi2perNDF ) continue;
+
+
+ fHistQA[6]->Fill(dxy);
+ fHistQA[7]->Fill(dz);
+ fHistQA[8]->Fill(pt);
+ fHistQA[9]->Fill(eta);
+ fHistQA[10]->Fill(aodTrack1->Phi());
+ fHistQA[11]->Fill(nclus);
+ fHistQA[12]->Fill(chi2ndf);
+ fHistDCA->Fill(dxy,dz);
+
+ Short_t gCharge = aodTrack1->Charge();
+
+ if( gCharge == 0 ) continue;
+
+ Int_t label = TMath::Abs(aodTrack1->GetLabel());
+ //fill the labels--------
+ fLabel[0][j] = label;
+
+ if(gCharge > 0) nPlusCharge += 1.;
+ if(gCharge < 0) nMinusCharge += 1.;
+
+ if( fUsePid ) {
+
+ gPid = (Int_t)fParticleSpecies;;
+ Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
+ Double_t tpcSignal = newAodTrack->GetTPCsignal();
+
+ if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
+
+ fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
+
+ Float_t nsigmaTPCPID = -999.;
+ Float_t nsigmaTOFPID = -999.;
+ //Float_t nsigmaTPCTOFPID = -999.;
+
+ nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
+ nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
+
+ if ( nsigmaTPCPID < fNSigmaCut ){
+
+ //Fill the labels----
+ fLabel[1][j] = label;
+
+ if (gCharge > 0) nPartile +=1.;
+ if( gCharge < 0 ) nAntiParticle +=1.;
+
+ }
+ }//fUsepid-----
+ //Check Contamination-------------------
+
+ if( fCheckCont ){
+ CheckContTrackAOD(label, gCharge, j);
+ }
+
+ }//--------- Track Loop to select with filterbit
+
+ if(fEff) FillEffSparse();
+
+ //cout << "Cent=" << fCentrality << " MC-PlusChrg=" << nPlusCharge << " MC-MinusChrg=" << nMinusCharge << endl;
+
+ //cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl;
+
+ Double_t fContainerCh[3] = { fCentrality, nPlusCharge, nMinusCharge};
+ Double_t fContainerPid[3] = { fCentrality, nPartile, nAntiParticle};
+
+ if( fUsePid ){
+
+ fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
+
+ }
+ else{
+
+ fTHnCentNplusNminusCh->Fill(fContainerCh);
+ }
+
+ fEventCounter->Fill(7);
+
+ return;
+
+}
+//---------------------------------------------------------------------------------------
+
+//---------------------------------------------------------------------------------------
+Bool_t AliEbyEHigherMomentsEffContTask::ProperVertex(){
+
+ 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);
+
+ fHistQA[3]->Fill(lvx);fHistQA[4]->Fill(lvy);fHistQA[5]->Fill(lvz);
+
+ if(TMath::Abs(lvx) < fVxMax) {
+ if(TMath::Abs(lvy) < fVyMax) {
+ if(TMath::Abs(lvz) < fVzMax) {
+
+ fEventCounter->Fill(6);
+ fHistQA[0]->Fill(lvx);fHistQA[1]->Fill(lvy);fHistQA[2]->Fill(lvz);
+
+ IsVtx = kTRUE;
+
+ }//Z-Vertex cut---
+ }//Y-vertex cut--
+ }//X-vertex cut---
+ }//Covariance------
+ }//Contributors check---
+ }//If vertex-----------
+
+ return IsVtx;
+}
+//---------------------------------------------------------------------------------------
+
+//---------------------------------------------------------------------------------------
+void AliEbyEHigherMomentsEffContTask::FillEffSparse(){
+
+ //For Efficiency-------------------------------------
+
+ Int_t nTracks = fAOD->GetNumberOfTracks();
+ Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies);
+
+ 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;
+ }
+
+ /* TDatabasePDG *db = TDatabasePDG::Instance();
+ TParticlePDG *part = 0x0;
+
+ if (!(part = db->GetParticle(partMC->PdgCode()))) continue;
+ if ((part = db->GetParticle(partMC->PdgCode()))->Charge() == 0.0) continue;
+ */
+ if(!(partMC->PdgCode())) continue;
+ if(partMC->Charge() == 0) continue;
+ //-(1) Check for primary----
+ if(!partMC->IsPhysicalPrimary()) continue;
+
+ //-(2) Basic track cuts--------
+
+ if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
+ if (partMC->Pt() < fPtLowerLimit || partMC->Pt() > fPtHigherLimit) continue;
+
+ Double_t rap = partMC->Y();
+
+ if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
+
+ if(TMath::Abs(partMC->GetPdgCode()) != gPdgCode) continue;
+
+ Short_t gCharge = partMC->Charge();
+
+ Float_t sign = (gCharge < 0 ) ? -1. : 1.;
+
+ Float_t findable = 1.;// Float_t(fHelper->IsParticleFindable(idxMC)); XXX
+
+ // (3) Initialize rec. variables-----
+ Float_t recStatus = 0.;
+ Float_t recPid = 0.;
+ // (4) Initialize track parameters------
+ Float_t etaRec = 0.;
+ Float_t phiRec = 0.;
+ Float_t ptRec = 0.;
+
+
+
+ // (5) Loop over all Labels--------
+
+
+ for( Int_t iRec = 0; iRec < nTracks; iRec++ ){
+
+ if( iMC != fLabel[0][iRec] ) continue;
+ else{
+ recStatus = 1.;
+ if( iMC == fLabel[1][iRec] )
+ recPid = 1.;
+
+ AliAODTrack* aodTrack = NULL;
+
+ if(fAOD)
+
+ aodTrack = fAOD->GetTrack(iRec);
+
+ if(aodTrack){
+
+ etaRec = aodTrack->Eta();
+ phiRec = aodTrack->Phi();
+ ptRec = aodTrack->Pt();
+
+ }//aodTrack
+
+ break;
+ }
+
+
+
+ }//all charged track rec. check---
+
+
+ //Fill ThnSparse----
+
+ Double_t EffContainer[15] = {fCentrality, partMC->Eta(), partMC->Y(), partMC->Phi(), partMC->Pt(), sign,findable, recStatus, recPid, etaRec, phiRec, ptRec, partMC->Eta()-etaRec, partMC->Phi()-phiRec, partMC->Pt()-ptRec};
+
+ fTHnEff->Fill(EffContainer);
+
+ }//iMC Track loop-----
+
+ return;
+}
+//---------------------------------------------------------------------------------------
+
+//---------------------------------------------------------------------------------------
+void AliEbyEHigherMomentsEffContTask::CheckContTrackAOD(Int_t label, Float_t sign, Int_t idxTrack){
+
+ AliAODMCParticle* particle = (AliAODMCParticle*)fArrayMC->At(label);
+
+ if (!particle)
+ return;
+
+ Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies);
+ Bool_t isSecondaryFromWeakDecay = kFALSE;
+ Bool_t isSecondaryFromMaterial = kFALSE;
+
+ // -- Check if correctly identified
+ if (particle->GetPdgCode() == (sign*gPdgCode)) {
+
+ // Check if is physical primary -> all ok
+ if(particle->IsPhysicalPrimary())
+ return;
+
+ // -- Check if secondaries from material or weak decay
+ isSecondaryFromWeakDecay = particle->IsSecondaryFromWeakDecay();
+ isSecondaryFromMaterial = particle->IsSecondaryFromMaterial();
+
+ }
+
+ // -- Get MC pdg
+ Float_t contSign = 0.;
+
+ // TDatabasePDG *db = TDatabasePDG::Instance();
+ //TParticlePDG *part = 0x0;
+ /*
+ if(( part = db->GetParticle(particle->PdgCode()))->Charge() == 0.) contSign = 0.;
+ else if((part = db->GetParticle(particle->PdgCode()))->Charge() < 0.) contSign = -1.;
+ else if((part = db->GetParticle(particle->PdgCode()))->Charge() > 0.) contSign = 1.;
+ */
+
+ if(particle->Charge() == 0.) contSign = 0.;
+ if(particle->Charge() < 0. ) contSign = -1;
+ if(particle->Charge() > 0. ) contSign = 1.;
+
+ // -- Get contaminating particle
+ Float_t contPart = 0;
+ if (isSecondaryFromWeakDecay) contPart = 7; // probeParticle from WeakDecay
+ else if (isSecondaryFromMaterial) contPart = 8; // probeParticle from Material
+ else {
+ if (TMath::Abs(particle->GetPdgCode()) == 211) contPart = 1; // pion
+ else if (TMath::Abs(particle->GetPdgCode()) == 321) contPart = 2; // kaon
+ else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton
+ else if (TMath::Abs(particle->GetPdgCode()) == 11) contPart = 4; // electron
+ else if (TMath::Abs(particle->GetPdgCode()) == 13) contPart = 5; // muon
+ else contPart = 6; // other
+ }
+
+ // -- Get Reconstructed values
+ Float_t etaRec = 0.;
+ Float_t phiRec = 0.;
+ Float_t ptRec = 0.;
+
+ AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(idxTrack));
+
+ if (aodTrack) {
+ // if no track present (which should not happen)
+ // -> pt = 0. , which is not in the looked at range
+
+ // -- Get Reconstructed eta/phi/pt
+ etaRec = aodTrack->Eta();
+ phiRec = aodTrack->Phi();
+ ptRec = aodTrack->Pt();
+ }
+
+ // -- Fill THnSparse
+ Double_t ContContainer[14] = { fCentrality, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, contPart, contSign, etaRec, phiRec, ptRec, particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec };
+
+ fTHnCont->Fill(ContContainer);
+
+ return;
+
+}
+//------------------------------------------------------------------------
+void AliEbyEHigherMomentsEffContTask::Terminate( Option_t * ){
+
+ Info("AliEbyEHigherMomentsEffContTask"," Task Successfully finished");
+
+}
+
+//------------------------------------------------------------------------
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff