[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALIsoPhoton.cxx

// $Id$

#include "AliAnalysisTaskEMCALIsoPhoton.h"

#include <TCanvas.h>
#include <TChain.h>
#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <THnSparse.h>
#include <TLorentzVector.h>
#include <TList.h>

#include "AliAnalysisManager.h"
#include "AliAnalysisTask.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALRecoUtils.h"
#include "AliESDCaloCells.h"
#include "AliESDCaloCluster.h"
#include "AliESDEvent.h"
#include "AliESDHeader.h"
#include "AliESDInputHandler.h"
#include "AliESDUtils.h"
#include "AliESDpid.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
#include "AliESDv0.h"
#include "AliKFParticle.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliStack.h"
#include "AliV0vertexer.h"
#include "AliVCluster.h"

ClassImp(AliAnalysisTaskEMCALIsoPhoton)

//________________________________________________________________________
AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() : 
  AliAnalysisTaskSE(), 
  fCaloClusters(0),
  fSelPrimTracks(0),
  fTracks(0),
  fEMCalCells(0),
  fPrTrCuts(0),
  fGeom(0x0),
  fGeoName("EMCAL_COMPLETEV1"),
  fPeriod("LHC11c"),
  fTrigBit("kEMC7"),
  fIsTrain(0),
  fExoticCut(0.97),
  fIsoConeR(0.4),
  fNDimensions(7),
  fECut(3.),
  fESD(0),
  fOutputList(0),
  fEvtSel(0),
  fPVtxZ(0),                 
  fCellAbsIdVsAmpl(0),       
  fNClusHighClusE(0),        
  fHnOutput(0)
{
  // Default constructor.
}

//________________________________________________________________________
AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) : 
  AliAnalysisTaskSE(name), 
  fCaloClusters(0),
  fSelPrimTracks(0),
  fTracks(0),
  fEMCalCells(0),
  fPrTrCuts(0),
  fGeom(0x0),
  fGeoName("EMCAL_COMPLETEV1"),
  fPeriod("LHC11c"),
  fTrigBit("kEMC7"),
  fIsTrain(0),
  fExoticCut(0.97),
  fIsoConeR(0.4),
  fNDimensions(7),
  fECut(3.),
  fESD(0),
  fOutputList(0),
  fEvtSel(0),
  fPVtxZ(0),            
  fCellAbsIdVsAmpl(0),  
  fNClusHighClusE(0),   
  fHnOutput(0)
{
  // 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, TList::Class());
}

//________________________________________________________________________
void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
{
  // Create histograms, called once.
    
  fCaloClusters = new TRefArray();
  fSelPrimTracks = new TObjArray();

  
  fOutputList = new TList();
  fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging) 
  
  fGeom = AliEMCALGeometry::GetInstance(fGeoName);
  
  fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
  fOutputList->Add(fEvtSel);
    
  
  fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
  fOutputList->Add(fPVtxZ);
  
  fCellAbsIdVsAmpl = new TH2F("hCellAbsIdVsAmpl","cell abs id vs cell amplitude (energy);E (GeV);ID",200,0,100,24*48*10,-0.5,24*48*10-0.5);
  fOutputList->Add(fCellAbsIdVsAmpl);

  fNClusHighClusE = new TH2F("hNClusHighClusE","total number of clusters vs. highest clus energy in the event;E (GeV);NClus",200,0,100,301,-0.5,300.5);
  fOutputList->Add(fNClusHighClusE);

  const Int_t ndims =   fNDimensions;
  Int_t nEt=1000, nM02=400, nCeIso=1000, nTrIso=1000,  nAllIso=1000, nTrClDphi=200, nTrClDeta=100;
  Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nTrClDphi, nTrClDeta};
  Double_t xmin[] = { 0.,   0.,  -10.,   -10., -10., -0.1,-0.05};
  Double_t xmax[] = { 100., 4., 190., 190., 190., 0.1, 0.05};
  fHnOutput =  new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, d#phi_{trk},d#eta_{trk}", ndims, bins, xmin, xmax);
  fOutputList->Add(fHnOutput);


  PostData(1, fOutputList);
}

