[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskElecV2.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.                  *               
 **************************************************************************/

// Class for heavy-flavour electron v2 with EMCal triggered events
// Author: Denise Godoy


#include "TChain.h"
#include "TTree.h"
#include "TH2F.h"
#include "TMath.h"
#include "TCanvas.h"
#include "THnSparse.h"
#include "TLorentzVector.h"
#include "TString.h"
#include "TFile.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDEvent.h"
#include "AliESDHandler.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"

#include "AliAnalysisTaskElecV2.h"
#include "TGeoGlobalMagField.h"
#include "AliLog.h"
#include "AliAnalysisTaskSE.h"
#include "TRefArray.h"
#include "TVector.h"
#include "AliESDInputHandler.h"
#include "AliESDpid.h"
#include "AliESDtrackCuts.h"
#include "AliPhysicsSelection.h"
#include "AliESDCaloCluster.h"
#include "AliAODCaloCluster.h"
#include "AliEMCALRecoUtils.h"
#include "AliEMCALGeometry.h"
#include "AliGeomManager.h"
#include "stdio.h"
#include "TGeoManager.h"
#include "iostream"
#include "fstream"

#include "AliEMCALTrack.h"
#include "AliMagF.h"

#include "AliKFParticle.h"
#include "AliKFVertex.h"

#include "AliPID.h"
#include "AliPIDResponse.h"
#include "AliHFEcontainer.h"
#include "AliHFEcuts.h"
#include "AliHFEpid.h"
#include "AliHFEpidBase.h"
#include "AliHFEpidQAmanager.h"
#include "AliHFEtools.h"
#include "AliCFContainer.h"
#include "AliCFManager.h"

#include "AliEventplane.h"
#include "AliCentrality.h"

ClassImp(AliAnalysisTaskElecV2)
//________________________________________________________________________
AliAnalysisTaskElecV2::AliAnalysisTaskElecV2(const char *name) 
  : AliAnalysisTaskSE(name)
  ,fESD(0)
  ,fOutputList(0)
  ,fTrackCuts(0)
  ,fCuts(0)
  ,fIdentifiedAsOutInz(kFALSE)
  ,fPassTheEventCut(kFALSE)
  ,fRejectKinkMother(kFALSE)
  ,fVz(0.0)
  ,fCFM(0)	
  ,fPID(0)
  ,fPIDqa(0)	       
  ,fOpeningAngleCut(0.1)
  ,fInvmassCut(0.01)	
  ,fNoEvents(0)
  ,fTrkpt(0)
  ,fTrkEovPBef(0)	 
  ,fTrkEovPAft(0)	
  ,fdEdxBef(0)	 
  ,fdEdxAft(0)	 
  ,fInvmassLS(0)		
  ,fInvmassULS(0)		
  ,fOpeningAngleLS(0)	
  ,fOpeningAngleULS(0)	
  ,fPhotoElecPt(0)
  ,fSemiInclElecPt(0)
  ,fTrackPtBefTrkCuts(0)	 
  ,fTrackPtAftTrkCuts(0)
  ,fTPCnsigma(0)
  ,fCent(0)
  ,fTPCsubEPres(0)
  ,fEPres(0)
  ,fCorr(0)
  ,feTPCV2(0)
  ,feV2(0)
  ,fphoteV2(0)
  ,fChargPartV2(0)
{
  //Named constructor
  
  fPID = new AliHFEpid("hfePid");
  fTrackCuts = new AliESDtrackCuts();
  
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 id reserved by the base class for AOD
  // Output slot #1 writes into a TH1 container
  // DefineOutput(1, TH1I::Class());
  DefineOutput(1, TList::Class());
  //  DefineOutput(3, TTree::Class());
}

//________________________________________________________________________
AliAnalysisTaskElecV2::AliAnalysisTaskElecV2() 
  : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
  ,fESD(0)
  ,fOutputList(0)
  ,fTrackCuts(0)
  ,fCuts(0)
  ,fIdentifiedAsOutInz(kFALSE)
  ,fPassTheEventCut(kFALSE)
  ,fRejectKinkMother(kFALSE)
  ,fVz(0.0)
  ,fCFM(0)	
  ,fPID(0)       
  ,fPIDqa(0)	       
  ,fOpeningAngleCut(0.1)
  ,fInvmassCut(0.01)	
  ,fNoEvents(0)
  ,fTrkpt(0)
  ,fTrkEovPBef(0)	 
  ,fTrkEovPAft(0)	 
  ,fdEdxBef(0)	 
  ,fdEdxAft(0)	 
  ,fInvmassLS(0)		
  ,fInvmassULS(0)		
  ,fOpeningAngleLS(0)	
  ,fOpeningAngleULS(0)	
  ,fPhotoElecPt(0)
  ,fSemiInclElecPt(0)
  ,fTrackPtBefTrkCuts(0)	 
  ,fTrackPtAftTrkCuts(0)	 	  
  ,fTPCnsigma(0)
  ,fCent(0)
  ,fTPCsubEPres(0)
  ,fEPres(0)
  ,fCorr(0)
  ,feTPCV2(0)
  ,feV2(0)
  ,fphoteV2(0)
  ,fChargPartV2(0)
{
	//Default constructor
	fPID = new AliHFEpid("hfePid");

	fTrackCuts = new AliESDtrackCuts();
	
	// Constructor
	// Define input and output slots here
	// Input slot #0 works with a TChain
	DefineInput(0, TChain::Class());
	// Output slot #0 id reserved by the base class for AOD
	// Output slot #1 writes into a TH1 container
	// DefineOutput(1, TH1I::Class());
	DefineOutput(1, TList::Class());
	//DefineOutput(3, TTree::Class());
}
//_________________________________________

AliAnalysisTaskElecV2::~AliAnalysisTaskElecV2()
{
  //Destructor 
  
  delete fOutputList;
  delete fPID;
  delete fCFM;
  delete fPIDqa;
  delete fTrackCuts;
}
//_________________________________________

void AliAnalysisTaskElecV2::UserExec(Option_t*)
{
  //Main loop
  //Called for each event
  
  // create pointer to event
  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!fESD) {
    printf("ERROR: fESD not available\n");
    return;
  }
  
  if(!fCuts){
    AliError("HFE cuts not available");
    return;
  }
  
  if(!fPID->IsInitialized()){ 
    // Initialize PID with the given run number
    AliWarning("PID not initialised, get from Run no");
    fPID->InitializePID(fESD->GetRunNumber());
  }
 
 
  Int_t fNOtrks =  fESD->GetNumberOfTracks();
  const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
  
  Double_t pVtxZ = -999;
  pVtxZ = pVtx->GetZ();
  
  if(TMath::Abs(pVtxZ)>10) return;
  fNoEvents->Fill(0);
  
  if(fNOtrks<2) return;
  
  AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
  if(!pidResponse){
    AliDebug(1, "Using default PID Response");
    pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class()); 
  }
  
  fPID->SetPIDResponse(pidResponse);
  
  fCFM->SetRecEventInfo(fESD);
  
  Float_t cent = -1.;
  AliCentrality *centrality = fESD->GetCentrality(); 
  cent = centrality->GetCentralityPercentile("V0M");
  fCent->Fill(cent);
  
  if(cent>90.) return;
	
  //Event planes
  
  Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
  if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
  
  Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
  if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
  
  AliEventplane* esdTPCep = fESD->GetEventplane();
  TVector2 *standardQ = esdTPCep->GetQVector(); 
  Double_t qx = -999., qy = -999.;
  if(standardQ)
  {
	  qx = standardQ->X();
	  qy = standardQ->Y();
  }
  TVector2 qVectorfortrack;
  qVectorfortrack.Set(qx,qy);
  Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;

  TVector2 *qsub1a = esdTPCep->GetQsub1();
  TVector2 *qsub2a = esdTPCep->GetQsub2();
  Double_t evPlaneResTPC = -999.;
  if(qsub1a && qsub2a)
  {
	  evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
  }
    
  fTPCsubEPres->Fill(evPlaneResTPC,cent);
  
  Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0A,evPlaneV0C),GetCos2DeltaPhi(evPlaneV0A,evPlaneTPC),GetCos2DeltaPhi(evPlaneV0C,evPlaneTPC),cent};
  fEPres->Fill(evPlaneRes);
  
  // Track loop 
  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
    AliESDtrack* track = fESD->GetTrack(iTracks);
    if (!track) {
      printf("ERROR: Could not receive track %d\n", iTracks);
      continue;
    }
    
    if(TMath::Abs(track->Eta())>0.7) continue;
    
    fTrackPtBefTrkCuts->Fill(track->Pt());		
    // RecKine: ITSTPC cuts  
    if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
    
    //RecKink
    if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
      if(track->GetKinkIndex(0) != 0) continue;
    } 
    
    // RecPrim
    if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
    
    // HFEcuts: ITS layers cuts
    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
    
    // HFE cuts: TPC PID cleanup
    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
    
    fTrackPtAftTrkCuts->Fill(track->Pt());		
    
    Double_t clsE = -999., p = -999., EovP=-999., pt = -999., dEdx=-999., fTPCnSigma=0, phi=-999.;
    
    // Track extrapolation
    
    pt = track->Pt();
    fTrkpt->Fill(pt);
    
    Int_t clsId = track->GetEMCALcluster();
    if (clsId>0){
      AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
      if(cluster && cluster->IsEMCAL()){
	clsE = cluster->E();
      }
    }
        
    p = track->P();
    phi = track->Phi();
    dEdx = track->GetTPCsignal();
    EovP = clsE/p;
    fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
    fdEdxBef->Fill(p,dEdx);
    fTPCnsigma->Fill(p,fTPCnSigma);
       
    Double_t corr[7]={fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C)};
    fCorr->Fill(corr);
       
    if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
    Int_t pidpassed = 1;
    
    //--- track accepted
    AliHFEpidObject hfetrack;
    hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
    hfetrack.SetRecTrack(track);
    hfetrack.SetPbPb();
    if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
    
    Double_t corrV2[6]={cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
    fChargPartV2->Fill(corrV2); 

    if(fTPCnSigma >= -0.5){
      	Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track); 
	Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track); 
	TVector2 newQVectorfortrack;
	newQVectorfortrack.Set(qX,qY);
	Double_t corrV2TPC = -999.;
	corrV2TPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2; 

	Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,corrV2TPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};

	feTPCV2->Fill(correctedV2);
    }
    
    if(pidpassed==0) continue;
    
    Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track); 
    Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track); 
    TVector2 newQVectorfortrack;
    newQVectorfortrack.Set(qX,qY);
    Double_t corrV2TPC = -999.;
    corrV2TPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2; 
        
    Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,corrV2TPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
    
    feV2->Fill(correctedV2);
    
    fTrkEovPAft->Fill(pt,EovP);
    fdEdxAft->Fill(p,dEdx);
    
    Bool_t fFlagPhotonicElec = kFALSE;
    SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
    
    if(fFlagPhotonicElec){
      fphoteV2->Fill(correctedV2);
      fPhotoElecPt->Fill(pt);
    }
    
    if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
 }
 PostData(1, fOutputList);
}
//_________________________________________
void AliAnalysisTaskElecV2::UserCreateOutputObjects()
{
  //--------Initialize PID
  fPID->SetHasMCData(kFALSE);
  if(!fPID->GetNumberOfPIDdetectors()) 
    {
      fPID->AddDetector("TPC", 0);
      fPID->AddDetector("EMCAL", 1);
    }
  
  fPID->SortDetectors(); 
  fPIDqa = new AliHFEpidQAmanager();
  fPIDqa->Initialize(fPID);
  
  //--------Initialize correction Framework and Cuts
  fCFM = new AliCFManager;
  const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
  fCFM->SetNStepParticle(kNcutSteps);
  for(Int_t istep = 0; istep < kNcutSteps; istep++)
    fCFM->SetParticleCutsList(istep, NULL);
  
  if(!fCuts){
    AliWarning("Cuts not available. Default cuts will be used");
    fCuts = new AliHFEcuts;
    fCuts->CreateStandardCuts();
  }
  fCuts->Initialize(fCFM);
  
  //---------Output Tlist
  fOutputList = new TList();
  fOutputList->SetOwner();
  fOutputList->Add(fPIDqa->MakeList("PIDQA"));
  
  fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
  fOutputList->Add(fNoEvents);
  
  fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
  fOutputList->Add(fTrkpt);
  
  fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
  fOutputList->Add(fTrackPtBefTrkCuts);
  
  fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
  fOutputList->Add(fTrackPtAftTrkCuts);
  
  fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
  fOutputList->Add(fTPCnsigma);
  
  fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
  fOutputList->Add(fTrkEovPBef);
  
  fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
  fOutputList->Add(fTrkEovPAft);
  
  fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
  fOutputList->Add(fdEdxBef);
  
  fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
  fOutputList->Add(fdEdxAft);
  
  fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
  fOutputList->Add(fInvmassLS);
  
  fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
  fOutputList->Add(fInvmassULS);
  
  fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
  fOutputList->Add(fOpeningAngleLS);
  
  fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
  fOutputList->Add(fOpeningAngleULS);
  
  fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
  fOutputList->Add(fPhotoElecPt);
  
  fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
  fOutputList->Add(fSemiInclElecPt);
  
  fCent = new TH1F("fCent","Centrality",100,0,100) ;
  fOutputList->Add(fCent);
  
  fTPCsubEPres = new TH2F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1,90,0,90);
  fOutputList->Add(fTPCsubEPres);
  
  Int_t binsv1[4]={100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, cent
  Double_t xminv1[4]={-1,-1,-1,0};
  Double_t xmaxv1[4]={1,1,1,90}; 
  fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
  fOutputList->Add(fEPres);
	
  Int_t binsv2[7]={100,90,100,100,100,100,100}; // fTPCnSigma,cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
  Double_t xminv2[7]={-3.5,0,0,0,0,0,0};
  Double_t xmaxv2[7]={3.5,90,50,3,TMath::Pi(),TMath::Pi(),TMath::Pi()}; 
  fCorr = new THnSparseD ("fCorr","Correlations",7,binsv2,xminv2,xmaxv2);
  fOutputList->Add(fCorr);
    
  Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
  Double_t xminv3[5]={0,0,-1,-1,-1};
  Double_t xmaxv3[5]={90,50,1,1,1}; 
  feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
  fOutputList->Add(feV2);
  
  Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
  Double_t xminv4[5]={0,0,-1,-1,-1};
  Double_t xmaxv4[5]={90,50,1,1,1}; 
  fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
  fOutputList->Add(fphoteV2);
  
  Int_t binsv5[6]={90,100,100,100,100,100}; // cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
  Double_t xminv5[6]={0,0,0,-1,-1,-1};
  Double_t xmaxv5[6]={90,50,3,1,1,1}; 
  fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",6,binsv5,xminv5,xmaxv5);
  fOutputList->Add(fChargPartV2);
  
  Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
  Double_t xminv6[5]={0,0,-1,-1,-1};
  Double_t xmaxv6[5]={90,50,1,1,1}; 
  feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
  fOutputList->Add(feTPCV2);
  
  PostData(1,fOutputList);
}

