[u/mrichter/AliRoot.git] / PWGHF / hfe / AliAnalysisTaskHFECal.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  with EMCal triggered events
// Author: Shingo Sakai


#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 "AliAnalysisTaskHFECal.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 "AliCentrality.h"

ClassImp(AliAnalysisTaskHFECal)
//________________________________________________________________________
AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(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)
  ,fEMCAccE(0)
  ,fTrkpt(0)
  ,fTrkEovPBef(0)	 
  ,fTrkEovPAft(0)	
  ,fdEdxBef(0)	 
  ,fdEdxAft(0)	 
  ,fIncpT(0)	
  ,fInvmassLS(0)		
  ,fInvmassULS(0)		
  ,fOpeningAngleLS(0)	
  ,fOpeningAngleULS(0)	
  ,fPhotoElecPt(0)
  ,fPhoElecPt(0)
  ,fSameElecPt(0)
  ,fTrackPtBefTrkCuts(0)	 
  ,fTrackPtAftTrkCuts(0)
  ,fTPCnsigma(0)
  ,fCent(0)
  ,fEleInfo(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());
}

//________________________________________________________________________
AliAnalysisTaskHFECal::AliAnalysisTaskHFECal() 
  : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
  ,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)
  ,fEMCAccE(0)
  ,fTrkpt(0)
  ,fTrkEovPBef(0)	 
  ,fTrkEovPAft(0)	 
  ,fdEdxBef(0)	 
  ,fdEdxAft(0)	 
  ,fIncpT(0)	 
  ,fInvmassLS(0)		
  ,fInvmassULS(0)		
  ,fOpeningAngleLS(0)	
  ,fOpeningAngleULS(0)	
  ,fPhotoElecPt(0)
  ,fPhoElecPt(0)
  ,fSameElecPt(0)
  ,fTrackPtBefTrkCuts(0)	 
  ,fTrackPtAftTrkCuts(0)	 	  
  ,fTPCnsigma(0)
  ,fCent(0)
  ,fEleInfo(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());
}
//_________________________________________

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