//________________________________________________________________________
void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) 
{
  // Main loop, called for each event.

  // event trigger selection
  Bool_t isSelected = 0;
  if(fPeriod.Contains("11a")){
    if(fTrigBit.Contains("kEMC"))
      isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
    else
      isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
  }
  else{
    if(fTrigBit.Contains("kEMC"))
      isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
    else
      isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
  }
  if(!isSelected )
        return; 

  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!fESD) {
    printf("ERROR: fESD not available\n");
    return;
  }
  
  fEvtSel->Fill(0);
  
  AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
  if(!pv) 
    return;
  fPVtxZ->Fill(pv->GetZ());
  if(TMath::Abs(pv->GetZ())>15)
    return;

  fEvtSel->Fill(1);

  if (!fTracks)  
    fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
  // Track loop to fill a pT spectrum
  const Int_t Ntracks = fTracks->GetEntriesFast();
  for (Int_t iTracks = 0;  iTracks < Ntracks; ++iTracks) {
    //  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
    //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
    AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(iTracks));
    if (!track)
      continue;
    if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
      fSelPrimTracks->Add(track);
      //printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi());
    }
  }

  if(!fIsTrain){
    for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
      if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
        break;
      fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
    }
  }
  fESD->GetEMCALClusters(fCaloClusters);
  fEMCalCells = fESD->GetEMCALCells();
  
    
  FillClusHists(); 

  fCaloClusters->Clear();
  fSelPrimTracks->Clear();
  PostData(1, fOutputList);
}      

//________________________________________________________________________
void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
{
  // Fill cluster histograms.

  if(!fCaloClusters)
    return;
  const Int_t nclus = fCaloClusters->GetEntries();
  if(nclus==0)
    return;
  Double_t maxE = 0;
  for(Int_t ic=0;ic<nclus;ic++){
    maxE=0;
    AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
    if(!c)
      continue;
    if(!c->IsEMCAL())
      continue;
    if(c->E()<fECut)
      continue;
    Short_t id;
    Double_t Emax = GetMaxCellEnergy( c, id);
    Double_t Ecross = GetCrossEnergy( c, id);
    if((1-Ecross/Emax)>fExoticCut)
      continue;
    Float_t clsPos[3] = {0,0,0};
    c->GetPosition(clsPos);
    TVector3 clsVec(clsPos);
    Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
    Float_t ceiso, cephiband, cecore;
    Float_t triso, trphiband, trcore;
    Float_t alliso, allphiband, allcore;
    Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR;
    Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04);
    GetCeIso(clsVec, ceiso, cephiband, cecore);
    GetTrIso(clsVec, triso, trphiband, trcore);
    alliso = ceiso + triso;
    allphiband = cephiband + trphiband;
    allcore = cecore + trcore;
    Float_t ceisoue =  cephiband/phibandArea*netConeArea;
    Float_t trisoue =  trphiband/phibandArea*netConeArea;
    Float_t allisoue =  allphiband/phibandArea*netConeArea;
    const Int_t ndims =   fNDimensions;
    Double_t outputValues[ndims];
    outputValues[0] = Et;
    outputValues[1] = c->GetM02();
    outputValues[2] = ceiso-cecore-ceisoue;
    outputValues[3] = triso-trisoue;
    outputValues[4] = alliso-cecore-allisoue;
    outputValues[5] = c->GetTrackDx();
    outputValues[6] = c->GetTrackDz();
    fHnOutput->Fill(outputValues);
    if(c->E()>maxE)
      maxE = c->E();
  }
  fNClusHighClusE->Fill(maxE,nclus);
} 