//________________________________________________________________________
void AliAnalysisTaskElecV2::Terminate(Option_t *)
{
  // Info("Terminate");
	AliAnalysisTaskSE::Terminate();
}

//________________________________________________________________________
Bool_t AliAnalysisTaskElecV2::ProcessCutStep(Int_t cutStep, AliVParticle *track)
{
  // Check single track cuts for a given cut step
  const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
  if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
  return kTRUE;
}
//_________________________________________
void AliAnalysisTaskElecV2::SelectPhotonicElectron(Int_t itrack, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
{
  //Identify non-heavy flavour electrons using Invariant mass method
  
  fTrackCuts->SetAcceptKinkDaughters(kFALSE);
  fTrackCuts->SetRequireTPCRefit(kTRUE);
  fTrackCuts->SetEtaRange(-0.7,0.7);
  fTrackCuts->SetRequireSigmaToVertex(kTRUE);
  fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
  fTrackCuts->SetMinNClustersTPC(100);
  
  const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
  
  Bool_t flagPhotonicElec = kFALSE;
  
  for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
    AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
    if (!trackAsso) {
      printf("ERROR: Could not receive track %d\n", jTracks);
      continue;
    }
    
    Double_t dEdxAsso = -999., ptAsso=-999., openingAngle = -999.;
    Double_t mass=999., width = -999;
    Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
    
    dEdxAsso = trackAsso->GetTPCsignal();
    ptAsso = trackAsso->Pt();
    Int_t chargeAsso = trackAsso->Charge();
    Int_t charge = track->Charge();
    
    if(ptAsso <0.3) continue;
    if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
    if(dEdxAsso <70 || dEdxAsso>100) continue; //11a pass1
    
    Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
    if(charge>0) fPDGe1 = -11;
    if(chargeAsso>0) fPDGe2 = -11;
    
    if(charge == chargeAsso) fFlagLS = kTRUE;
    if(charge != chargeAsso) fFlagULS = kTRUE;
    
    AliKFParticle ge1(*track, fPDGe1);
    AliKFParticle ge2(*trackAsso, fPDGe2);
    AliKFParticle recg(ge1, ge2);
    
    if(recg.GetNDF()<1) continue;
    Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
    if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
    
    AliKFVertex primV(*pVtx);
    primV += recg;
    recg.SetProductionVertex(primV);
    
    recg.SetMassConstraint(0,0.0001);
    
    openingAngle = ge1.GetAngle(ge2);
    if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
    if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
    
    if(openingAngle > fOpeningAngleCut) continue;
    
    recg.GetMass(mass,width);
    
    if(fFlagLS) fInvmassLS->Fill(mass);
    if(fFlagULS) fInvmassULS->Fill(mass);
          
    if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
      flagPhotonicElec = kTRUE;
    }
    
  }
  fFlagPhotonicElec = flagPhotonicElec;
  
}
//_________________________________________
Double_t AliAnalysisTaskElecV2::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
{
  //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
  Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB); 
  if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
  Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
  
  return cos2DeltaPhi;
}