void AliAnalysisTaskHFECal::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());
  }
 
  fNoEvents->Fill(0);

  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(1);
  
  if(fNOtrks<2) return;
  fNoEvents->Fill(2);
  
  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;
	
 // Calorimeter info.
 
  // make EMCAL array 
  for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
     {
      AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
      if(clust->IsEMCAL())
        {
         double clustE = clust->E();
         float  emcx[3]; // cluster pos
         clust->GetPosition(emcx);
         TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
         double emcphi = clustpos.Phi(); 
         double emceta = clustpos.Eta();
         double calInfo[4];
         calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; 
         if(clustE>1.5)fEMCAccE->Fill(calInfo); 
        }
   }

  // 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;
    if(TMath::Abs(track->Pt()<2.0)) 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 mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
    pt = track->Pt();
    if(pt<2.0)continue;
    
    // Track extrapolation
    
    Int_t charge = track->Charge();
    fTrkpt->Fill(pt);
    mom = track->P();
    phi = track->Phi();
    eta = track->Eta();
    dEdx = track->GetTPCsignal();
    fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
    
        double ncells = -1.0;
        double m20 = -1.0;
        double m02 = -1.0;
        double disp = -1.0;
        double rmatch = -1.0;  
        double nmatch = -1.0;
        double oppstatus = 0.0;
        double samestatus = 0.0;

    Bool_t fFlagPhotonicElec = kFALSE;
    Bool_t fFlagConvinatElec = kFALSE;
    SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec);
    if(fFlagPhotonicElec)oppstatus = 1.0;
    if(fFlagConvinatElec)samestatus = 1.0;

    Int_t clsId = track->GetEMCALcluster();
    if (clsId>0){
      AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
      if(clust && clust->IsEMCAL()){

	        double clustE = clust->E();
                eop = clustE/fabs(mom);
                //double clustT = clust->GetTOF();
                ncells = clust->GetNCells();
                m02 = clust->GetM02();
                m20 = clust->GetM20();
                disp = clust->GetDispersion();
		double delphi = clust->GetTrackDx(); 
		double deleta = clust->GetTrackDz(); 
		rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
		nmatch = clust->GetNTracksMatched();

		  double valdedx[16];
		  valdedx[0] = mom; valdedx[1] = pt; valdedx[2] = dEdx; valdedx[3] = phi; valdedx[4] = eta; valdedx[5] = fTPCnSigma;
		  valdedx[6] = eop; valdedx[7] = rmatch; valdedx[8] = ncells,  valdedx[9] = m02; valdedx[10] = m20; valdedx[11] = disp;
		  valdedx[12] = cent; valdedx[13] = charge; valdedx[14] = oppstatus; valdedx[15] = samestatus;
                  fEleInfo->Fill(valdedx);
                 

      }
    }
        
    fdEdxBef->Fill(mom,dEdx);
    fTPCnsigma->Fill(mom,fTPCnSigma);
    if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);

    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;

    if(pidpassed==0) continue;
    
    fTrkEovPAft->Fill(pt,eop);
    fdEdxAft->Fill(mom,dEdx);
    fIncpT->Fill(cent,pt);    

    
    if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
    if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
 }
 PostData(1, fOutputList);
}
//_________________________________________
void AliAnalysisTaskHFECal::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","",4,-0.5,3.5) ;
  fOutputList->Add(fNoEvents);
  
  Int_t binsE[4] =    {250, 100,  1000, 100};
  Double_t xminE[4] = {1.0,  -1,   0.0,   0}; 
  Double_t xmaxE[4] = {3.5,   1, 100.0, 100}; 
  fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality",4,binsE,xminE,xmaxE);
  fOutputList->Add(fEMCAccE);

  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);
  
  fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",100,0,100,100,0,50);
  fOutputList->Add(fIncpT);


  Int_t nBinspho[3] =  { 100, 100, 500};
  Double_t minpho[3] = {  0.,  0., 0.};   
  Double_t maxpho[3] = {100., 50., 0.5};   

  fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2);", 3, nBinspho,minpho, maxpho);
  fOutputList->Add(fInvmassLS);
  
  fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2);", 3, nBinspho,minpho, maxpho);
  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);
  
  fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",100,0,100,100,0,50);
  fOutputList->Add(fPhoElecPt);
  
  fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",100,0,100,100,0,50);
  fOutputList->Add(fSameElecPt);

  fCent = new TH1F("fCent","Centrality",100,0,100) ;
  fOutputList->Add(fCent);
 
  // Make common binning
  const Double_t kMinP = 2.;
  const Double_t kMaxP = 50.;
  const Double_t kTPCSigMim = 40.;
  const Double_t kTPCSigMax = 140.;

  // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
  Int_t nBins[16] =  {  480,   480,        200,   60,    20,   600,  300, 100,   40,   200, 200, 200, 100,    3,    3,    3};
  Double_t min[16] = {kMinP, kMinP,  kTPCSigMim, 1.0,  -1.0,  -8.0,    0,   0,    0,   0.0, 0.0, 0.0,   0, -1.5, -0.5, -0.5};
  Double_t max[16] = {kMaxP, kMaxP,  kTPCSigMax, 4.0,   1.0,   4.0,  3.0, 0.1,   40,   2.0, 2.0, 2.0, 100,  1.5,  2.5,  2.5};
  fEleInfo = new THnSparseD("fEleInfo", "Electron Info; p [GeV/c]; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;M02;M20;Disp; Centrality; charge", 16, nBins, min, max);
  fOutputList->Add(fEleInfo);

//_________________________________________________________
 
  PostData(1,fOutputList);
}

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