//________________________________________________________________________
void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
{
  // Get cell isolation.

  if(!fEMCalCells)
    return;
  const Int_t ncells = fEMCalCells->GetNumberOfCells();
  Float_t totiso=0;
  Float_t totband=0;
  Float_t totcore=0;
  Float_t etacl = vec.Eta();
  Float_t phicl = vec.Phi();
  Float_t thetacl = vec.Theta();
  if(phicl<0)
    phicl+=TMath::TwoPi();
  Int_t absid = -1;
  Float_t eta=-1, phi=-1;  
  for(int icell=0;icell<ncells;icell++){
    absid = TMath::Abs(fEMCalCells->GetCellNumber(icell));
    if(!fGeom)
      return;
    fGeom->EtaPhiFromIndex(absid,eta,phi);
    Float_t dphi = TMath::Abs(phi-phicl);
    Float_t deta = TMath::Abs(eta-etacl);
    Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
    Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
    if(R<fIsoConeR){
      totiso += etcell;
      if(R<0.04)
	totcore += etcell;
    }
    else{
      if(dphi>fIsoConeR)
	continue;
      if(deta<fIsoConeR)
	continue;
      totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
    }
  }
  iso = totiso;
  phiband = totband;
  core = totcore;
} 

//________________________________________________________________________
void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
{
  // Get track isolation.

  if(!fSelPrimTracks)
    return;
  const Int_t ntracks = fSelPrimTracks->GetEntries();
  Double_t totiso=0;
  Double_t totband=0;
  Double_t totcore=0;
  Double_t etacl = vec.Eta();
  Double_t phicl = vec.Phi();
  if(phicl<0)
    phicl+=TMath::TwoPi();
  for(int itrack=0;itrack<ntracks;itrack++){
    AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack));
    if(!track)
      continue;
    Double_t dphi = TMath::Abs(track->Phi()-phicl);
    Double_t deta = TMath::Abs(track->Eta()-etacl);
    Double_t R = TMath::Sqrt(deta*deta + dphi*dphi);
    Double_t pt = track->Pt();
    if(R<fIsoConeR){
      totiso += track->Pt();
      if(R<0.04)
	totcore += pt;
    }
    else{
      if(dphi>fIsoConeR)
	continue;
      if(deta<fIsoConeR)
	continue;
      totband += track->Pt();
    }
  }
  iso = totiso;
  phiband = totband;
  core = totcore;
} 

//________________________________________________________________________
Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
{
  // Calculate the energy of cross cells around the leading cell.

  AliVCaloCells *cells = 0;
  cells = fESD->GetEMCALCells();
  if (!cells)
    return 0;

  if (!fGeom)
    return 0;

  Int_t iSupMod = -1;
  Int_t iTower  = -1;
  Int_t iIphi   = -1;
  Int_t iIeta   = -1;
  Int_t iphi    = -1;
  Int_t ieta    = -1;
  Int_t iphis   = -1;
  Int_t ietas   = -1;

  Double_t crossEnergy = 0;

  fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
  fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);

  Int_t ncells = cluster->GetNCells();
  for (Int_t i=0; i<ncells; i++) {
    Int_t cellAbsId = cluster->GetCellAbsId(i);
    fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
    fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
    Int_t aphidiff = TMath::Abs(iphi-iphis);
    if (aphidiff>1)
      continue;
    Int_t aetadiff = TMath::Abs(ieta-ietas);
    if (aetadiff>1)
      continue;
    if ( (aphidiff==1 && aetadiff==0) ||
	(aphidiff==0 && aetadiff==1) ) {
      crossEnergy += cells->GetCellAmplitude(cellAbsId);
    }
  }

  return crossEnergy;
}

//________________________________________________________________________
Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
{
  // Get maximum energy of attached cell.

  id = -1;

  AliVCaloCells *cells = 0;
  cells = fESD->GetEMCALCells();
  if (!cells)
    return 0;

  Double_t maxe = 0;
  Int_t ncells = cluster->GetNCells();
  for (Int_t i=0; i<ncells; i++) {
    Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
    if (e>maxe) {
      maxe = e;
      id   = cluster->GetCellAbsId(i);
    }
  }
  return maxe;
}

//________________________________________________________________________
void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *) 
{
  // Called once at the end of the query.
}
Commit	Line	Data
bd092f0f	1	// $Id$
30159e6f	2
bd092f0f	3	#include "AliAnalysisTaskEMCALIsoPhoton.h"
30159e6f	4
bd092f0f	5	#include <TCanvas.h>
	6	#include <TChain.h>
	7	#include <TFile.h>
	8	#include <TH1F.h>
	9	#include <TH2F.h>