//_________________________________________
Double_t AliAnalysisTaskElecV2::GetDeltaPhi(Double_t phiA,Double_t phiB) const
{
  //Get phi-psi_EP
  Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB); 
  if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
  
  return dPhi;
}
Commit	Line	Data
1d08b6e8	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	// Class for heavy-flavour electron v2 with EMCal triggered events
	17	// Author: Denise Godoy
	18
	19
	20	#include "TChain.h"
	21	#include "TTree.h"
	22	#include "TH2F.h"
	23	#include "TMath.h"
	24	#include "TCanvas.h"
	25	#include "THnSparse.h"
	26	#include "TLorentzVector.h"
	27	#include "TString.h"
	28	#include "TFile.h"
	29
	30	#include "AliAnalysisTask.h"
	31	#include "AliAnalysisManager.h"
	32
	33	#include "AliESDEvent.h"
	34	#include "AliESDHandler.h"
	35	#include "AliAODEvent.h"
	36	#include "AliAODHandler.h"
	37
	38	#include "AliAnalysisTaskElecV2.h"
	39	#include "TGeoGlobalMagField.h"
	40	#include "AliLog.h"
	41	#include "AliAnalysisTaskSE.h"
	42	#include "TRefArray.h"
	43	#include "TVector.h"
	44	#include "AliESDInputHandler.h"
	45	#include "AliESDpid.h"
	46	#include "AliESDtrackCuts.h"
	47	#include "AliPhysicsSelection.h"
	48	#include "AliESDCaloCluster.h"
	49	#include "AliAODCaloCluster.h"
	50	#include "AliEMCALRecoUtils.h"
	51	#include "AliEMCALGeometry.h"
	52	#include "AliGeomManager.h"
	53	#include "stdio.h"
	54	#include "TGeoManager.h"
	55	#include "iostream"
	56	#include "fstream"
	57
	58	#include "AliEMCALTrack.h"
	59	#include "AliMagF.h"
	60
	61	#include "AliKFParticle.h"
	62	#include "AliKFVertex.h"
	63
	64	#include "AliPID.h"