//________________________________________________________________________
Bool_t AliAnalysisTaskHFECal::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 AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec)
{
  //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(70);
  
  const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
  
  Bool_t flagPhotonicElec = kFALSE;
  Bool_t flagConvinatElec = kFALSE;
  
  //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
  for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
    AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
    if (!trackAsso) {
      printf("ERROR: Could not receive track %d\n", jTracks);
      continue;
    }
    if(itrack==jTracks)continue;    

    Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
    Double_t mass=999., width = -999;
    Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
    
    ptPrim = track->Pt();

    dEdxAsso = trackAsso->GetTPCsignal();
    ptAsso = trackAsso->Pt();
    Int_t chargeAsso = trackAsso->Charge();
    Int_t charge = track->Charge();
    
    //if(ptAsso <0.3) continue;
    if(ptAsso <0.5) continue;
    if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
    if(dEdxAsso <65 || 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);
    
    double phoinfo[3];
    phoinfo[0] = cent;
    phoinfo[1] = ptPrim;
    phoinfo[2] = mass;

    if(fFlagLS) fInvmassLS->Fill(phoinfo);
    if(fFlagULS) fInvmassULS->Fill(phoinfo);
          
    if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
      flagPhotonicElec = kTRUE;
    }
    if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
      flagConvinatElec = kTRUE;
    }
    
  }
  fFlagPhotonicElec = flagPhotonicElec;
  fFlagConvinatElec = flagConvinatElec;
  
}
Commit	Line	Data
bfc7c23b	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 with EMCal triggered events
	17	// Author: Shingo Sakai
	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 "AliAnalysisTaskHFECal.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 "AliCentrality.h"
76
77	ClassImp(AliAnalysisTaskHFECal)
78	//________________________________________________________________________
79	AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
80	: AliAnalysisTaskSE(name)
81	,fESD(0)
82	,fOutputList(0)
83	,fTrackCuts(0)
84	,fCuts(0)
85	,fIdentifiedAsOutInz(kFALSE)
86	,fPassTheEventCut(kFALSE)
87	,fRejectKinkMother(kFALSE)
88	,fVz(0.0)
89	,fCFM(0)
90	,fPID(0)
91	,fPIDqa(0)
92	,fOpeningAngleCut(0.1)
93	,fInvmassCut(0.01)
94	,fNoEvents(0)
95	,fEMCAccE(0)
96	,fTrkpt(0)
97	,fTrkEovPBef(0)
98	,fTrkEovPAft(0)
99	,fdEdxBef(0)
100	,fdEdxAft(0)
101	,fIncpT(0)
102	,fInvmassLS(0)
103	,fInvmassULS(0)
104	,fOpeningAngleLS(0)
105	,fOpeningAngleULS(0)
106	,fPhotoElecPt(0)
107	,fPhoElecPt(0)
f09766dd	108	,fSameElecPt(0)
bfc7c23b	109	,fTrackPtBefTrkCuts(0)
	110	,fTrackPtAftTrkCuts(0)
	111	,fTPCnsigma(0)
	112	,fCent(0)
	113	,fEleInfo(0)
	114	{
	115	//Named constructor
	116
	117	fPID = new AliHFEpid("hfePid");
	118	fTrackCuts = new AliESDtrackCuts();
	119
	120	// Define input and output slots here
	121	// Input slot #0 works with a TChain
	122	DefineInput(0, TChain::Class());
	123	// Output slot #0 id reserved by the base class for AOD
	124	// Output slot #1 writes into a TH1 container
	125	// DefineOutput(1, TH1I::Class());
	126	DefineOutput(1, TList::Class());
	127	// DefineOutput(3, TTree::Class());
	128	}
	129
	130	//________________________________________________________________________
	131	AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
	132	: AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
	133	,fESD(0)
	134	,fOutputList(0)
	135	,fTrackCuts(0)
	136	,fCuts(0)
	137	,fIdentifiedAsOutInz(kFALSE)
	138	,fPassTheEventCut(kFALSE)
	139	,fRejectKinkMother(kFALSE)
	140	,fVz(0.0)
	141	,fCFM(0)
	142	,fPID(0)
	143	,fPIDqa(0)
	144	,fOpeningAngleCut(0.1)
	145	,fInvmassCut(0.01)
	146	,fNoEvents(0)
	147	,fEMCAccE(0)
	148	,fTrkpt(0)
	149	,fTrkEovPBef(0)
	150	,fTrkEovPAft(0)
	151	,fdEdxBef(0)
	152	,fdEdxAft(0)
	153	,fIncpT(0)
	154	,fInvmassLS(0)
	155	,fInvmassULS(0)
	156	,fOpeningAngleLS(0)
	157	,fOpeningAngleULS(0)
	158	,fPhotoElecPt(0)
	159	,fPhoElecPt(0)
