[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, nClEta=140, nClPhi=128;
  Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nTrClDphi, nTrClDeta,nClEta,nClPhi};
  Double_t xmin[] = { 0.,   0.,  -10.,   -10., -10., -0.1,-0.05, -0.7, 1.4};
  Double_t xmax[] = { 100., 4., 190., 190., 190., 0.1, 0.05, 0.7, 3.192};
  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();
    outputValues[7] = clsVec.Eta();
    outputValues[8] = clsVec.Phi();
    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;
0cbf2a1e	128	Int_t nEt=1000, nM02=400, nCeIso=1000, nTrIso=1000, nAllIso=1000, nTrClDphi=200, nTrClDeta=100, nClEta=140, nClPhi=128;
	129	Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nTrClDphi, nTrClDeta,nClEta,nClPhi};
	130	Double_t xmin[] = { 0., 0., -10., -10., -10., -0.1,-0.05, -0.7, 1.4};
	131	Double_t xmax[] = { 100., 4., 190., 190., 190., 0.1, 0.05, 0.7, 3.192};
965c985f	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();
0cbf2a1e	264	outputValues[7] = clsVec.Eta();
0cbf2a1e	265	outputValues[8] = clsVec.Phi();
965c985f	266	fHnOutput->Fill(outputValues);
30159e6f	267	if(c->E()>maxE)
	268	maxE = c->E();
	269	}
	270	fNClusHighClusE->Fill(maxE,nclus);
	271	}
bd092f0f	272
30159e6f	273	//________________________________________________________________________
	274	void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
	275	{
bd092f0f	276	// Get cell isolation.
bd092f0f	277
30159e6f	278	if(!fEMCalCells)
	279	return;
	280	const Int_t ncells = fEMCalCells->GetNumberOfCells();
	281	Float_t totiso=0;
	282	Float_t totband=0;
	283	Float_t totcore=0;
	284	Float_t etacl = vec.Eta();
	285	Float_t phicl = vec.Phi();
	286	Float_t thetacl = vec.Theta();
	287	if(phicl<0)
	288	phicl+=TMath::TwoPi();
	289	Int_t absid = -1;
	290	Float_t eta=-1, phi=-1;
	291	for(int icell=0;icell<ncells;icell++){
	292	absid = TMath::Abs(fEMCalCells->GetCellNumber(icell));
	293	if(!fGeom)
	294	return;
	295	fGeom->EtaPhiFromIndex(absid,eta,phi);
	296	Float_t dphi = TMath::Abs(phi-phicl);
	297	Float_t deta = TMath::Abs(eta-etacl);
	298	Float_t R = TMath::Sqrt(detadeta + dphidphi);
	299	Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
	300	if(R<fIsoConeR){
	301	totiso += etcell;
	302	if(R<0.04)
	303	totcore += etcell;
	304	}
	305	else{
	306	if(dphi>fIsoConeR)
	307	continue;
	308	if(deta<fIsoConeR)
	309	continue;
	310	totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
	311	}
	312	}
	313	iso = totiso;
	314	phiband = totband;
	315	core = totcore;
	316	}
bd092f0f	317
30159e6f	318	//________________________________________________________________________
	319	void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
	320	{
bd092f0f	321	// Get track isolation.
bd092f0f	322
30159e6f	323	if(!fSelPrimTracks)
	324	return;
	325	const Int_t ntracks = fSelPrimTracks->GetEntries();
	326	Double_t totiso=0;
	327	Double_t totband=0;
	328	Double_t totcore=0;
	329	Double_t etacl = vec.Eta();
	330	Double_t phicl = vec.Phi();
	331	if(phicl<0)
	332	phicl+=TMath::TwoPi();
	333	for(int itrack=0;itrack<ntracks;itrack++){
3f4073ba	334	AliVTrack track = static_cast<AliVTrack> (fSelPrimTracks->At(itrack));
30159e6f	335	if(!track)
	336	continue;
	337	Double_t dphi = TMath::Abs(track->Phi()-phicl);
	338	Double_t deta = TMath::Abs(track->Eta()-etacl);
	339	Double_t R = TMath::Sqrt(detadeta + dphidphi);
	340	Double_t pt = track->Pt();
	341	if(R<fIsoConeR){
	342	totiso += track->Pt();
	343	if(R<0.04)
	344	totcore += pt;
	345	}
	346	else{
	347	if(dphi>fIsoConeR)
	348	continue;
	349	if(deta<fIsoConeR)
	350	continue;
	351	totband += track->Pt();
	352	}
	353	}
	354	iso = totiso;
	355	phiband = totband;
	356	core = totcore;
	357	}
bd092f0f	358
30159e6f	359	//________________________________________________________________________
	360	Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
	361	{
	362	// Calculate the energy of cross cells around the leading cell.
	363
	364	AliVCaloCells *cells = 0;
	365	cells = fESD->GetEMCALCells();
	366	if (!cells)
	367	return 0;
	368
30159e6f	369	if (!fGeom)
	370	return 0;
	371
	372	Int_t iSupMod = -1;
	373	Int_t iTower = -1;
	374	Int_t iIphi = -1;
	375	Int_t iIeta = -1;
	376	Int_t iphi = -1;
	377	Int_t ieta = -1;
	378	Int_t iphis = -1;
	379	Int_t ietas = -1;
	380
	381	Double_t crossEnergy = 0;
	382
	383	fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
	384	fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
	385
	386	Int_t ncells = cluster->GetNCells();
	387	for (Int_t i=0; i<ncells; i++) {
	388	Int_t cellAbsId = cluster->GetCellAbsId(i);
	389	fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
	390	fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
	391	Int_t aphidiff = TMath::Abs(iphi-iphis);
	392	if (aphidiff>1)
	393	continue;
	394	Int_t aetadiff = TMath::Abs(ieta-ietas);
	395	if (aetadiff>1)
	396	continue;
	397	if ( (aphidiff==1 && aetadiff==0) \|\|
	398	(aphidiff==0 && aetadiff==1) ) {
	399	crossEnergy += cells->GetCellAmplitude(cellAbsId);
	400	}
	401	}
	402
	403	return crossEnergy;
	404	}
	405
30159e6f	406	//________________________________________________________________________
	407	Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
	408	{
	409	// Get maximum energy of attached cell.
	410
	411	id = -1;
	412
	413	AliVCaloCells *cells = 0;
	414	cells = fESD->GetEMCALCells();
	415	if (!cells)
	416	return 0;
	417
	418	Double_t maxe = 0;
	419	Int_t ncells = cluster->GetNCells();
	420	for (Int_t i=0; i<ncells; i++) {
	421	Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
	422	if (e>maxe) {
	423	maxe = e;
	424	id = cluster->GetCellAbsId(i);
	425	}
	426	}
	427	return maxe;
	428	}
bd092f0f	429
30159e6f	430	//________________________________________________________________________
	431	void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
	432	{
bd092f0f	433	// Called once at the end of the query.
30159e6f	434	}