[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)
  ,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)
  ,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[3];
         calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; 
         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;
    
    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;
    
    
    fTrackPtAftTrkCuts->Fill(track->Pt());		
    
    Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
    
    // Track extrapolation
    
    pt = track->Pt();
    if(pt<2.0)continue;
    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;


    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[14];
		  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;
                  fEleInfo->Fill(valdedx);
                 

      }
    }
        
    fdEdxBef->Fill(mom,dEdx);
    fTPCnsigma->Fill(mom,fTPCnSigma);
 

    // HFE cuts: TPC PID cleanup
    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
      
    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);    

    Bool_t fFlagPhotonicElec = kFALSE;
    SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec);
    
    if(fFlagPhotonicElec) fPhoElecPt->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);

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