f09766dd	160	,fSameElecPt(0)
bfc7c23b	161	,fTrackPtBefTrkCuts(0)
	162	,fTrackPtAftTrkCuts(0)
	163	,fTPCnsigma(0)
	164	,fCent(0)
	165	,fEleInfo(0)
	166	{
	167	//Default constructor
	168	fPID = new AliHFEpid("hfePid");
	169
	170	fTrackCuts = new AliESDtrackCuts();
	171
	172	// Constructor
	173	// Define input and output slots here
	174	// Input slot #0 works with a TChain
	175	DefineInput(0, TChain::Class());
	176	// Output slot #0 id reserved by the base class for AOD
	177	// Output slot #1 writes into a TH1 container
	178	// DefineOutput(1, TH1I::Class());
	179	DefineOutput(1, TList::Class());
	180	//DefineOutput(3, TTree::Class());
	181	}
	182	//_________________________________________
	183
	184	AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
	185	{
	186	//Destructor
	187
	188	delete fOutputList;
	189	delete fPID;
	190	delete fCFM;
	191	delete fPIDqa;
	192	delete fTrackCuts;
	193	}
	194	//_________________________________________
	195
	196	void AliAnalysisTaskHFECal::UserExec(Option_t*)
	197	{
	198	//Main loop
	199	//Called for each event
	200
	201	// create pointer to event
	202	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
	203	if (!fESD) {
	204	printf("ERROR: fESD not available\n");
	205	return;
	206	}
	207
	208	if(!fCuts){
	209	AliError("HFE cuts not available");
	210	return;
	211	}
	212
	213	if(!fPID->IsInitialized()){
	214	// Initialize PID with the given run number
	215	AliWarning("PID not initialised, get from Run no");
	216	fPID->InitializePID(fESD->GetRunNumber());
	217	}
	218
	219	fNoEvents->Fill(0);
	220
	221	Int_t fNOtrks = fESD->GetNumberOfTracks();
	222	const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
	223
	224	Double_t pVtxZ = -999;
225	pVtxZ = pVtx->GetZ();
226
227	if(TMath::Abs(pVtxZ)>10) return;
228	fNoEvents->Fill(1);
229
230	if(fNOtrks<2) return;
231	fNoEvents->Fill(2);
232
233	AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
234	if(!pidResponse){
235	AliDebug(1, "Using default PID Response");
236	pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
237	}
238
239	fPID->SetPIDResponse(pidResponse);
240
241	fCFM->SetRecEventInfo(fESD);
242
243	Float_t cent = -1.;
244	AliCentrality *centrality = fESD->GetCentrality();
245	cent = centrality->GetCentralityPercentile("V0M");
246	fCent->Fill(cent);
247
248	if(cent>90.) return;
249
250	// Calorimeter info.
251
252	// make EMCAL array
253	for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
254	{
255	AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
256	if(clust->IsEMCAL())
257	{
258	double clustE = clust->E();
259	float emcx[3]; // cluster pos
260	clust->GetPosition(emcx);
261	TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
262	double emcphi = clustpos.Phi();
263	double emceta = clustpos.Eta();
5e65985c	264	double calInfo[4];
bfc7c23b	265	calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent;
f09766dd	266	if(clustE>1.5)fEMCAccE->Fill(calInfo);
bfc7c23b	267	}
	268	}
	269
	270	// Track loop
	271	for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
	272	AliESDtrack* track = fESD->GetTrack(iTracks);
	273	if (!track) {
	274	printf("ERROR: Could not receive track %d\n", iTracks);
	275	continue;
	276	}
	277
	278	if(TMath::Abs(track->Eta())>0.7) continue;
f09766dd	279	if(TMath::Abs(track->Pt()<2.0)) continue;
bfc7c23b	280
	281	fTrackPtBefTrkCuts->Fill(track->Pt());
	282	// RecKine: ITSTPC cuts
	283	if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
	284
	285	//RecKink
	286	if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
	287	if(track->GetKinkIndex(0) != 0) continue;
	288	}
	289
	290	// RecPrim
	291	if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
	292
	293	// HFEcuts: ITS layers cuts
	294	if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
	295
