[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); 
            
    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)};
    
    if(fTPCnSigma >= -0.5) feTPCV2->Fill(correctedV2);
    
    if(pidpassed==0) continue;
    
    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);
356
357	Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
358	Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
359	TVector2 newQVectorfortrack;
360	newQVectorfortrack.Set(qX,qY);
361	Double_t corrV2TPC = -999.;
362	corrV2TPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
363
364	Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,corrV2TPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
365
366	if(fTPCnSigma >= -0.5) feTPCV2->Fill(correctedV2);
367
368	if(pidpassed==0) continue;
369
370	feV2->Fill(correctedV2);
371
372	fTrkEovPAft->Fill(pt,EovP);
373	fdEdxAft->Fill(p,dEdx);
374
375	Bool_t fFlagPhotonicElec = kFALSE;
376	SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
377
378	if(fFlagPhotonicElec){
379	fphoteV2->Fill(correctedV2);
380	fPhotoElecPt->Fill(pt);
381	}
382
383	if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
384	}
385	PostData(1, fOutputList);
386	}
387	//_________________________________________
388	void AliAnalysisTaskElecV2::UserCreateOutputObjects()
389	{
390	//--------Initialize PID
391	fPID->SetHasMCData(kFALSE);
392	if(!fPID->GetNumberOfPIDdetectors())
393	{
394	fPID->AddDetector("TPC", 0);
395	fPID->AddDetector("EMCAL", 1);
396	}
397
398	fPID->SortDetectors();
399	fPIDqa = new AliHFEpidQAmanager();
400	fPIDqa->Initialize(fPID);
401
402	//--------Initialize correction Framework and Cuts
403	fCFM = new AliCFManager;
404	const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
405	fCFM->SetNStepParticle(kNcutSteps);
406	for(Int_t istep = 0; istep < kNcutSteps; istep++)
407	fCFM->SetParticleCutsList(istep, NULL);
408
409	if(!fCuts){
410	AliWarning("Cuts not available. Default cuts will be used");
411	fCuts = new AliHFEcuts;
412	fCuts->CreateStandardCuts();
413	}
414	fCuts->Initialize(fCFM);
415
416	//---------Output Tlist
417	fOutputList = new TList();
418	fOutputList->SetOwner();
419	fOutputList->Add(fPIDqa->MakeList("PIDQA"));
420
421	fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
422	fOutputList->Add(fNoEvents);
423
424	fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
425	fOutputList->Add(fTrkpt);
426
427	fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
428	fOutputList->Add(fTrackPtBefTrkCuts);
429
430	fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
431	fOutputList->Add(fTrackPtAftTrkCuts);
432
433	fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
434	fOutputList->Add(fTPCnsigma);
435
436	fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
437	fOutputList->Add(fTrkEovPBef);
438
439	fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
440	fOutputList->Add(fTrkEovPAft);
441
442	fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
443	fOutputList->Add(fdEdxBef);
444
445	fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
446	fOutputList->Add(fdEdxAft);
447
448	fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
449	fOutputList->Add(fInvmassLS);
450
451	fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
452	fOutputList->Add(fInvmassULS);
453
454	fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
455	fOutputList->Add(fOpeningAngleLS);
456
457	fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
458	fOutputList->Add(fOpeningAngleULS);
459
460	fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
461	fOutputList->Add(fPhotoElecPt);
462
463	fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
464	fOutputList->Add(fSemiInclElecPt);
465
466	fCent = new TH1F("fCent","Centrality",100,0,100) ;
467	fOutputList->Add(fCent);
468
469	fTPCsubEPres = new TH2F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1,90,0,90);
470	fOutputList->Add(fTPCsubEPres);
471
472	Int_t binsv1[4]={100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, cent
473	Double_t xminv1[4]={-1,-1,-1,0};
474	Double_t xmaxv1[4]={1,1,1,90};
475	fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
476	fOutputList->Add(fEPres);
477
478	Int_t binsv2[7]={100,90,100,100,100,100,100}; // fTPCnSigma,cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
479	Double_t xminv2[7]={-3.5,0,0,0,0,0,0};
480	Double_t xmaxv2[7]={3.5,90,50,3,TMath::Pi(),TMath::Pi(),TMath::Pi()};
481	fCorr = new THnSparseD ("fCorr","Correlations",7,binsv2,xminv2,xmaxv2);
482	fOutputList->Add(fCorr);
483
484	Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
485	Double_t xminv3[5]={0,0,-1,-1,-1};
486	Double_t xmaxv3[5]={90,50,1,1,1};
487	feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
488	fOutputList->Add(feV2);
489
490	Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
491	Double_t xminv4[5]={0,0,-1,-1,-1};
492	Double_t xmaxv4[5]={90,50,1,1,1};
493	fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
494	fOutputList->Add(fphoteV2);
495
496	Int_t binsv5[6]={90,100,100,100,100,100}; // cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
497	Double_t xminv5[6]={0,0,0,-1,-1,-1};
498	Double_t xmaxv5[6]={90,50,3,1,1,1};
499	fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",6,binsv5,xminv5,xmaxv5);
500	fOutputList->Add(fChargPartV2);
501
502	Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
503	Double_t xminv6[5]={0,0,-1,-1,-1};
504	Double_t xmaxv6[5]={90,50,1,1,1};
505	feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
506	fOutputList->Add(feTPCV2);
507
508	PostData(1,fOutputList);
509	}
510
511	//________________________________________________________________________
512	void AliAnalysisTaskElecV2::Terminate(Option_t *)
513	{
514	// Info("Terminate");
515	AliAnalysisTaskSE::Terminate();
516	}
517
518	//________________________________________________________________________
519	Bool_t AliAnalysisTaskElecV2::ProcessCutStep(Int_t cutStep, AliVParticle *track)
520	{
521	// Check single track cuts for a given cut step
522	const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
523	if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
524	return kTRUE;
525	}
526	//_________________________________________
527	void AliAnalysisTaskElecV2::SelectPhotonicElectron(Int_t itrack, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
528	{
529	//Identify non-heavy flavour electrons using Invariant mass method
530
531	fTrackCuts->SetAcceptKinkDaughters(kFALSE);
532	fTrackCuts->SetRequireTPCRefit(kTRUE);
533	fTrackCuts->SetEtaRange(-0.7,0.7);
534	fTrackCuts->SetRequireSigmaToVertex(kTRUE);
535	fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
536	fTrackCuts->SetMinNClustersTPC(100);
537
538	const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
539
540	Bool_t flagPhotonicElec = kFALSE;
541
542	for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
543	AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
544	if (!trackAsso) {
545	printf("ERROR: Could not receive track %d\n", jTracks);
546	continue;
547	}
548
549	Double_t dEdxAsso = -999., ptAsso=-999., openingAngle = -999.;
550	Double_t mass=999., width = -999;
551	Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
552
553	dEdxAsso = trackAsso->GetTPCsignal();
554	ptAsso = trackAsso->Pt();
555	Int_t chargeAsso = trackAsso->Charge();
556	Int_t charge = track->Charge();
557
558	if(ptAsso <0.3) continue;
559	if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
560	if(dEdxAsso <70 \|\| dEdxAsso>100) continue; //11a pass1
561
562	Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
563	if(charge>0) fPDGe1 = -11;
564	if(chargeAsso>0) fPDGe2 = -11;
565
566	if(charge == chargeAsso) fFlagLS = kTRUE;
567	if(charge != chargeAsso) fFlagULS = kTRUE;
568
569	AliKFParticle ge1(*track, fPDGe1);
570	AliKFParticle ge2(*trackAsso, fPDGe2);
571	AliKFParticle recg(ge1, ge2);
572
573	if(recg.GetNDF()<1) continue;
574	Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
575	if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
576
577	AliKFVertex primV(*pVtx);
578	primV += recg;
579	recg.SetProductionVertex(primV);
580
581	recg.SetMassConstraint(0,0.0001);
582
583	openingAngle = ge1.GetAngle(ge2);
584	if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
585	if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
586
587	if(openingAngle > fOpeningAngleCut) continue;
588
589	recg.GetMass(mass,width);
590
591	if(fFlagLS) fInvmassLS->Fill(mass);
592	if(fFlagULS) fInvmassULS->Fill(mass);
593
594	if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
595	flagPhotonicElec = kTRUE;
596	}
597
598	}
599	fFlagPhotonicElec = flagPhotonicElec;
600
601	}
602	//_________________________________________
603	Double_t AliAnalysisTaskElecV2::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
604	{
605	//Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
606	Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
607	if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
608	Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
609
610	return cos2DeltaPhi;
611	}
612
613	//_________________________________________
614	Double_t AliAnalysisTaskElecV2::GetDeltaPhi(Double_t phiA,Double_t phiB) const
615	{
616	//Get phi-psi_EP
617	Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
618	if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
619
620	return dPhi;
621	}