[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskPIDconfig.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * 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.                  *
 **************************************************************************/

/* $Id: AliAnalysisTaskPIDconfig.cxx 43811 2014-10-11 Naghmeh Mohammadi $ */

#include "TChain.h"
#include "TTree.h"
#include "TList.h"
#include "TMath.h"
#include "TObjArray.h"
#include "TCanvas.h"
#include "TGraphErrors.h"
#include "TString.h"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h" 
#include "TH2D.h"                  
#include "TH3D.h"
#include "TArrayF.h"
#include "TF1.h"
#include "TROOT.h"
#include "stdio.h"
#include "TCutG.h"


#include "AliTHn.h"
#include "AliLog.h"
#include "AliAnalysisManager.h"
#include "AliESDEvent.h"
#include "AliAODInputHandler.h"
#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliAODInputHandler.h"
#include "AliCollisionGeometry.h"
#include "AliGenEventHeader.h"
#include "AliAnalysisUtils.h"
#include "AliPIDCombined.h"        
#include "AliAnalysisTask.h"
#include "AliAODHandler.h"
#include <AliInputEventHandler.h>
#include <AliVEventHandler.h>
#include <AliVParticle.h>
#include <AliVTrack.h>
#include <AliTPCPIDResponse.h>
#include <AliTOFPIDResponse.h>
#include "AliAnalysisTaskPIDconfig.h"
#include "AliAnalysisTaskSE.h"
#include "AliAODPid.h"
#include "AliPhysicsSelection.h"
#include "AliCentralitySelectionTask.h"
#include "AliCentrality.h"
#include "AliKFParticle.h"
#include "AliKFVertex.h"
#include "AliPID.h"
#include "AliPIDResponse.h"
#include "AliCFContainer.h"
#include "AliCFManager.h"
#include "AliVEvent.h"
#include "AliAODVZERO.h"

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

ClassImp(AliAnalysisTaskPIDconfig)
//ClassImp()
//___________________________________________________________________
AliAnalysisTaskPIDconfig::AliAnalysisTaskPIDconfig():
AliAnalysisTaskSE(),
fVevent(0),
fESD(0),
fAOD(0),
fPIDResponse(0),
fTriggerSelection(0),
fCentralityPercentileMin(0.),
fCentralityPercentileMax(5.),
fFilterBit(128),
fDCAxyCut(-1),
fDCAzCut(-1),
fData2011(kFALSE),
fTriggerMB(kTRUE),
fTriggerCentral(kFALSE),
fUseCentrality(kTRUE),
fCutTPCmultiplicityOutliersAOD(kFALSE),
fPIDcuts(kFALSE),
fCentralityEstimator("V0M"),
fContourCutList(0),
fListQA(0x0),
fListQAtpctof(0x0),
fListQAInfo(0x0),
fhistCentralityPass(0),
fNoEvents(0),
fpVtxZ(0),
fhistDCABefore(0),
fhistDCAAfter(0),
fhistPhiDistBefore(0),
fhistPhiDistAfter(0),
fhistEtaDistBefore(0),
fhistEtaDistAfter(0),
fTPCvsGlobalMultBeforeOutliers(0),
fTPCvsGlobalMultAfterOutliers(0),
fTPCvsGlobalMultAfter(0),
fHistBetavsPTOFbeforePID(0),
fHistdEdxVsPTPCbeforePID(0),
fHistBetavsPTOFafterPID(0),
fHistdEdxVsPTPCafterPID(0),
fhistNsigmaP(0),
fhistNsigmaPt(0)
//fSparseSpecies(0),
//fvalueSpecies(0)
{
    //Dummy Constructor
}


//

AliAnalysisTaskPIDconfig::AliAnalysisTaskPIDconfig(const char *name):
AliAnalysisTaskSE(name),
fVevent(0),
fESD(0),
fAOD(0),
fPIDResponse(0),
fTriggerSelection(0),
fCentralityPercentileMin(0.),
fCentralityPercentileMax(5.),
fFilterBit(1),
fDCAxyCut(-1),
fDCAzCut(-1),
fData2011(kFALSE),
fTriggerMB(kTRUE),
fTriggerCentral(kFALSE),
fUseCentrality(kTRUE),
fCutTPCmultiplicityOutliersAOD(kFALSE),
fPIDcuts(kFALSE),
fCentralityEstimator("V0M"),
fContourCutList(0),
fListQA(0x0),
fListQAtpctof(0x0),
fListQAInfo(0x0),
fhistCentralityPass(0),
fNoEvents(0),
fpVtxZ(0),
fhistDCABefore(0),
fhistDCAAfter(0),
fhistPhiDistBefore(0),
fhistPhiDistAfter(0),
fhistEtaDistBefore(0),
fhistEtaDistAfter(0),
fTPCvsGlobalMultBeforeOutliers(0),
fTPCvsGlobalMultAfterOutliers(0),
fTPCvsGlobalMultAfter(0),
fHistBetavsPTOFbeforePID(0),
fHistdEdxVsPTPCbeforePID(0),
fHistBetavsPTOFafterPID(0),
fHistdEdxVsPTPCafterPID(0),
fhistNsigmaP(0),
fhistNsigmaPt(0)
//fSparseSpecies(0),
//fvalueSpecies(0)
{
    //fvalueSpecies = new Double_t[9];
    //Default Constructor
    DefineInput(0,TChain::Class());
    DefineOutput(1,TList::Class());
}

//_____________________________________________________________________
AliAnalysisTaskPIDconfig::~AliAnalysisTaskPIDconfig()
{
    //Destructor
    
 //   delete fPID;
  //  delete fPIDqa;
 //   delete fTrackCuts;
   // delete fSparseSpecies;
    //delete []fvalueSpecies;


}
//______________________________________________________________________
void AliAnalysisTaskPIDconfig::UserCreateOutputObjects()
{
    //
    // Create the output QA objects
    //
    
    AliLog::SetClassDebugLevel("AliAnalysisTaskPIDconfig",10);
    
    //input hander
    AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
    AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
    if (!inputHandler) AliFatal("Input handler needed");
    
    //pid response object
    fPIDResponse=inputHandler->GetPIDResponse();
    if (!fPIDResponse) AliError("PIDResponse object was not created");
    
    //
    fListQA=new TList;
    fListQA->SetOwner();
    
    fListQAtpctof=new TList;
    fListQAtpctof->SetOwner();
    fListQAtpctof->SetName("PID_TPC_TOF");
    
    fListQAInfo=new TList;
    fListQAInfo->SetOwner();
    fListQAInfo->SetName("Event_Track_Info");
    
    fListQA->Add(fListQAtpctof);
    fListQA->Add(fListQAInfo);
    
    SetupTPCTOFqa();
    SetupEventInfo();
    
    PostData(1,fListQA);
}
//______________________________________________________________________
void AliAnalysisTaskPIDconfig::UserExec(Option_t*){
    //Main loop
    //Called for each event
    
    // create pointer to event
    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
    
    Int_t ntracks=fAOD->GetNumberOfTracks();


    if(!(fESD || fAOD)){
        printf("ERROR: fESD & fAOD not available\n");
        return;
    }
    fVevent = dynamic_cast<AliVEvent*>(InputEvent());
    if (!fVevent) {
        printf("ERROR: fVevent not available\n");
        return;
    }
    
    Bool_t pass = kFALSE;
    
    CheckCentrality(fVevent,pass);
    
    if(!pass){ return;}
    
    const AliVVertex *pVtx = fVevent->GetPrimaryVertex();
    
    Double_t pVtxZ = -999;
    pVtxZ = pVtx->GetZ();
    
    if(TMath::Abs(pVtxZ)>10) return;

    TH1F *hNoEvents = (TH1F*)fListQAInfo->At(1);
    TH1F *histpVtxZ = (TH1F*)fListQAInfo->At(2);

    if(hNoEvents) hNoEvents->Fill(0);
    if(histpVtxZ) histpVtxZ->Fill(pVtxZ);

    if(ntracks<2) return;
    
   // if(!pass) return;
    
    TH2F *HistTPCvsGlobalMultBeforeOutliers = (TH2F*)fListQAInfo->At(3);

    TH2F *HistTPCvsGlobalMultAfterOutliers = (TH2F*)fListQAInfo->At(4);

    
    Float_t multTPC(0.); // tpc mult estimate
    Float_t multGlobal(0.); // global multiplicity
    
    const Int_t nGoodTracks = fVevent->GetNumberOfTracks();
    if(!fData2011) { // cut on outliers
      
        for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill tpc mult
            AliAODTrack* AODtrack =dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
            if (!AODtrack) continue;
            if (!(AODtrack->TestFilterBit(1))) continue;
            if ((AODtrack->Pt() < .2) || (AODtrack->Pt() > 5.0) || (TMath::Abs(AODtrack->Eta()) > .8) || (AODtrack->GetTPCNcls() < 70)  || (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) || (AODtrack->Chi2perNDF() < 0.2)) continue;
            multTPC++;
        }//track loop

        for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill global mult
            AliAODTrack *AODtrack=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
            if (!AODtrack) continue;
            if (!(AODtrack->TestFilterBit(16))) continue;
            if ((AODtrack->Pt() < .2) || (AODtrack->Pt() > 5.0) || (TMath::Abs(AODtrack->Eta()) > .8) || (AODtrack->GetTPCNcls() < 70) || (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) || (AODtrack->Chi2perNDF() < 0.1)) continue;
            Double_t b[2] = {-99., -99.};
            Double_t bCov[3] = {-99., -99., -99.};
            AliAODTrack copy(*AODtrack);
            if (!(copy.PropagateToDCA(fVevent->GetPrimaryVertex(), fVevent->GetMagneticField(), 100., b, bCov))) continue;
            if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
            multGlobal++;
        } //track loop
        
        HistTPCvsGlobalMultBeforeOutliers->Fill(multGlobal,multTPC);

        if(multTPC < (-40.3+1.22*multGlobal) || multTPC > (32.1+1.59*multGlobal)){ pass = kFALSE;}
        
        if(!pass) return;
        HistTPCvsGlobalMultAfterOutliers->Fill(multGlobal,multTPC);

    }
    
    
    if(fData2011) { // cut on outliers
        //Float_t multTPC(0.); // tpc mult estimate
        //Float_t multGlob(0.); // global multiplicity
        for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill tpc mult
            AliAODTrack *AODtrack=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
            if (!AODtrack) continue;
            if (!(AODtrack->TestFilterBit(1))) continue;
            if ((AODtrack->Pt() < .2) || (AODtrack->Pt() > 5.0) || (TMath::Abs(AODtrack->Eta()) > .8) || (AODtrack->GetTPCNcls() < 70)  || (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) || (AODtrack->Chi2perNDF() < 0.2)) continue;
            multTPC++;
        }
        for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill global mult
            AliAODTrack *AODtrack=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
            if (!AODtrack) continue;
            if (!(AODtrack->TestFilterBit(16))) continue;
            if ((AODtrack->Pt() < .2) || (AODtrack->Pt() > 5.0) || (TMath::Abs(AODtrack->Eta()) > .8) || (AODtrack->GetTPCNcls() < 70) || (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) || (AODtrack->Chi2perNDF() < 0.1)) continue;
            Double_t b[2] = {-99., -99.};
            Double_t bCov[3] = {-99., -99., -99.};
            AliAODTrack copy(*AODtrack);
            if (!(copy.PropagateToDCA(fVevent->GetPrimaryVertex(), fVevent->GetMagneticField(), 100., b, bCov))) continue;
            if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
            multGlobal++;
            
        } //track loop
        
        HistTPCvsGlobalMultBeforeOutliers->Fill(multGlobal,multTPC);

        if(multTPC < (-36.73 + 1.48*multGlobal) || multTPC > (62.87 + 1.78*multGlobal)){pass = kFALSE;}
        
        if(!pass) return;
        HistTPCvsGlobalMultAfterOutliers->Fill(multGlobal,multTPC);

    }
 

     for(Int_t itrack = 0; itrack < ntracks; itrack++){

        AliAODTrack *track=dynamic_cast<AliAODTrack*>(fVevent->GetTrack(itrack));
        if(!track) continue;
        
        Float_t dcaXY = track->DCA();
        Float_t dcaZ  = track->ZAtDCA();
         
        TH2F* HistDCAbefore =(TH2F*)fListQAInfo->At(5);
        HistDCAbefore->Fill(dcaZ,dcaXY);
        
        Double_t p = -999, /*pTPC = -999,*/ pT = -999, phi = -999, eta = -999, dEdx =-999;
        Double_t length = -999., beta =-999, tofTime = -999., tof = -999.;
        Double_t c = TMath::C()*1.E-9;// m/ns

        //cout<<"track->GetFilterMap()= "<<track->GetFilterMap()<<endl;
        if(!track->TestFilterBit(fFilterBit)) continue;

        //Float_t dcaXY = -999, dcaZ = -999;
        p=track->P();
        /*pTPC=track->GetTPCmomentum();*/
        pT=track->Pt();
        phi=track->Phi();
        eta=track->Eta();
        dEdx=track->GetDetPid()->GetTPCsignal();
            
        if ( (track->IsOn(AliAODTrack::kTOFin)) &&
        (track->IsOn(AliAODTrack::kTIME)) &&  (track->IsOn(AliAODTrack::kTOFout))) {
      //  if ( (track->IsOn(AliAODTrack::kTOFin)) &&
      //              (track->IsOn(AliAODTrack::kTOFout))  ) {
                
                tofTime = track->GetTOFsignal();//in ps
                length = track->GetIntegratedLength();
                
                tof = tofTime*1E-3; // ns
                //cout<<"tof = "<<tof<<endl;
                if (tof <= 0)continue;
                //cout<<"length = "<<length<<endl;
                if (length <= 0) continue;
                
                length = length*0.01; // in meters
                tof = tof*c;
                beta = length/tof;
                
                TH2F *HistBetavsPTOFbeforePID = (TH2F*)fListQAInfo->At(6);
                HistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
            }//TOF signal

         TH2F *HistdEdxVsPTPCbeforePID = (TH2F*)fListQAInfo->At(7);
         HistdEdxVsPTPCbeforePID -> Fill(p*track->Charge(),dEdx); //TPC signal
         
         
            //QA plot
            TH1F *HistPhiDistBefore = (TH1F*)fListQAInfo->At(8);
            HistPhiDistBefore->Fill(phi);
            //
            TH1F *HistEtaDistBefore = (TH1F*)fListQAInfo->At(9);
            HistEtaDistBefore->Fill(eta);
         

            if(pT<0.1) continue;
            if(TMath::Abs(eta)>0.8) continue;
         
            Int_t TPCNcls = track->GetTPCNcls();
            
            if(TPCNcls<70 || dEdx<10) continue;
         
            // fill QA histograms

            TH2F* HistDCAAfter =(TH2F*)fListQAInfo->At(10);
            HistDCAAfter->Fill(dcaZ,dcaXY);
        
            TH1F *HistPhiDistAfter = (TH1F*)fListQAInfo->At(11);
            HistPhiDistAfter->Fill(phi);
                
            TH1F *HistEtaDistAfter = (TH1F*)fListQAInfo->At(12);
            HistEtaDistAfter->Fill(eta);
         
            Bool_t pWithinRange = kFALSE;
            Int_t pRange = -999;
            TCutG *cut[3][10];
            if(fPIDcuts){
                if(fContourCutList){cout<<"The contour file has been retrieved"<<endl;}
                
                TGraph *ContourCut[3][10];
                //Double_t plow[10] = {0.2,0.5,1,1.5,2,2.5,3,3.5,4,4.5};
                //Double_t phigh[10] = {0.5,1,1.5,2,2.5,3,3.5,4,4.5,5};
                Double_t phigh=0 ,plow=0;
                Double_t pBins[11] = {0.2,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5};
                TString Graph_Name = "contourlines_";
                TString species[3] = {"pion","kaon","proton"};
                for(int i=0;i<3;i++){
                    for(int j=0;j<10;j++){
                        
                        if(p>pBins[j] && p<pBins[j+1]){

                        //if(p>plow[j] && p<phigh[j]){
                            pWithinRange = kTRUE;
                            pRange = j;
                            
                            TList *Species_contours = (TList*)fContourCutList->Get(species[i]);
                            phigh = pBins[j+1]*10;
                            plow = pBins[j]*10;

                            Graph_Name += species[i];
                            Graph_Name += Form("%.f%.f-%i%icent",plow,phigh,fCentralityPercentileMin,fCentralityPercentileMax);

                            //TList *Pcontours = (TList*)Scontours->At(j);
                            //if (!Scontours) return;

                           // if (!Pcontours || !Scontours) return;
                            ContourCut[i][j] = (TGraph*)Species_contours->FindObject(Graph_Name);
                            cut[i][j] = new TCutG("cut",ContourCut[i][j]->GetN(),ContourCut[i][j]->GetX(),ContourCut[i][j]->GetY());
                            //ContourCut[i][j] = (TGraph *)fContourCutList->At(10*i+j);
                        }
                    }
                }
            }

            for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
                //TOF nSigma
                Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
                Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
                if(fPIDcuts && ispecie>1 && ispecie<5 && pWithinRange){// for pions, kaons and protons only
                    if(cut[ispecie-2][pRange]->IsInside(nSigmaTOF,nSigmaTPC)){
                        pass = kTRUE;
                    }
                    else{
                        pass = kFALSE;
                        continue;
                    }
                }
                
                //TPC and TOF cuts, TPC TOF nsigma vs. momentum
                if(pass){
                    TH3 *hist1 = (TH3*)fListQAtpctof->At(ispecie);
                    if (hist1){
                        hist1->Fill(nSigmaTPC,nSigmaTOF,p);}
                
                    TH3 *hist2 = (TH3*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);
                    if (hist2){
                        hist2->Fill(nSigmaTPC,nSigmaTOF,pT);}
                }
                
            }
         

     }//track loop
    
    TH2F *HistTPCvsGlobalMultAfter = (TH2F*) fListQAInfo->At(13);
    HistTPCvsGlobalMultAfter->Fill(multGlobal,multTPC);
    
}
//_________________________________________
void AliAnalysisTaskPIDconfig::CheckCentrality(AliVEvent* event, Bool_t &centralitypass)
{
    // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
    if (!fUseCentrality) AliFatal("No centrality method set! FATAL ERROR!");
    Double_t centvalue = event->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
     //cout << "Centrality evaluated-------------------------: " << centvalue <<endl;
    if ((centvalue >= fCentralityPercentileMin) && (centvalue < fCentralityPercentileMax))
    {
        TH1F *hCentralityPass = (TH1F*)fListQAInfo->At(0);
        hCentralityPass->Fill(centvalue);
        //cout << "--------------Fill pass-------------------------"<<endl;
        centralitypass = kTRUE;
    }
    
}
//______________________________________________________________________________
void AliAnalysisTaskPIDconfig::SetupTPCTOFqa()
{
    //
    // Create the qa objects for TPC + TOF combination
    

    //TPC and TOF signal vs. momentum
    for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
        fhistNsigmaP = new TH3F(Form("NsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),Form("TPC n#sigma vs. TOF n#sigma %s vs. p ;TPC n#sigma;TOF n#sigma;p [GeV]",AliPID::ParticleName(ispecie)),200,-20,20,200,-20,20,60,0.1,6);
        fListQAtpctof->Add(fhistNsigmaP);
    }
    //TPC and TOF signal vs. transverse momentum
    for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
        fhistNsigmaPt = new TH3F(Form("NsigmaPt_TPC_TOF_%s",AliPID::ParticleName(ispecie)),Form("TPC n#sigma vs. TOF n#sigma %s vs. Pt ;TPC n#sigma;TOF n#sigma;pT [GeV]",AliPID::ParticleName(ispecie)),200,-20,20,200,-20,20,60,0.1,6);
        fListQAtpctof->Add(fhistNsigmaPt);
    }

}
//______________________________________________________________________________
void AliAnalysisTaskPIDconfig::SetupEventInfo()
{
    //event and track info
    
    fhistCentralityPass = new TH1F("fcentralityPass","centralityPass", 100,0,100);
    fListQAInfo->Add(fhistCentralityPass);

    fNoEvents = new TH1F("number of events","no. of events",1,0,1);
    fListQAInfo->Add(fNoEvents);

    fpVtxZ = new TH1F("pVtxZ","pVtxZ",100,-20,20);
    fListQAInfo->Add(fpVtxZ);
    
    fTPCvsGlobalMultBeforeOutliers = new TH2F("TPC vs. Global Multiplicity Before","TPC vs. Global Multiplicity Before",500,0,6000,500,0,6000);
    fListQAInfo->Add(fTPCvsGlobalMultBeforeOutliers);
    
    fTPCvsGlobalMultAfterOutliers = new TH2F("TPC vs. Global Multiplicity After outliers","TPC vs. Global Multiplicity After outliers",500,0,6000,500,0,6000);
    fListQAInfo->Add(fTPCvsGlobalMultAfterOutliers);
    
    fhistDCABefore = new TH2F("DCA xy vs z (before)","DCA before",200,0,10,200,0,10);
    fListQAInfo->Add(fhistDCABefore);
    
    fHistBetavsPTOFbeforePID = new TH2F("momentum vs beta before PID","momentum vs beta before PID",1000,-10.,10.,1000,0,1.2);
    fListQAInfo->Add(fHistBetavsPTOFbeforePID);
    
    fHistdEdxVsPTPCbeforePID = new TH2F("momentum vs dEdx before PID","momentum vs dEdx before PID",1000,-10.,10.,1000,0,1000);
    fListQAInfo->Add(fHistdEdxVsPTPCbeforePID);
    
    fhistPhiDistBefore = new TH1F("Phi Distribution Before Cuts","Phi Distribution Before Cuts",200,0,6.4);
    fListQAInfo->Add(fhistPhiDistBefore);

    fhistEtaDistBefore = new TH1F("Eta Distribution Before Cuts","Eta Distribution Before Cuts",200,-10,10);
    fListQAInfo->Add(fhistEtaDistBefore);
    
    fhistDCAAfter = new TH2F("DCA xy vs z (after)","DCA after",200,0,10,200,0,10);
    fListQAInfo->Add(fhistDCAAfter);
    
    fhistPhiDistAfter = new TH1F("Phi Distribution After Cuts","Phi Distribution After Cuts",200,0,6.4);
    fListQAInfo->Add(fhistPhiDistAfter);
    
    fhistEtaDistAfter = new TH1F("Eta Distribution After Cuts","Eta Distribution After Cuts",200,-10,10);
    fListQAInfo->Add(fhistEtaDistAfter);
    
    fTPCvsGlobalMultAfter = new TH2F("TPC vs. Global Multiplicity After","TPC vs. Global Multiplicity After",500,0,6000,500,0,6000);
    fListQAInfo->Add(fTPCvsGlobalMultAfter);
    
//    fHistBetavsPTOFafterPID = new TH2F("momentum vs beta after PID","momentum vs beta after PID",1000,-10.,10.,1000,0,1.2);
//    fListQAInfo->Add(fHistBetavsPTOFafterPID);
    
//    fHistdEdxVsPTPCafterPID = new TH2F("momentum vs dEdx after PID","momentum vs dEdx after PID",1000,-10.,10.,1000,0,1000);
//    fListQAInfo->Add(fHistdEdxVsPTPCafterPID);

}
Commit	Line	Data
2c03c485	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	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	/* $Id: AliAnalysisTaskPIDconfig.cxx 43811 2014-10-11 Naghmeh Mohammadi $ */
	17
	18	#include "TChain.h"
	19	#include "TTree.h"
	20	#include "TList.h"
	21	#include "TMath.h"
	22	#include "TObjArray.h"
	23	#include "TCanvas.h"
	24	#include "TGraphErrors.h"
	25	#include "TString.h"
	26	#include "TFile.h"
	27	#include "TH1F.h"
	28	#include "TH2F.h"
	29	#include "TH3F.h"
	30	#include "TH2D.h"
	31	#include "TH3D.h"
	32	#include "TArrayF.h"
	33	#include "TF1.h"
	34	#include "TROOT.h"
	35	#include "stdio.h"
	36	#include "TCutG.h"
	37
	38
	39	#include "AliTHn.h"
	40	#include "AliLog.h"
	41	#include "AliAnalysisManager.h"
	42	#include "AliESDEvent.h"
	43	#include "AliAODInputHandler.h"
	44	#include "AliAODEvent.h"
	45	#include "AliAODTrack.h"
	46	#include "AliAODInputHandler.h"
	47	#include "AliCollisionGeometry.h"
	48	#include "AliGenEventHeader.h"
	49	#include "AliAnalysisUtils.h"
	50	#include "AliPIDCombined.h"
	51	#include "AliAnalysisTask.h"
	52	#include "AliAODHandler.h"
	53	#include <AliInputEventHandler.h>
	54	#include <AliVEventHandler.h>
	55	#include <AliVParticle.h>
	56	#include <AliVTrack.h>
	57	#include <AliTPCPIDResponse.h>
	58	#include <AliTOFPIDResponse.h>
	59	#include "AliAnalysisTaskPIDconfig.h"
	60	#include "AliAnalysisTaskSE.h"
	61	#include "AliAODPid.h"
	62	#include "AliPhysicsSelection.h"
	63	#include "AliCentralitySelectionTask.h"
	64	#include "AliCentrality.h"