f09766dd	296	// HFE cuts: TPC PID cleanup
	297	if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
	298
bfc7c23b	299
	300	fTrackPtAftTrkCuts->Fill(track->Pt());
	301
	302	Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
f09766dd	303	pt = track->Pt();
f09766dd	304	if(pt<2.0)continue;
bfc7c23b	305
	306	// Track extrapolation
	307
bfc7c23b	308	Int_t charge = track->Charge();
	309	fTrkpt->Fill(pt);
	310	mom = track->P();
	311	phi = track->Phi();
	312	eta = track->Eta();
	313	dEdx = track->GetTPCsignal();
	314	fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
	315
	316	double ncells = -1.0;
	317	double m20 = -1.0;
	318	double m02 = -1.0;
	319	double disp = -1.0;
	320	double rmatch = -1.0;
	321	double nmatch = -1.0;
f09766dd	322	double oppstatus = 0.0;
f09766dd	323	double samestatus = 0.0;
bfc7c23b	324
f09766dd	325	Bool_t fFlagPhotonicElec = kFALSE;
	326	Bool_t fFlagConvinatElec = kFALSE;
	327	SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec);
	328	if(fFlagPhotonicElec)oppstatus = 1.0;
	329	if(fFlagConvinatElec)samestatus = 1.0;
bfc7c23b	330
	331	Int_t clsId = track->GetEMCALcluster();
	332	if (clsId>0){
	333	AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
	334	if(clust && clust->IsEMCAL()){
	335
	336	double clustE = clust->E();
	337	eop = clustE/fabs(mom);
	338	//double clustT = clust->GetTOF();
	339	ncells = clust->GetNCells();
	340	m02 = clust->GetM02();
	341	m20 = clust->GetM20();
	342	disp = clust->GetDispersion();
	343	double delphi = clust->GetTrackDx();
	344	double deleta = clust->GetTrackDz();
	345	rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
	346	nmatch = clust->GetNTracksMatched();
	347
f09766dd	348	double valdedx[16];
bfc7c23b	349	valdedx[0] = mom; valdedx[1] = pt; valdedx[2] = dEdx; valdedx[3] = phi; valdedx[4] = eta; valdedx[5] = fTPCnSigma;
bfc7c23b	350	valdedx[6] = eop; valdedx[7] = rmatch; valdedx[8] = ncells, valdedx[9] = m02; valdedx[10] = m20; valdedx[11] = disp;
f09766dd	351	valdedx[12] = cent; valdedx[13] = charge; valdedx[14] = oppstatus; valdedx[15] = samestatus;
bfc7c23b	352	fEleInfo->Fill(valdedx);
	353
	354
	355	}
	356	}
	357
	358	fdEdxBef->Fill(mom,dEdx);
	359	fTPCnsigma->Fill(mom,fTPCnSigma);
f09766dd	360	if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
bfc7c23b	361
bfc7c23b	362	Int_t pidpassed = 1;
	363
	364
	365	//--- track accepted
	366	AliHFEpidObject hfetrack;
	367	hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
	368	hfetrack.SetRecTrack(track);
	369	hfetrack.SetPbPb();
	370	if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
	371
	372	if(pidpassed==0) continue;
	373
	374	fTrkEovPAft->Fill(pt,eop);
	375	fdEdxAft->Fill(mom,dEdx);
	376	fIncpT->Fill(cent,pt);
	377