965c985f	10	#include <THnSparse.h>
bd092f0f	11	#include <TLorentzVector.h>
a62631a9	12	#include <TList.h>
bd092f0f	13
	14	#include "AliAnalysisManager.h"
	15	#include "AliAnalysisTask.h"
	16	#include "AliEMCALGeometry.h"
	17	#include "AliEMCALRecoUtils.h"
	18	#include "AliESDCaloCells.h"
	19	#include "AliESDCaloCluster.h"
30159e6f	20	#include "AliESDEvent.h"
30159e6f	21	#include "AliESDHeader.h"
30159e6f	22	#include "AliESDInputHandler.h"
bd092f0f	23	#include "AliESDUtils.h"
30159e6f	24	#include "AliESDpid.h"
30159e6f	25	#include "AliESDtrack.h"
	26	#include "AliESDtrackCuts.h"
	27	#include "AliESDv0.h"
bd092f0f	28	#include "AliKFParticle.h"
	29	#include "AliMCEvent.h"
	30	#include "AliMCEventHandler.h"
	31	#include "AliStack.h"
30159e6f	32	#include "AliV0vertexer.h"
30159e6f	33	#include "AliVCluster.h"
30159e6f	34
30159e6f	35	ClassImp(AliAnalysisTaskEMCALIsoPhoton)
	36
	37	//________________________________________________________________________
bd092f0f	38	AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() :
	39	AliAnalysisTaskSE(),
	40	fCaloClusters(0),
	41	fSelPrimTracks(0),
3f4073ba	42	fTracks(0),
bd092f0f	43	fEMCalCells(0),
	44	fPrTrCuts(0),
	45	fGeom(0x0),
9148fa89	46	fGeoName("EMCAL_COMPLETEV1"),
	47	fPeriod("LHC11c"),
	48	fTrigBit("kEMC7"),
	49	fIsTrain(0),
	50	fExoticCut(0.97),
	51	fIsoConeR(0.4),
965c985f	52	fNDimensions(7),
965c985f	53	fECut(3.),
bd092f0f	54	fESD(0),
	55	fOutputList(0),
	56	fEvtSel(0),
	57	fPVtxZ(0),
	58	fCellAbsIdVsAmpl(0),
	59	fNClusHighClusE(0),
965c985f	60	fHnOutput(0)
bd092f0f	61	{
	62	// Default constructor.
	63	}
	64
	65	//________________________________________________________________________
	66	AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
	67	AliAnalysisTaskSE(name),
30159e6f	68	fCaloClusters(0),
30159e6f	69	fSelPrimTracks(0),
3f4073ba	70	fTracks(0),
30159e6f	71	fEMCalCells(0),
	72	fPrTrCuts(0),
	73	fGeom(0x0),
	74	fGeoName("EMCAL_COMPLETEV1"),
	75	fPeriod("LHC11c"),
751194e8	76	fTrigBit("kEMC7"),
30159e6f	77	fIsTrain(0),
	78	fExoticCut(0.97),
	79	fIsoConeR(0.4),
965c985f	80	fNDimensions(7),
965c985f	81	fECut(3.),
30159e6f	82	fESD(0),
30159e6f	83	fOutputList(0),
30159e6f	84	fEvtSel(0),
bd092f0f	85	fPVtxZ(0),
	86	fCellAbsIdVsAmpl(0),
	87	fNClusHighClusE(0),
965c985f	88	fHnOutput(0)
30159e6f	89	{
	90	// Constructor
	91
	92	// Define input and output slots here
	93	// Input slot #0 works with a TChain
	94	DefineInput(0, TChain::Class());
	95	// Output slot #0 id reserved by the base class for AOD
	96	// Output slot #1 writes into a TH1 container
	97	DefineOutput(1, TList::Class());
	98	}