65	#include "AliKFParticle.h"
66	#include "AliKFVertex.h"
67	#include "AliPID.h"
68	#include "AliPIDResponse.h"
69	#include "AliCFContainer.h"
70	#include "AliCFManager.h"
71	#include "AliVEvent.h"
72	#include "AliAODVZERO.h"
73
5a1b363f	74	using std::cout;
5a1b363f	75	using std::endl;
2c03c485	76
	77	ClassImp(AliAnalysisTaskPIDconfig)
	78	//ClassImp()
	79	//___________________________________________________________________
	80	AliAnalysisTaskPIDconfig::AliAnalysisTaskPIDconfig():
	81	AliAnalysisTaskSE(),
	82	fVevent(0),
	83	fESD(0),
	84	fAOD(0),
	85	fPIDResponse(0),
	86	fTriggerSelection(0),
	87	fCentralityPercentileMin(0.),
	88	fCentralityPercentileMax(5.),
	89	fFilterBit(128),
	90	fDCAxyCut(-1),
	91	fDCAzCut(-1),
	92	fData2011(kFALSE),
	93	fTriggerMB(kTRUE),
	94	fTriggerCentral(kFALSE),
	95	fUseCentrality(kTRUE),
	96	fCutTPCmultiplicityOutliersAOD(kFALSE),
	97	fPIDcuts(kFALSE),
	98	fCentralityEstimator("V0M"),
	99	fContourCutList(0),
	100	fListQA(0x0),
	101	fListQAtpctof(0x0),
	102	fListQAInfo(0x0),
	103	fhistCentralityPass(0),
	104	fNoEvents(0),
	105	fpVtxZ(0),
	106	fhistDCABefore(0),
	107	fhistDCAAfter(0),
	108	fhistPhiDistBefore(0),
	109	fhistPhiDistAfter(0),
	110	fhistEtaDistBefore(0),
	111	fhistEtaDistAfter(0),
	112	fTPCvsGlobalMultBeforeOutliers(0),
	113	fTPCvsGlobalMultAfterOutliers(0),
	114	fTPCvsGlobalMultAfter(0),
	115	fHistBetavsPTOFbeforePID(0),
	116	fHistdEdxVsPTPCbeforePID(0),
	117	fHistBetavsPTOFafterPID(0),
	118	fHistdEdxVsPTPCafterPID(0),
	119	fhistNsigmaP(0),
	120	fhistNsigmaPt(0)
	121	//fSparseSpecies(0),
	122	//fvalueSpecies(0)
	123	{
	124	//Dummy Constructor
	125	}
	126
	127
	128	//
	129
	130	AliAnalysisTaskPIDconfig::AliAnalysisTaskPIDconfig(const char *name):
	131	AliAnalysisTaskSE(name),
	132	fVevent(0),
	133	fESD(0),
	134	fAOD(0),
	135	fPIDResponse(0),
	136	fTriggerSelection(0),
	137	fCentralityPercentileMin(0.),
	138	fCentralityPercentileMax(5.),
	139	fFilterBit(1),