65	#include "AliPIDResponse.h"
66	#include "AliHFEcontainer.h"
67	#include "AliHFEcuts.h"
68	#include "AliHFEpid.h"
69	#include "AliHFEpidBase.h"
70	#include "AliHFEpidQAmanager.h"
71	#include "AliHFEtools.h"
72	#include "AliCFContainer.h"
73	#include "AliCFManager.h"
74
75	#include "AliEventplane.h"
76	#include "AliCentrality.h"
77
78	ClassImp(AliAnalysisTaskElecV2)
79	//________________________________________________________________________
80	AliAnalysisTaskElecV2::AliAnalysisTaskElecV2(const char *name)
81	: AliAnalysisTaskSE(name)
82	,fESD(0)
83	,fOutputList(0)
84	,fTrackCuts(0)
85	,fCuts(0)
86	,fIdentifiedAsOutInz(kFALSE)
87	,fPassTheEventCut(kFALSE)
88	,fRejectKinkMother(kFALSE)
89	,fVz(0.0)
90	,fCFM(0)
91	,fPID(0)
92	,fPIDqa(0)
93	,fOpeningAngleCut(0.1)
94	,fInvmassCut(0.01)
95	,fNoEvents(0)
96	,fTrkpt(0)
97	,fTrkEovPBef(0)
98	,fTrkEovPAft(0)
99	,fdEdxBef(0)
100	,fdEdxAft(0)
101	,fInvmassLS(0)
102	,fInvmassULS(0)
103	,fOpeningAngleLS(0)
104	,fOpeningAngleULS(0)
105	,fPhotoElecPt(0)
106	,fSemiInclElecPt(0)
107	,fTrackPtBefTrkCuts(0)
108	,fTrackPtAftTrkCuts(0)
109	,fTPCnsigma(0)
110	,fCent(0)
111	,fTPCsubEPres(0)
112	,fEPres(0)
113	,fCorr(0)
114	,feTPCV2(0)
115	,feV2(0)
116	,fphoteV2(0)
117	,fChargPartV2(0)
118	{
119	//Named constructor
120
121	fPID = new AliHFEpid("hfePid");
122	fTrackCuts = new AliESDtrackCuts();
123
124	// Define input and output slots here
125	// Input slot #0 works with a TChain
126	DefineInput(0, TChain::Class());
127	// Output slot #0 id reserved by the base class for AOD
128	// Output slot #1 writes into a TH1 container
129	// DefineOutput(1, TH1I::Class());
130	DefineOutput(1, TList::Class());
131	// DefineOutput(3, TTree::Class());
132	}
133
134	//________________________________________________________________________
135	AliAnalysisTaskElecV2::AliAnalysisTaskElecV2()
136	: AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
137	,fESD(0)
138	,fOutputList(0)
139	,fTrackCuts(0)
140	,fCuts(0)
141	,fIdentifiedAsOutInz(kFALSE)
142	,fPassTheEventCut(kFALSE)
143	,fRejectKinkMother(kFALSE)
144	,fVz(0.0)
145	,fCFM(0)
146	,fPID(0)
147	,fPIDqa(0)
148	,fOpeningAngleCut(0.1)
149	,fInvmassCut(0.01)
150	,fNoEvents(0)
151	,fTrkpt(0)
152	,fTrkEovPBef(0)
153	,fTrkEovPAft(0)
154	,fdEdxBef(0)
155	,fdEdxAft(0)
156	,fInvmassLS(0)
157	,fInvmassULS(0)
158	,fOpeningAngleLS(0)
159	,fOpeningAngleULS(0)
160	,fPhotoElecPt(0)
161	,fSemiInclElecPt(0)
162	,fTrackPtBefTrkCuts(0)
163	,fTrackPtAftTrkCuts(0)
164	,fTPCnsigma(0)
165	,fCent(0)
166	,fTPCsubEPres(0)
167	,fEPres(0)
168	,fCorr(0)
169	,feTPCV2(0)
170	,feV2(0)
171	,fphoteV2(0)
172	,fChargPartV2(0)
173	{
174	//Default constructor
175	fPID = new AliHFEpid("hfePid");
176
177	fTrackCuts = new AliESDtrackCuts();
178
179	// Constructor
180	// Define input and output slots here
181	// Input slot #0 works with a TChain
182	DefineInput(0, TChain::Class());
183	// Output slot #0 id reserved by the base class for AOD
184	// Output slot #1 writes into a TH1 container
185	// DefineOutput(1, TH1I::Class());
186	DefineOutput(1, TList::Class());
187	//DefineOutput(3, TTree::Class());
188	}
189	//_________________________________________
190
191	AliAnalysisTaskElecV2::~AliAnalysisTaskElecV2()
192	{
193	//Destructor
194
195	delete fOutputList;
196	delete fPID;
197	delete fCFM;
198	delete fPIDqa;
199	delete fTrackCuts;
200	}
201	//_________________________________________
202
203	void AliAnalysisTaskElecV2::UserExec(Option_t*)
204	{
205	//Main loop
206	//Called for each event
207
208	// create pointer to event
209	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
210	if (!fESD) {
211	printf("ERROR: fESD not available\n");
212	return;
213	}
214
215	if(!fCuts){
216	AliError("HFE cuts not available");
217	return;
218	}
219
220	if(!fPID->IsInitialized()){
221	// Initialize PID with the given run number
222	AliWarning("PID not initialised, get from Run no");
223	fPID->InitializePID(fESD->GetRunNumber());
224	}
225
226
227	Int_t fNOtrks = fESD->GetNumberOfTracks();
228	const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
229
230	Double_t pVtxZ = -999;
231	pVtxZ = pVtx->GetZ();
232
233	if(TMath::Abs(pVtxZ)>10) return;
234	fNoEvents->Fill(0);
235
236	if(fNOtrks<2) return;
237
238	AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
239	if(!pidResponse){
240	AliDebug(1, "Using default PID Response");
241	pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
242	}
243
244	fPID->SetPIDResponse(pidResponse);
245
246	fCFM->SetRecEventInfo(fESD);
247
248	Float_t cent = -1.;
249	AliCentrality *centrality = fESD->GetCentrality();
250	cent = centrality->GetCentralityPercentile("V0M");
251	fCent->Fill(cent);
252
253	if(cent>90.) return;
254
255	//Event planes
256
257	Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
258	if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
259
260	Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
261	if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
262
263	AliEventplane* esdTPCep = fESD->GetEventplane();
264	TVector2 *standardQ = esdTPCep->GetQVector();
265	Double_t qx = -999., qy = -999.;
266	if(standardQ)
267	{
268	qx = standardQ->X();
269	qy = standardQ->Y();
270	}
271	TVector2 qVectorfortrack;
272	qVectorfortrack.Set(qx,qy);
273	Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
274
275	TVector2 *qsub1a = esdTPCep->GetQsub1();
276	TVector2 *qsub2a = esdTPCep->GetQsub2();
277	Double_t evPlaneResTPC = -999.;
278	if(qsub1a && qsub2a)
279	{
280	evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
281	}
282
283	fTPCsubEPres->Fill(evPlaneResTPC,cent);
284
285	Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0A,evPlaneV0C),GetCos2DeltaPhi(evPlaneV0A,evPlaneTPC),GetCos2DeltaPhi(evPlaneV0C,evPlaneTPC),cent};
286	fEPres->Fill(evPlaneRes);
287
288	// Track loop
289	for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
290	AliESDtrack* track = fESD->GetTrack(iTracks);
291	if (!track) {
292	printf("ERROR: Could not receive track %d\n", iTracks);
293	continue;
294	}
295
296	if(TMath::Abs(track->Eta())>0.7) continue;
297
298	fTrackPtBefTrkCuts->Fill(track->Pt());
299	// RecKine: ITSTPC cuts
300	if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
301
302	//RecKink
303	if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
304	if(track->GetKinkIndex(0) != 0) continue;
305	}
306
307	// RecPrim
308	if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
309
310	// HFEcuts: ITS layers cuts
311	if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
312
313	// HFE cuts: TPC PID cleanup
314	if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
315
316	fTrackPtAftTrkCuts->Fill(track->Pt());
317
318	Double_t clsE = -999., p = -999., EovP=-999., pt = -999., dEdx=-999., fTPCnSigma=0, phi=-999.;
319
320	// Track extrapolation
321
322	pt = track->Pt();
323	fTrkpt->Fill(pt);
324
325	Int_t clsId = track->GetEMCALcluster();
326	if (clsId>0){
327	AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
328	if(cluster && cluster->IsEMCAL()){
329	clsE = cluster->E();
330	}
331	}
332
333	p = track->P();
334	phi = track->Phi();
335	dEdx = track->GetTPCsignal();
336	EovP = clsE/p;
337	fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
338	fdEdxBef->Fill(p,dEdx);
339	fTPCnsigma->Fill(p,fTPCnSigma);
340
341	Double_t corr[7]={fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C)};
342	fCorr->Fill(corr);
343
344	if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
345	Int_t pidpassed = 1;
346
347	//--- track accepted
348	AliHFEpidObject hfetrack;
349	hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
350	hfetrack.SetRecTrack(track);
351	hfetrack.SetPbPb();
352	if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
353
354	Double_t corrV2[6]={cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
355	fChargPartV2->Fill(corrV2);
b3bd01b3	356
	357	if(fTPCnSigma >= -0.5){
	358	Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
	359	Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
	360	TVector2 newQVectorfortrack;
	361	newQVectorfortrack.Set(qX,qY);
	362	Double_t corrV2TPC = -999.;
	363	corrV2TPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
	364
	365	Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,corrV2TPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
	366
	367	feTPCV2->Fill(correctedV2);
	368	}
	369
	370	if(pidpassed==0) continue;
	371
1d08b6e8	372	Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
	373	Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
	374	TVector2 newQVectorfortrack;
	375	newQVectorfortrack.Set(qX,qY);
	376	Double_t corrV2TPC = -999.;
	377	corrV2TPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
	378
	379	Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,corrV2TPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
	380
1d08b6e8	381	feV2->Fill(correctedV2);
	382
	383	fTrkEovPAft->Fill(pt,EovP);
	384	fdEdxAft->Fill(p,dEdx);
	385
	386	Bool_t fFlagPhotonicElec = kFALSE;
	387	SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
	388
	389	if(fFlagPhotonicElec){
	390	fphoteV2->Fill(correctedV2);
	391	fPhotoElecPt->Fill(pt);
	392	}
	393
	394	if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
	395	}
	396	PostData(1, fOutputList);
	397	}
	398	//_________________________________________
	399	void AliAnalysisTaskElecV2::UserCreateOutputObjects()
	400	{
	401	//--------Initialize PID
	402	fPID->SetHasMCData(kFALSE);
	403	if(!fPID->GetNumberOfPIDdetectors())
	404	{
	405	fPID->AddDetector("TPC", 0);
	406	fPID->AddDetector("EMCAL", 1);
	407	}
	408
	409	fPID->SortDetectors();
	410	fPIDqa = new AliHFEpidQAmanager();
	411	fPIDqa->Initialize(fPID);
	412
	413	//--------Initialize correction Framework and Cuts
	414	fCFM = new AliCFManager;
	415	const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
	416	fCFM->SetNStepParticle(kNcutSteps);
	417	for(Int_t istep = 0; istep < kNcutSteps; istep++)
	418	fCFM->SetParticleCutsList(istep, NULL);
	419
	420	if(!fCuts){
	421	AliWarning("Cuts not available. Default cuts will be used");
	422	fCuts = new AliHFEcuts;
	423	fCuts->CreateStandardCuts();
	424	}
	425	fCuts->Initialize(fCFM);
	426
	427	//---------Output Tlist
	428	fOutputList = new TList();
	429	fOutputList->SetOwner();
	430	fOutputList->Add(fPIDqa->MakeList("PIDQA"));
	431
	432	fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
	433	fOutputList->Add(fNoEvents);
	434
	435	fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
	436	fOutputList->Add(fTrkpt);
	437
	438	fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
	439	fOutputList->Add(fTrackPtBefTrkCuts);
	440
	441	fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
	442	fOutputList->Add(fTrackPtAftTrkCuts);
	443
	444	fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
445	fOutputList->Add(fTPCnsigma);
446
447	fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
448	fOutputList->Add(fTrkEovPBef);
449
450	fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
451	fOutputList->Add(fTrkEovPAft);
452
453	fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
454	fOutputList->Add(fdEdxBef);
455
456	fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
457	fOutputList->Add(fdEdxAft);
458
459	fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
460	fOutputList->Add(fInvmassLS);
461
462	fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
463	fOutputList->Add(fInvmassULS);
464
465	fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
466	fOutputList->Add(fOpeningAngleLS);
467
468	fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
469	fOutputList->Add(fOpeningAngleULS);
470
471	fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
472	fOutputList->Add(fPhotoElecPt);
473
474	fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
475	fOutputList->Add(fSemiInclElecPt);
476
477	fCent = new TH1F("fCent","Centrality",100,0,100) ;
478	fOutputList->Add(fCent);
479
480	fTPCsubEPres = new TH2F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1,90,0,90);
481	fOutputList->Add(fTPCsubEPres);
482
483	Int_t binsv1[4]={100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, cent
484	Double_t xminv1[4]={-1,-1,-1,0};
485	Double_t xmaxv1[4]={1,1,1,90};
486	fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
487	fOutputList->Add(fEPres);
488
489	Int_t binsv2[7]={100,90,100,100,100,100,100}; // fTPCnSigma,cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
490	Double_t xminv2[7]={-3.5,0,0,0,0,0,0};
491	Double_t xmaxv2[7]={3.5,90,50,3,TMath::Pi(),TMath::Pi(),TMath::Pi()};
492	fCorr = new THnSparseD ("fCorr","Correlations",7,binsv2,xminv2,xmaxv2);
493	fOutputList->Add(fCorr);
494
495	Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
496	Double_t xminv3[5]={0,0,-1,-1,-1};
497	Double_t xmaxv3[5]={90,50,1,1,1};
498	feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
499	fOutputList->Add(feV2);
500
501	Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
502	Double_t xminv4[5]={0,0,-1,-1,-1};
503	Double_t xmaxv4[5]={90,50,1,1,1};
504	fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
505	fOutputList->Add(fphoteV2);
506
507	Int_t binsv5[6]={90,100,100,100,100,100}; // cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
508	Double_t xminv5[6]={0,0,0,-1,-1,-1};
509	Double_t xmaxv5[6]={90,50,3,1,1,1};
510	fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",6,binsv5,xminv5,xmaxv5);
511	fOutputList->Add(fChargPartV2);
512
513	Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
514	Double_t xminv6[5]={0,0,-1,-1,-1};
515	Double_t xmaxv6[5]={90,50,1,1,1};
516	feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
517	fOutputList->Add(feTPCV2);
518
519	PostData(1,fOutputList);
520	}
521
522	//________________________________________________________________________
523	void AliAnalysisTaskElecV2::Terminate(Option_t *)
524	{
525	// Info("Terminate");
526	AliAnalysisTaskSE::Terminate();
527	}
528
529	//________________________________________________________________________
530	Bool_t AliAnalysisTaskElecV2::ProcessCutStep(Int_t cutStep, AliVParticle *track)
531	{
532	// Check single track cuts for a given cut step
533	const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
534	if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
535	return kTRUE;
536	}
537	//_________________________________________
538	void AliAnalysisTaskElecV2::SelectPhotonicElectron(Int_t itrack, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
539	{
540	//Identify non-heavy flavour electrons using Invariant mass method
541
542	fTrackCuts->SetAcceptKinkDaughters(kFALSE);
543	fTrackCuts->SetRequireTPCRefit(kTRUE);
544	fTrackCuts->SetEtaRange(-0.7,0.7);
545	fTrackCuts->SetRequireSigmaToVertex(kTRUE);
546	fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
547	fTrackCuts->SetMinNClustersTPC(100);
548
549	const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
550
551	Bool_t flagPhotonicElec = kFALSE;
552
553	for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
554	AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
555	if (!trackAsso) {
556	printf("ERROR: Could not receive track %d\n", jTracks);
557	continue;
558	}
559
560	Double_t dEdxAsso = -999., ptAsso=-999., openingAngle = -999.;
561	Double_t mass=999., width = -999;
562	Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
563
564	dEdxAsso = trackAsso->GetTPCsignal();
565	ptAsso = trackAsso->Pt();
566	Int_t chargeAsso = trackAsso->Charge();
567	Int_t charge = track->Charge();
568
569	if(ptAsso <0.3) continue;
570	if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
571	if(dEdxAsso <70 \|\| dEdxAsso>100) continue; //11a pass1
572
573	Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
574	if(charge>0) fPDGe1 = -11;
575	if(chargeAsso>0) fPDGe2 = -11;
576
577	if(charge == chargeAsso) fFlagLS = kTRUE;
578	if(charge != chargeAsso) fFlagULS = kTRUE;
579
580	AliKFParticle ge1(*track, fPDGe1);
581	AliKFParticle ge2(*trackAsso, fPDGe2);
582	AliKFParticle recg(ge1, ge2);
583
584	if(recg.GetNDF()<1) continue;
585	Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
586	if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
587
588	AliKFVertex primV(*pVtx);
589	primV += recg;
590	recg.SetProductionVertex(primV);
591
592	recg.SetMassConstraint(0,0.0001);
593
594	openingAngle = ge1.GetAngle(ge2);
595	if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
596	if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
597
598	if(openingAngle > fOpeningAngleCut) continue;
599
600	recg.GetMass(mass,width);
601
602	if(fFlagLS) fInvmassLS->Fill(mass);
603	if(fFlagULS) fInvmassULS->Fill(mass);
604
605	if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
606	flagPhotonicElec = kTRUE;
607	}
608
609	}
610	fFlagPhotonicElec = flagPhotonicElec;
611
612	}
613	//_________________________________________
614	Double_t AliAnalysisTaskElecV2::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
615	{
616	//Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
617	Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
618	if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
619	Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
620
621	return cos2DeltaPhi;
622	}
623
624	//_________________________________________
625	Double_t AliAnalysisTaskElecV2::GetDeltaPhi(Double_t phiA,Double_t phiB) const
626	{
627	//Get phi-psi_EP
628	Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
629	if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
630
631	return dPhi;
632	}