bd092f0f	99
30159e6f	100	//________________________________________________________________________
	101	void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
	102	{
bd092f0f	103	// Create histograms, called once.
30159e6f	104
	105	fCaloClusters = new TRefArray();
	106	fSelPrimTracks = new TObjArray();
	107
	108
	109	fOutputList = new TList();
a62631a9	110	fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
30159e6f	111
	112	fGeom = AliEMCALGeometry::GetInstance(fGeoName);
	113
	114	fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
	115	fOutputList->Add(fEvtSel);
	116
	117
	118	fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
	119	fOutputList->Add(fPVtxZ);
	120
	121	fCellAbsIdVsAmpl = new TH2F("hCellAbsIdVsAmpl","cell abs id vs cell amplitude (energy);E (GeV);ID",200,0,100,244810,-0.5,244810-0.5);
a62631a9	122	fOutputList->Add(fCellAbsIdVsAmpl);
30159e6f	123
	124	fNClusHighClusE = new TH2F("hNClusHighClusE","total number of clusters vs. highest clus energy in the event;E (GeV);NClus",200,0,100,301,-0.5,300.5);
	125	fOutputList->Add(fNClusHighClusE);
	126
965c985f	127	const Int_t ndims = fNDimensions;
	128	Int_t nEt=1000, nM02=400, nCeIso=1000, nTrIso=1000, nAllIso=1000, nTrClDphi=200, nTrClDeta=100;
	129	Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nTrClDphi, nTrClDeta};
	130	Double_t xmin[] = { 0., 0., -10., -10., -10., -0.1,-0.05};
	131	Double_t xmax[] = { 100., 4., 190., 190., 190., 0.1, 0.05};
	132	fHnOutput = new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, d#phi_{trk},d#eta_{trk}", ndims, bins, xmin, xmax);
	133	fOutputList->Add(fHnOutput);
	134
	135
	136
30159e6f	137	PostData(1, fOutputList);
	138	}
	139
	140	//________________________________________________________________________
	141	void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
	142	{
bd092f0f	143	// Main loop, called for each event.
	144
	145	// event trigger selection
30159e6f	146	Bool_t isSelected = 0;
751194e8	147	if(fPeriod.Contains("11a")){
	148	if(fTrigBit.Contains("kEMC"))
	149	isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
	150	else
	151	isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
	152	}
	153	else{
	154	if(fTrigBit.Contains("kEMC"))
	155	isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
	156	else
	157	isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
	158	}
30159e6f	159	if(!isSelected )
30159e6f	160	return;
30159e6f	161
30159e6f	162	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
	163	if (!fESD) {
	164	printf("ERROR: fESD not available\n");
	165	return;
	166	}
	167
	168	fEvtSel->Fill(0);
	169
	170	AliESDVertex pv = (AliESDVertex)fESD->GetPrimaryVertex();
	171	if(!pv)
	172	return;
	173	fPVtxZ->Fill(pv->GetZ());
	174	if(TMath::Abs(pv->GetZ())>15)
	175	return;
	176
	177	fEvtSel->Fill(1);
bd092f0f	178
3f4073ba	179	if (!fTracks)
3f4073ba	180	fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
30159e6f	181	// Track loop to fill a pT spectrum
3f4073ba	182	const Int_t Ntracks = fTracks->GetEntriesFast();
	183	for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
	184	// for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
	185	//AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
	186	AliVTrack track = static_cast<AliVTrack>(fTracks->At(iTracks));
30159e6f	187	if (!track)
	188	continue;
	189	if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
	190	fSelPrimTracks->Add(track);
	191	//printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi());
	192	}
bd092f0f	193	}
bd092f0f	194
30159e6f	195	if(!fIsTrain){
	196	for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
	197	if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
	198	break;
	199	fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
	200	}
	201	}
	202	fESD->GetEMCALClusters(fCaloClusters);
	203	fEMCalCells = fESD->GetEMCALCells();
	204
	205
	206	FillClusHists();
	207
	208	fCaloClusters->Clear();
	209	fSelPrimTracks->Clear();
	210	PostData(1, fOutputList);
	211	}