140	fDCAxyCut(-1),
141	fDCAzCut(-1),
142	fData2011(kFALSE),
143	fTriggerMB(kTRUE),
144	fTriggerCentral(kFALSE),
145	fUseCentrality(kTRUE),
146	fCutTPCmultiplicityOutliersAOD(kFALSE),
147	fPIDcuts(kFALSE),
148	fCentralityEstimator("V0M"),
149	fContourCutList(0),
150	fListQA(0x0),
151	fListQAtpctof(0x0),
152	fListQAInfo(0x0),
153	fhistCentralityPass(0),
154	fNoEvents(0),
155	fpVtxZ(0),
156	fhistDCABefore(0),
157	fhistDCAAfter(0),
158	fhistPhiDistBefore(0),
159	fhistPhiDistAfter(0),
160	fhistEtaDistBefore(0),
161	fhistEtaDistAfter(0),
162	fTPCvsGlobalMultBeforeOutliers(0),
163	fTPCvsGlobalMultAfterOutliers(0),
164	fTPCvsGlobalMultAfter(0),
165	fHistBetavsPTOFbeforePID(0),
166	fHistdEdxVsPTPCbeforePID(0),
167	fHistBetavsPTOFafterPID(0),
168	fHistdEdxVsPTPCafterPID(0),
169	fhistNsigmaP(0),
170	fhistNsigmaPt(0)
171	//fSparseSpecies(0),
172	//fvalueSpecies(0)
173	{
174	//fvalueSpecies = new Double_t[9];
175	//Default Constructor
176	DefineInput(0,TChain::Class());
177	DefineOutput(1,TList::Class());
178	}
179
180	//_____________________________________________________________________
181	AliAnalysisTaskPIDconfig::~AliAnalysisTaskPIDconfig()
182	{
183	//Destructor
184
185	// delete fPID;
186	// delete fPIDqa;
187	// delete fTrackCuts;
188	// delete fSparseSpecies;
189	//delete []fvalueSpecies;
190
191
192	}
193	//______________________________________________________________________
194	void AliAnalysisTaskPIDconfig::UserCreateOutputObjects()
195	{
196	//
197	// Create the output QA objects
198	//
199
200	AliLog::SetClassDebugLevel("AliAnalysisTaskPIDconfig",10);
201
202	//input hander
203	AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
204	AliInputEventHandler inputHandler=dynamic_cast<AliInputEventHandler>(mgr->GetInputEventHandler());
205	if (!inputHandler) AliFatal("Input handler needed");
206
207	//pid response object
208	fPIDResponse=inputHandler->GetPIDResponse();
209	if (!fPIDResponse) AliError("PIDResponse object was not created");
210
211	//
212	fListQA=new TList;
213	fListQA->SetOwner();
214
215	fListQAtpctof=new TList;
216	fListQAtpctof->SetOwner();
217	fListQAtpctof->SetName("PID_TPC_TOF");
218
219	fListQAInfo=new TList;
220	fListQAInfo->SetOwner();
221	fListQAInfo->SetName("Event_Track_Info");
222
223	fListQA->Add(fListQAtpctof);
224	fListQA->Add(fListQAInfo);
225
226	SetupTPCTOFqa();
227	SetupEventInfo();
228
229	PostData(1,fListQA);
230	}
231	//______________________________________________________________________
232	void AliAnalysisTaskPIDconfig::UserExec(Option_t*){
233	//Main loop
234	//Called for each event
235
236	// create pointer to event
237	fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
238	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
239
240	Int_t ntracks=fAOD->GetNumberOfTracks();
241
242
243	if(!(fESD \|\| fAOD)){
244	printf("ERROR: fESD & fAOD not available\n");
245	return;
246	}
247	fVevent = dynamic_cast<AliVEvent*>(InputEvent());
248	if (!fVevent) {
249	printf("ERROR: fVevent not available\n");
250	return;
251	}
252
253	Bool_t pass = kFALSE;
254
255	CheckCentrality(fVevent,pass);
256
257	if(!pass){ return;}
258
259	const AliVVertex *pVtx = fVevent->GetPrimaryVertex();
260
261	Double_t pVtxZ = -999;
262	pVtxZ = pVtx->GetZ();
263
264	if(TMath::Abs(pVtxZ)>10) return;
265
266	TH1F hNoEvents = (TH1F)fListQAInfo->At(1);
267	TH1F histpVtxZ = (TH1F)fListQAInfo->At(2);
268
269	if(hNoEvents) hNoEvents->Fill(0);
270	if(histpVtxZ) histpVtxZ->Fill(pVtxZ);
271
272	if(ntracks<2) return;
273
274	// if(!pass) return;
275
276	TH2F HistTPCvsGlobalMultBeforeOutliers = (TH2F)fListQAInfo->At(3);
277
278	TH2F HistTPCvsGlobalMultAfterOutliers = (TH2F)fListQAInfo->At(4);
279
280
281	Float_t multTPC(0.); // tpc mult estimate
282	Float_t multGlobal(0.); // global multiplicity
283
284	const Int_t nGoodTracks = fVevent->GetNumberOfTracks();
285	if(!fData2011) { // cut on outliers
286
287	for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill tpc mult
288	AliAODTrack* AODtrack =dynamic_cast<AliAODTrack*>(fVevent->GetTrack(iTrack));
289	if (!AODtrack) continue;
290	if (!(AODtrack->TestFilterBit(1))) continue;
291	if ((AODtrack->Pt() < .2) \|\| (AODtrack->Pt() > 5.0) \|\| (TMath::Abs(AODtrack->Eta()) > .8) \|\| (AODtrack->GetTPCNcls() < 70) \|\| (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) \|\| (AODtrack->Chi2perNDF() < 0.2)) continue;
292	multTPC++;
293	}//track loop
294
295	for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill global mult
296	AliAODTrack AODtrack=dynamic_cast<AliAODTrack>(fVevent->GetTrack(iTrack));
297	if (!AODtrack) continue;
298	if (!(AODtrack->TestFilterBit(16))) continue;
299	if ((AODtrack->Pt() < .2) \|\| (AODtrack->Pt() > 5.0) \|\| (TMath::Abs(AODtrack->Eta()) > .8) \|\| (AODtrack->GetTPCNcls() < 70) \|\| (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) \|\| (AODtrack->Chi2perNDF() < 0.1)) continue;
300	Double_t b[2] = {-99., -99.};
301	Double_t bCov[3] = {-99., -99., -99.};
302	AliAODTrack copy(*AODtrack);
303	if (!(copy.PropagateToDCA(fVevent->GetPrimaryVertex(), fVevent->GetMagneticField(), 100., b, bCov))) continue;
304	if ((TMath::Abs(b[0]) > 0.3) \|\| (TMath::Abs(b[1]) > 0.3)) continue;
305	multGlobal++;
306	} //track loop
307
308	HistTPCvsGlobalMultBeforeOutliers->Fill(multGlobal,multTPC);
309
310	if(multTPC < (-40.3+1.22multGlobal) \|\| multTPC > (32.1+1.59multGlobal)){ pass = kFALSE;}
311
312	if(!pass) return;
313	HistTPCvsGlobalMultAfterOutliers->Fill(multGlobal,multTPC);
314
315	}
316
317
318	if(fData2011) { // cut on outliers
319	//Float_t multTPC(0.); // tpc mult estimate
320	//Float_t multGlob(0.); // global multiplicity
321	for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill tpc mult
322	AliAODTrack AODtrack=dynamic_cast<AliAODTrack>(fVevent->GetTrack(iTrack));
323	if (!AODtrack) continue;
324	if (!(AODtrack->TestFilterBit(1))) continue;
325	if ((AODtrack->Pt() < .2) \|\| (AODtrack->Pt() > 5.0) \|\| (TMath::Abs(AODtrack->Eta()) > .8) \|\| (AODtrack->GetTPCNcls() < 70) \|\| (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) \|\| (AODtrack->Chi2perNDF() < 0.2)) continue;
326	multTPC++;
327	}
328	for(Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { // fill global mult
329	AliAODTrack AODtrack=dynamic_cast<AliAODTrack>(fVevent->GetTrack(iTrack));
330	if (!AODtrack) continue;
331	if (!(AODtrack->TestFilterBit(16))) continue;
332	if ((AODtrack->Pt() < .2) \|\| (AODtrack->Pt() > 5.0) \|\| (TMath::Abs(AODtrack->Eta()) > .8) \|\| (AODtrack->GetTPCNcls() < 70) \|\| (AODtrack->GetDetPid()->GetTPCsignal() < 10.0) \|\| (AODtrack->Chi2perNDF() < 0.1)) continue;
333	Double_t b[2] = {-99., -99.};
334	Double_t bCov[3] = {-99., -99., -99.};
335	AliAODTrack copy(*AODtrack);
336	if (!(copy.PropagateToDCA(fVevent->GetPrimaryVertex(), fVevent->GetMagneticField(), 100., b, bCov))) continue;
337	if ((TMath::Abs(b[0]) > 0.3) \|\| (TMath::Abs(b[1]) > 0.3)) continue;
338	multGlobal++;
339
340	} //track loop
341
342	HistTPCvsGlobalMultBeforeOutliers->Fill(multGlobal,multTPC);
343
344	if(multTPC < (-36.73 + 1.48multGlobal) \|\| multTPC > (62.87 + 1.78multGlobal)){pass = kFALSE;}
345
346	if(!pass) return;
347	HistTPCvsGlobalMultAfterOutliers->Fill(multGlobal,multTPC);
348
349	}
350
351
352
353	for(Int_t itrack = 0; itrack < ntracks; itrack++){
354
355	AliAODTrack track=dynamic_cast<AliAODTrack>(fVevent->GetTrack(itrack));
356	if(!track) continue;
357
358	Float_t dcaXY = track->DCA();
359	Float_t dcaZ = track->ZAtDCA();
360
361	TH2F* HistDCAbefore =(TH2F*)fListQAInfo->At(5);
362	HistDCAbefore->Fill(dcaZ,dcaXY);
363
d802b6ea	364	Double_t p = -999, /pTPC = -999,/ pT = -999, phi = -999, eta = -999, dEdx =-999;
2c03c485	365	Double_t length = -999., beta =-999, tofTime = -999., tof = -999.;
	366	Double_t c = TMath::C()*1.E-9;// m/ns
	367
	368	//cout<<"track->GetFilterMap()= "<<track->GetFilterMap()<<endl;
	369	if(!track->TestFilterBit(fFilterBit)) continue;
	370
	371	//Float_t dcaXY = -999, dcaZ = -999;
	372	p=track->P();
d802b6ea	373	/pTPC=track->GetTPCmomentum();/
2c03c485	374	pT=track->Pt();
	375	phi=track->Phi();
	376	eta=track->Eta();
	377	dEdx=track->GetDetPid()->GetTPCsignal();
	378
	379	if ( (track->IsOn(AliAODTrack::kTOFin)) &&
	380	(track->IsOn(AliAODTrack::kTIME)) && (track->IsOn(AliAODTrack::kTOFout))) {
	381	// if ( (track->IsOn(AliAODTrack::kTOFin)) &&
	382	// (track->IsOn(AliAODTrack::kTOFout)) ) {
	383
	384	tofTime = track->GetTOFsignal();//in ps
	385	length = track->GetIntegratedLength();
	386
	387	tof = tofTime*1E-3; // ns
	388	//cout<<"tof = "<<tof<<endl;
	389	if (tof <= 0)continue;
	390	//cout<<"length = "<<length<<endl;
	391	if (length <= 0) continue;
	392
	393	length = length*0.01; // in meters
	394	tof = tof*c;
	395	beta = length/tof;
	396
	397	TH2F HistBetavsPTOFbeforePID = (TH2F)fListQAInfo->At(6);
	398	HistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);
	399	}//TOF signal
	400
	401	TH2F HistdEdxVsPTPCbeforePID = (TH2F)fListQAInfo->At(7);
	402	HistdEdxVsPTPCbeforePID -> Fill(p*track->Charge(),dEdx); //TPC signal
	403
	404
	405	//QA plot
	406	TH1F HistPhiDistBefore = (TH1F)fListQAInfo->At(8);
	407	HistPhiDistBefore->Fill(phi);
	408	//
	409	TH1F HistEtaDistBefore = (TH1F)fListQAInfo->At(9);
	410	HistEtaDistBefore->Fill(eta);
	411
	412
	413	if(pT<0.1) continue;
	414	if(TMath::Abs(eta)>0.8) continue;
	415
	416	Int_t TPCNcls = track->GetTPCNcls();
	417
	418	if(TPCNcls<70 \|\| dEdx<10) continue;
	419
	420	// fill QA histograms
	421
	422	TH2F* HistDCAAfter =(TH2F*)fListQAInfo->At(10);
	423	HistDCAAfter->Fill(dcaZ,dcaXY);
	424
	425	TH1F HistPhiDistAfter = (TH1F)fListQAInfo->At(11);
	426	HistPhiDistAfter->Fill(phi);
	427
	428	TH1F HistEtaDistAfter = (TH1F)fListQAInfo->At(12);
	429	HistEtaDistAfter->Fill(eta);
	430
	431	Bool_t pWithinRange = kFALSE;
	432	Int_t pRange = -999;
	433	TCutG *cut[3][10];
	434	if(fPIDcuts){
d802b6ea	435	if(fContourCutList){cout<<"The contour file has been retrieved"<<endl;}
d802b6ea	436
2c03c485	437	TGraph *ContourCut[3][10];
d802b6ea	438	//Double_t plow[10] = {0.2,0.5,1,1.5,2,2.5,3,3.5,4,4.5};
	439	//Double_t phigh[10] = {0.5,1,1.5,2,2.5,3,3.5,4,4.5,5};
	440	Double_t phigh=0 ,plow=0;
	441	Double_t pBins[11] = {0.2,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5};
	442	TString Graph_Name = "contourlines_";
	443	TString species[3] = {"pion","kaon","proton"};
2c03c485	444	for(int i=0;i<3;i++){
2c03c485	445	for(int j=0;j<10;j++){
d802b6ea	446
	447	if(p>pBins[j] && p<pBins[j+1]){
	448
	449	//if(p>plow[j] && p<phigh[j]){
2c03c485	450	pWithinRange = kTRUE;
2c03c485	451	pRange = j;
d802b6ea	452
	453	TList Species_contours = (TList)fContourCutList->Get(species[i]);
	454	phigh = pBins[j+1]*10;
	455	plow = pBins[j]*10;
	456
	457	Graph_Name += species[i];
	458	Graph_Name += Form("%.f%.f-%i%icent",plow,phigh,fCentralityPercentileMin,fCentralityPercentileMax);
	459
	460	//TList Pcontours = (TList)Scontours->At(j);
	461	//if (!Scontours) return;
	462
	463	// if (!Pcontours \|\| !Scontours) return;
	464	ContourCut[i][j] = (TGraph*)Species_contours->FindObject(Graph_Name);
2c03c485	465	cut[i][j] = new TCutG("cut",ContourCut[i][j]->GetN(),ContourCut[i][j]->GetX(),ContourCut[i][j]->GetY());
	466	//ContourCut[i][j] = (TGraph )fContourCutList->At(10i+j);
	467	}
	468	}
	469	}
	470	}
	471
	472	for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
	473	//TOF nSigma
	474	Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
	475	Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
	476	if(fPIDcuts && ispecie>1 && ispecie<5 && pWithinRange){// for pions, kaons and protons only
	477	if(cut[ispecie-2][pRange]->IsInside(nSigmaTOF,nSigmaTPC)){
	478	pass = kTRUE;
	479	}
	480	else{
	481	pass = kFALSE;
	482	continue;
	483	}
	484	}
	485
	486	//TPC and TOF cuts, TPC TOF nsigma vs. momentum
	487	if(pass){
	488	TH3 hist1 = (TH3)fListQAtpctof->At(ispecie);
	489	if (hist1){
	490	hist1->Fill(nSigmaTPC,nSigmaTOF,p);}
	491
	492	TH3 hist2 = (TH3)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);
	493	if (hist2){
	494	hist2->Fill(nSigmaTPC,nSigmaTOF,pT);}
	495	}
	496
	497	}
	498
	499
	500	}//track loop
	501
	502	TH2F HistTPCvsGlobalMultAfter = (TH2F) fListQAInfo->At(13);
	503	HistTPCvsGlobalMultAfter->Fill(multGlobal,multTPC);
	504
	505	}
	506	//_________________________________________
	507	void AliAnalysisTaskPIDconfig::CheckCentrality(AliVEvent* event, Bool_t &centralitypass)
	508	{
	509	// Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
	510	if (!fUseCentrality) AliFatal("No centrality method set! FATAL ERROR!");
	511	Double_t centvalue = event->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
	512	//cout << "Centrality evaluated-------------------------: " << centvalue <<endl;
	513	if ((centvalue >= fCentralityPercentileMin) && (centvalue < fCentralityPercentileMax))
	514	{
	515	TH1F hCentralityPass = (TH1F)fListQAInfo->At(0);
	516	hCentralityPass->Fill(centvalue);
	517	//cout << "--------------Fill pass-------------------------"<<endl;
	518	centralitypass = kTRUE;
	519	}
	520
	521	}
	522	//______________________________________________________________________________
	523	void AliAnalysisTaskPIDconfig::SetupTPCTOFqa()
	524	{
	525	//
	526	// Create the qa objects for TPC + TOF combination
	527
	528
529	//TPC and TOF signal vs. momentum
530	for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
531	fhistNsigmaP = new TH3F(Form("NsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),Form("TPC n#sigma vs. TOF n#sigma %s vs. p ;TPC n#sigma;TOF n#sigma;p [GeV]",AliPID::ParticleName(ispecie)),200,-20,20,200,-20,20,60,0.1,6);
532	fListQAtpctof->Add(fhistNsigmaP);
533	}
534	//TPC and TOF signal vs. transverse momentum
535	for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
536	fhistNsigmaPt = new TH3F(Form("NsigmaPt_TPC_TOF_%s",AliPID::ParticleName(ispecie)),Form("TPC n#sigma vs. TOF n#sigma %s vs. Pt ;TPC n#sigma;TOF n#sigma;pT [GeV]",AliPID::ParticleName(ispecie)),200,-20,20,200,-20,20,60,0.1,6);
537	fListQAtpctof->Add(fhistNsigmaPt);
538	}
539
540	}
541	//______________________________________________________________________________
542	void AliAnalysisTaskPIDconfig::SetupEventInfo()
543	{
544	//event and track info
545
546	fhistCentralityPass = new TH1F("fcentralityPass","centralityPass", 100,0,100);
547	fListQAInfo->Add(fhistCentralityPass);
548
549	fNoEvents = new TH1F("number of events","no. of events",1,0,1);
550	fListQAInfo->Add(fNoEvents);
551
552	fpVtxZ = new TH1F("pVtxZ","pVtxZ",100,-20,20);
553	fListQAInfo->Add(fpVtxZ);
554
555	fTPCvsGlobalMultBeforeOutliers = new TH2F("TPC vs. Global Multiplicity Before","TPC vs. Global Multiplicity Before",500,0,6000,500,0,6000);
556	fListQAInfo->Add(fTPCvsGlobalMultBeforeOutliers);
557
558	fTPCvsGlobalMultAfterOutliers = new TH2F("TPC vs. Global Multiplicity After outliers","TPC vs. Global Multiplicity After outliers",500,0,6000,500,0,6000);
559	fListQAInfo->Add(fTPCvsGlobalMultAfterOutliers);
560
561	fhistDCABefore = new TH2F("DCA xy vs z (before)","DCA before",200,0,10,200,0,10);
562	fListQAInfo->Add(fhistDCABefore);
563
564	fHistBetavsPTOFbeforePID = new TH2F("momentum vs beta before PID","momentum vs beta before PID",1000,-10.,10.,1000,0,1.2);
565	fListQAInfo->Add(fHistBetavsPTOFbeforePID);
566
567	fHistdEdxVsPTPCbeforePID = new TH2F("momentum vs dEdx before PID","momentum vs dEdx before PID",1000,-10.,10.,1000,0,1000);
568	fListQAInfo->Add(fHistdEdxVsPTPCbeforePID);
569
570	fhistPhiDistBefore = new TH1F("Phi Distribution Before Cuts","Phi Distribution Before Cuts",200,0,6.4);
571	fListQAInfo->Add(fhistPhiDistBefore);
572
573	fhistEtaDistBefore = new TH1F("Eta Distribution Before Cuts","Eta Distribution Before Cuts",200,-10,10);
574	fListQAInfo->Add(fhistEtaDistBefore);
575
576	fhistDCAAfter = new TH2F("DCA xy vs z (after)","DCA after",200,0,10,200,0,10);
577	fListQAInfo->Add(fhistDCAAfter);
578
579	fhistPhiDistAfter = new TH1F("Phi Distribution After Cuts","Phi Distribution After Cuts",200,0,6.4);
580	fListQAInfo->Add(fhistPhiDistAfter);
581
582	fhistEtaDistAfter = new TH1F("Eta Distribution After Cuts","Eta Distribution After Cuts",200,-10,10);
583	fListQAInfo->Add(fhistEtaDistAfter);
584
585	fTPCvsGlobalMultAfter = new TH2F("TPC vs. Global Multiplicity After","TPC vs. Global Multiplicity After",500,0,6000,500,0,6000);
586	fListQAInfo->Add(fTPCvsGlobalMultAfter);
587
588	// fHistBetavsPTOFafterPID = new TH2F("momentum vs beta after PID","momentum vs beta after PID",1000,-10.,10.,1000,0,1.2);
589	// fListQAInfo->Add(fHistBetavsPTOFafterPID);
590
591	// fHistdEdxVsPTPCafterPID = new TH2F("momentum vs dEdx after PID","momentum vs dEdx after PID",1000,-10.,10.,1000,0,1000);
592	// fListQAInfo->Add(fHistdEdxVsPTPCafterPID);
593
594	}
595
596