bfc7c23b	378
bfc7c23b	379	if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
f09766dd	380	if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
bfc7c23b	381	}
	382	PostData(1, fOutputList);
	383	}
	384	//_________________________________________
	385	void AliAnalysisTaskHFECal::UserCreateOutputObjects()
	386	{
	387	//--------Initialize PID
	388	fPID->SetHasMCData(kFALSE);
	389	if(!fPID->GetNumberOfPIDdetectors())
	390	{
	391	fPID->AddDetector("TPC", 0);
	392	fPID->AddDetector("EMCAL", 1);
	393	}
	394
	395	fPID->SortDetectors();
	396	fPIDqa = new AliHFEpidQAmanager();
	397	fPIDqa->Initialize(fPID);
	398
	399	//--------Initialize correction Framework and Cuts
	400	fCFM = new AliCFManager;
	401	const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
	402	fCFM->SetNStepParticle(kNcutSteps);
	403	for(Int_t istep = 0; istep < kNcutSteps; istep++)
	404	fCFM->SetParticleCutsList(istep, NULL);
	405
	406	if(!fCuts){
	407	AliWarning("Cuts not available. Default cuts will be used");
	408	fCuts = new AliHFEcuts;
	409	fCuts->CreateStandardCuts();
	410	}
	411	fCuts->Initialize(fCFM);
	412
	413	//---------Output Tlist
	414	fOutputList = new TList();
	415	fOutputList->SetOwner();
	416	fOutputList->Add(fPIDqa->MakeList("PIDQA"));
	417
	418	fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
	419	fOutputList->Add(fNoEvents);
	420
	421	Int_t binsE[4] = {250, 100, 1000, 100};
	422	Double_t xminE[4] = {1.0, -1, 0.0, 0};
	423	Double_t xmaxE[4] = {3.5, 1, 100.0, 100};
	424	fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality",4,binsE,xminE,xmaxE);
	425	fOutputList->Add(fEMCAccE);
	426
	427	fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
	428	fOutputList->Add(fTrkpt);
	429
	430	fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
	431	fOutputList->Add(fTrackPtBefTrkCuts);
	432
	433	fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
	434	fOutputList->Add(fTrackPtAftTrkCuts);
	435
	436	fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
	437	fOutputList->Add(fTPCnsigma);
	438
	439	fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
	440	fOutputList->Add(fTrkEovPBef);
	441
	442	fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
	443	fOutputList->Add(fTrkEovPAft);
	444
445	fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
446	fOutputList->Add(fdEdxBef);
447
448	fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
449	fOutputList->Add(fdEdxAft);
450
451	fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",100,0,100,100,0,50);
452	fOutputList->Add(fIncpT);
453
f09766dd	454
	455	Int_t nBinspho[3] = { 100, 100, 500};
	456	Double_t minpho[3] = { 0., 0., 0.};
	457	Double_t maxpho[3] = {100., 50., 0.5};
	458
	459	fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2);", 3, nBinspho,minpho, maxpho);
bfc7c23b	460	fOutputList->Add(fInvmassLS);
bfc7c23b	461
f09766dd	462	fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2);", 3, nBinspho,minpho, maxpho);
bfc7c23b	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
f09766dd	474	fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",100,0,100,100,0,50);
bfc7c23b	475	fOutputList->Add(fPhoElecPt);
bfc7c23b	476
f09766dd	477	fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",100,0,100,100,0,50);
	478	fOutputList->Add(fSameElecPt);
	479
bfc7c23b	480	fCent = new TH1F("fCent","Centrality",100,0,100) ;
	481	fOutputList->Add(fCent);
	482
	483	// Make common binning
f09766dd	484	const Double_t kMinP = 2.;
bfc7c23b	485	const Double_t kMaxP = 50.;
	486	const Double_t kTPCSigMim = 40.;
	487	const Double_t kTPCSigMax = 140.;
	488
	489	// 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select)
f09766dd	490	Int_t nBins[16] = { 480, 480, 200, 60, 20, 600, 300, 100, 40, 200, 200, 200, 100, 3, 3, 3};
	491	Double_t min[16] = {kMinP, kMinP, kTPCSigMim, 1.0, -1.0, -8.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, -1.5, -0.5, -0.5};
	492	Double_t max[16] = {kMaxP, kMaxP, kTPCSigMax, 4.0, 1.0, 4.0, 3.0, 0.1, 40, 2.0, 2.0, 2.0, 100, 1.5, 2.5, 2.5};
	493	fEleInfo = new THnSparseD("fEleInfo", "Electron Info; p [GeV/c]; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;M02;M20;Disp; Centrality; charge", 16, nBins, min, max);