bd092f0f	212
30159e6f	213	//________________________________________________________________________
	214	void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
	215	{
bd092f0f	216	// Fill cluster histograms.
bd092f0f	217
30159e6f	218	if(!fCaloClusters)
	219	return;
	220	const Int_t nclus = fCaloClusters->GetEntries();
	221	if(nclus==0)
	222	return;
e681d9ce	223	Double_t maxE = 0;
30159e6f	224	for(Int_t ic=0;ic<nclus;ic++){
	225	maxE=0;
	226	AliESDCaloCluster c = static_cast<AliESDCaloCluster>(fCaloClusters->At(ic));
	227	if(!c)
	228	continue;
	229	if(!c->IsEMCAL())
	230	continue;
965c985f	231	if(c->E()<fECut)
965c985f	232	continue;
30159e6f	233	Short_t id;
	234	Double_t Emax = GetMaxCellEnergy( c, id);
	235	Double_t Ecross = GetCrossEnergy( c, id);
	236	if((1-Ecross/Emax)>fExoticCut)
	237	continue;
	238	Float_t clsPos[3] = {0,0,0};
	239	c->GetPosition(clsPos);
	240	TVector3 clsVec(clsPos);
3dc2825e	241	Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
30159e6f	242	Float_t ceiso, cephiband, cecore;
	243	Float_t triso, trphiband, trcore;
	244	Float_t alliso, allphiband, allcore;
	245	Float_t phibandArea = (1.4 - 2fIsoConeR)2*fIsoConeR;
	246	Float_t netConeArea = TMath::Pi()(fIsoConeRfIsoConeR - 0.04*0.04);
	247	GetCeIso(clsVec, ceiso, cephiband, cecore);
	248	GetTrIso(clsVec, triso, trphiband, trcore);
	249	alliso = ceiso + triso;
	250	allphiband = cephiband + trphiband;
	251	allcore = cecore + trcore;
	252	Float_t ceisoue = cephiband/phibandArea*netConeArea;
	253	Float_t trisoue = trphiband/phibandArea*netConeArea;
	254	Float_t allisoue = allphiband/phibandArea*netConeArea;
965c985f	255	const Int_t ndims = fNDimensions;
	256	Double_t outputValues[ndims];
	257	outputValues[0] = Et;
	258	outputValues[1] = c->GetM02();
	259	outputValues[2] = ceiso-cecore-ceisoue;
	260	outputValues[3] = triso-trisoue;
	261	outputValues[4] = alliso-cecore-allisoue;
	262	outputValues[5] = c->GetTrackDx();
	263	outputValues[6] = c->GetTrackDz();
	264	fHnOutput->Fill(outputValues);
30159e6f	265	if(c->E()>maxE)
	266	maxE = c->E();
	267	}
	268	fNClusHighClusE->Fill(maxE,nclus);
	269	}
bd092f0f	270
30159e6f	271	//________________________________________________________________________
	272	void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
	273	{
bd092f0f	274	// Get cell isolation.
bd092f0f	275
30159e6f	276	if(!fEMCalCells)
	277	return;
	278	const Int_t ncells = fEMCalCells->GetNumberOfCells();
	279	Float_t totiso=0;
	280	Float_t totband=0;
	281	Float_t totcore=0;
	282	Float_t etacl = vec.Eta();
	283	Float_t phicl = vec.Phi();
	284	Float_t thetacl = vec.Theta();
	285	if(phicl<0)
	286	phicl+=TMath::TwoPi();
	287	Int_t absid = -1;
	288	Float_t eta=-1, phi=-1;
	289	for(int icell=0;icell<ncells;icell++){
	290	absid = TMath::Abs(fEMCalCells->GetCellNumber(icell));
	291	if(!fGeom)
	292	return;
	293	fGeom->EtaPhiFromIndex(absid,eta,phi);
	294	Float_t dphi = TMath::Abs(phi-phicl);
	295	Float_t deta = TMath::Abs(eta-etacl);
	296	Float_t R = TMath::Sqrt(detadeta + dphidphi);
	297	Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
	298	if(R<fIsoConeR){
	299	totiso += etcell;
	300	if(R<0.04)
	301	totcore += etcell;
	302	}
	303	else{
	304	if(dphi>fIsoConeR)
	305	continue;
	306	if(deta<fIsoConeR)
	307	continue;
	308	totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
	309	}
	310	}
	311	iso = totiso;
	312	phiband = totband;
	313	core = totcore;
	314	}
bd092f0f	315
30159e6f	316	//________________________________________________________________________
	317	void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
	318	{
bd092f0f	319	// Get track isolation.
bd092f0f	320
30159e6f	321	if(!fSelPrimTracks)
	322	return;
	323	const Int_t ntracks = fSelPrimTracks->GetEntries();
	324	Double_t totiso=0;
	325	Double_t totband=0;
	326	Double_t totcore=0;
	327	Double_t etacl = vec.Eta();
	328	Double_t phicl = vec.Phi();
	329	if(phicl<0)
	330	phicl+=TMath::TwoPi();
	331	for(int itrack=0;itrack<ntracks;itrack++){
3f4073ba	332	AliVTrack track = static_cast<AliVTrack> (fSelPrimTracks->At(itrack));
30159e6f	333	if(!track)
	334	continue;
	335	Double_t dphi = TMath::Abs(track->Phi()-phicl);
	336	Double_t deta = TMath::Abs(track->Eta()-etacl);
	337	Double_t R = TMath::Sqrt(detadeta + dphidphi);
	338	Double_t pt = track->Pt();
	339	if(R<fIsoConeR){
	340	totiso += track->Pt();
	341	if(R<0.04)
	342	totcore += pt;
	343	}
	344	else{
	345	if(dphi>fIsoConeR)
	346	continue;
	347	if(deta<fIsoConeR)
	348	continue;
	349	totband += track->Pt();
	350	}
	351	}
	352	iso = totiso;
	353	phiband = totband;
	354	core = totcore;
	355	}
bd092f0f	356
30159e6f	357	//________________________________________________________________________
	358	Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
	359	{
	360	// Calculate the energy of cross cells around the leading cell.
	361
	362	AliVCaloCells *cells = 0;
	363	cells = fESD->GetEMCALCells();
	364	if (!cells)
	365	return 0;
	366
30159e6f	367	if (!fGeom)
	368	return 0;
	369
	370	Int_t iSupMod = -1;
	371	Int_t iTower = -1;
	372	Int_t iIphi = -1;
	373	Int_t iIeta = -1;
	374	Int_t iphi = -1;
	375	Int_t ieta = -1;
	376	Int_t iphis = -1;
	377	Int_t ietas = -1;
	378
	379	Double_t crossEnergy = 0;
	380
	381	fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
	382	fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
	383
	384	Int_t ncells = cluster->GetNCells();
	385	for (Int_t i=0; i<ncells; i++) {
	386	Int_t cellAbsId = cluster->GetCellAbsId(i);
	387	fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
	388	fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
	389	Int_t aphidiff = TMath::Abs(iphi-iphis);
	390	if (aphidiff>1)
	391	continue;
	392	Int_t aetadiff = TMath::Abs(ieta-ietas);
	393	if (aetadiff>1)
	394	continue;
	395	if ( (aphidiff==1 && aetadiff==0) \|\|
	396	(aphidiff==0 && aetadiff==1) ) {
	397	crossEnergy += cells->GetCellAmplitude(cellAbsId);
	398	}
	399	}
	400
	401	return crossEnergy;
	402	}
	403
30159e6f	404	//________________________________________________________________________
	405	Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
	406	{
	407	// Get maximum energy of attached cell.
	408
	409	id = -1;
	410
	411	AliVCaloCells *cells = 0;
	412	cells = fESD->GetEMCALCells();
	413	if (!cells)
	414	return 0;
	415
	416	Double_t maxe = 0;
	417	Int_t ncells = cluster->GetNCells();
	418	for (Int_t i=0; i<ncells; i++) {
	419	Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
	420	if (e>maxe) {
	421	maxe = e;
	422	id = cluster->GetCellAbsId(i);
	423	}
	424	}
	425	return maxe;
	426	}
bd092f0f	427
30159e6f	428	//________________________________________________________________________
	429	void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
	430	{
bd092f0f	431	// Called once at the end of the query.
30159e6f	432	}