bfc7c23b	494	fOutputList->Add(fEleInfo);
	495
	496	//_________________________________________________________
	497
	498	PostData(1,fOutputList);
	499	}
	500
	501	//________________________________________________________________________
	502	void AliAnalysisTaskHFECal::Terminate(Option_t *)
	503	{
	504	// Info("Terminate");
	505	AliAnalysisTaskSE::Terminate();
	506	}
	507
	508	//________________________________________________________________________
	509	Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
	510	{
	511	// Check single track cuts for a given cut step
	512	const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
	513	if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
	514	return kTRUE;
	515	}
	516	//_________________________________________
f09766dd	517	//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
f09766dd	518	void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec)
bfc7c23b	519	{
	520	//Identify non-heavy flavour electrons using Invariant mass method
	521
f09766dd	522	//fTrackCuts->SetAcceptKinkDaughters(kFALSE);
	523	//fTrackCuts->SetRequireTPCRefit(kTRUE);
	524	//fTrackCuts->SetEtaRange(-0.7,0.7);
	525	//fTrackCuts->SetRequireSigmaToVertex(kTRUE);
bfc7c23b	526	fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
f09766dd	527	fTrackCuts->SetMinNClustersTPC(70);
bfc7c23b	528
	529	const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
	530
	531	Bool_t flagPhotonicElec = kFALSE;
f09766dd	532	Bool_t flagConvinatElec = kFALSE;
bfc7c23b	533
5e65985c	534	//for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
5e65985c	535	for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
bfc7c23b	536	AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
	537	if (!trackAsso) {
	538	printf("ERROR: Could not receive track %d\n", jTracks);
	539	continue;
	540	}
5e65985c	541	if(itrack==jTracks)continue;
5e65985c	542
f09766dd	543	Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
bfc7c23b	544	Double_t mass=999., width = -999;
	545	Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
	546
f09766dd	547	ptPrim = track->Pt();
f09766dd	548
bfc7c23b	549	dEdxAsso = trackAsso->GetTPCsignal();
	550	ptAsso = trackAsso->Pt();
	551	Int_t chargeAsso = trackAsso->Charge();
	552	Int_t charge = track->Charge();
	553
5e65985c	554	//if(ptAsso <0.3) continue;
5e65985c	555	if(ptAsso <0.5) continue;
bfc7c23b	556	if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
f09766dd	557	if(dEdxAsso <65 \|\| dEdxAsso>100) continue; //11a pass1
bfc7c23b	558
	559	Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
	560	if(charge>0) fPDGe1 = -11;
	561	if(chargeAsso>0) fPDGe2 = -11;
	562
	563	if(charge == chargeAsso) fFlagLS = kTRUE;
	564	if(charge != chargeAsso) fFlagULS = kTRUE;
	565
	566	AliKFParticle ge1(*track, fPDGe1);
	567	AliKFParticle ge2(*trackAsso, fPDGe2);
	568	AliKFParticle recg(ge1, ge2);
	569
	570	if(recg.GetNDF()<1) continue;
	571	Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
	572	if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
	573
	574	AliKFVertex primV(*pVtx);
	575	primV += recg;
	576	recg.SetProductionVertex(primV);
	577
	578	recg.SetMassConstraint(0,0.0001);
	579
	580	openingAngle = ge1.GetAngle(ge2);
	581	if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
	582	if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
	583
	584	if(openingAngle > fOpeningAngleCut) continue;
	585
	586	recg.GetMass(mass,width);
	587
f09766dd	588	double phoinfo[3];
	589	phoinfo[0] = cent;
	590	phoinfo[1] = ptPrim;
	591	phoinfo[2] = mass;
	592
	593	if(fFlagLS) fInvmassLS->Fill(phoinfo);
	594	if(fFlagULS) fInvmassULS->Fill(phoinfo);
bfc7c23b	595
	596	if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
	597	flagPhotonicElec = kTRUE;
	598	}
f09766dd	599	if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
	600	flagConvinatElec = kTRUE;
	601	}
bfc7c23b	602
	603	}
	604	fFlagPhotonicElec = flagPhotonicElec;
f09766dd	605	fFlagConvinatElec = flagConvinatElec;
bfc7c23b	606
	607	}
	608