[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 <TLorentzVector.h>
#include <TTree.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),
  fEMCalCells(0),
  fPrTrCuts(0),
  fGeom(0x0),
  fESD(0),
  fOutputList(0),
  fEvtSel(0),
  fPVtxZ(0),                 
  fCellAbsIdVsAmpl(0),       
  fNClusHighClusE(0),        
  fM02Et(0),                 
  fM02EtTM(0),               
  fM02EtCeIso1(0),           
  fM02EtCeIso2(0),           
  fM02EtCeIso5(0),           
  fM02EtTrIso1(0),           
  fM02EtTrIso2(0),           
  fM02EtTrIso5(0),           
  fM02EtAllIso1(0),          
  fM02EtAllIso2(0),          
  fM02EtAllIso5(0),          
  fCeIsoVsEtPho(0),          
  fTrIsoVsEtPho(0),          
  fAllIsoVsEtPho(0),         
  fM02EtCeIso1TM(0),         
  fM02EtCeIso2TM(0),         
  fM02EtCeIso5TM(0),         
  fM02EtTrIso1TM(0),         
  fM02EtTrIso2TM(0),         
  fM02EtTrIso5TM(0),         
  fM02EtAllIso1TM(0),        
  fM02EtAllIso2TM(0),        
  fM02EtAllIso5TM(0),        
  fCeIsoVsEtPhoTM(0),        
  fTrIsoVsEtPhoTM(0),        
  fAllIsoVsEtPhoTM(0)    
{
  // Default constructor.
}

//________________________________________________________________________
AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) : 
  AliAnalysisTaskSE(name), 
  fCaloClusters(0),
  fSelPrimTracks(0),
  fEMCalCells(0),
  fPrTrCuts(0),
  fGeom(0x0),
  fGeoName("EMCAL_COMPLETEV1"),
  fPeriod("LHC11c"),
  fIsTrain(0),
  fExoticCut(0.97),
  fIsoConeR(0.4),
  fESD(0),
  fOutputList(0),
  fEvtSel(0),
  fPVtxZ(0),            
  fCellAbsIdVsAmpl(0),  
  fNClusHighClusE(0),   
  fM02Et(0),            
  fM02EtTM(0),          
  fM02EtCeIso1(0),      
  fM02EtCeIso2(0),      
  fM02EtCeIso5(0),      
  fM02EtTrIso1(0),      
  fM02EtTrIso2(0),      
  fM02EtTrIso5(0),      
  fM02EtAllIso1(0),     
  fM02EtAllIso2(0),     
  fM02EtAllIso5(0),     
  fCeIsoVsEtPho(0),     
  fTrIsoVsEtPho(0),     
  fAllIsoVsEtPho(0),    
  fM02EtCeIso1TM(0),
  fM02EtCeIso2TM(0),        
  fM02EtCeIso5TM(0),        
  fM02EtTrIso1TM(0),        
  fM02EtTrIso2TM(0),        
  fM02EtTrIso5TM(0),        
  fM02EtAllIso1TM(0),       
  fM02EtAllIso2TM(0),       
  fM02EtAllIso5TM(0),       
  fCeIsoVsEtPhoTM(0),       
  fTrIsoVsEtPhoTM(0),       
  fAllIsoVsEtPhoTM(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(fPVtxZ);

  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);

  fM02Et = new TH2F("fM02Et","M02 vs Et for all clusters;E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02Et->Sumw2();
  fOutputList->Add(fM02Et);

  fM02EtTM = new TH2F("fM02EtTM","M02 vs Et for all track-matched clusters;E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtTM->Sumw2();
  fOutputList->Add(fM02EtTM);

  fM02EtCeIso1 = new TH2F("fM02EtCeIso1","M02 vs Et for all clusters (ISO_{EMC}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtCeIso1->Sumw2();
  fOutputList->Add(fM02EtCeIso1);

  fM02EtCeIso2 = new TH2F("fM02EtCeIso2","M02 vs Et for all clusters (ISO_{EMC}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtCeIso2->Sumw2();
  fOutputList->Add(fM02EtCeIso2);

  fM02EtCeIso5 = new TH2F("fM02EtCeIso5","M02 vs Et for all clusters (ISO_{EMC}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtCeIso5->Sumw2();
  fOutputList->Add(fM02EtCeIso5);

  fM02EtTrIso1 = new TH2F("fM02EtTrIso1","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtTrIso1->Sumw2();
  fOutputList->Add(fM02EtTrIso1);

  fM02EtTrIso2 = new TH2F("fM02EtTrIso2","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtTrIso2->Sumw2();
  fOutputList->Add(fM02EtTrIso2);

  fM02EtTrIso5 = new TH2F("fM02EtTrIso5","M02 vs Et for all clusters (ISO_{TRK}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtTrIso5->Sumw2();
  fOutputList->Add(fM02EtTrIso5);

  fM02EtAllIso1 = new TH2F("fM02EtAllIso1","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtAllIso1->Sumw2();
  fOutputList->Add(fM02EtAllIso1);

  fM02EtAllIso2 = new TH2F("fM02EtAllIso2","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtAllIso2->Sumw2();
  fOutputList->Add(fM02EtAllIso2);

  fM02EtAllIso5 = new TH2F("fM02EtAllIso5","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtAllIso5->Sumw2();
  fOutputList->Add(fM02EtAllIso5);

  fCeIsoVsEtPho = new TH2F("fCeIsoVsEtPho","ISO_{EMC} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC}",1000,0,100,1000,-10,190);
  fCeIsoVsEtPho->Sumw2();
  fOutputList->Add(fCeIsoVsEtPho);

  fTrIsoVsEtPho = new TH2F("fTrIsoVsEtPho","ISO_{TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{TRK}",1000,0,100,1000,-10,190);
  fTrIsoVsEtPho->Sumw2();
  fOutputList->Add(fTrIsoVsEtPho);

  fAllIsoVsEtPho = new TH2F("fAllIsoVsEtPho","ISO_{EMC+TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC+TRK}",1000,0,100,1000,-10,190);
  fAllIsoVsEtPho->Sumw2();
  fOutputList->Add(fAllIsoVsEtPho);

  fM02EtCeIso1TM = new TH2F("fM02EtCeIso1TM","M02 vs Et for all clusters (ISO_{EMC}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtCeIso1TM->Sumw2();
  fOutputList->Add(fM02EtCeIso1TM);

  fM02EtCeIso2TM = new TH2F("fM02EtCeIso2TM","M02 vs Et for all clusters (ISO_{EMC}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtCeIso2TM->Sumw2();
  fOutputList->Add(fM02EtCeIso2TM);

  fM02EtCeIso5TM = new TH2F("fM02EtCeIso5TM","M02 vs Et for all clusters (ISO_{EMC}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtCeIso5TM->Sumw2();
  fOutputList->Add(fM02EtCeIso5TM);

  fM02EtTrIso1TM = new TH2F("fM02EtTrIso1TM","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtTrIso1TM->Sumw2();
  fOutputList->Add(fM02EtTrIso1TM);

  fM02EtTrIso2TM = new TH2F("fM02EtTrIso2TM","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtTrIso2TM->Sumw2();
  fOutputList->Add(fM02EtTrIso2TM);

  fM02EtTrIso5TM = new TH2F("fM02EtTrIso5TM","M02 vs Et for all clusters (ISO_{TRK}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtTrIso5TM->Sumw2();
  fOutputList->Add(fM02EtTrIso5TM);

  fM02EtAllIso1TM = new TH2F("fM02EtAllIso1TM","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtAllIso1TM->Sumw2();
  fOutputList->Add(fM02EtAllIso1TM);

  fM02EtAllIso2TM = new TH2F("fM02EtAllIso2TM","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtAllIso2TM->Sumw2();
  fOutputList->Add(fM02EtAllIso2TM);

  fM02EtAllIso5TM = new TH2F("fM02EtAllIso5TM","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
  fM02EtAllIso5TM->Sumw2();
  fOutputList->Add(fM02EtAllIso5TM);

  fCeIsoVsEtPhoTM = new TH2F("fCeIsoVsEtPhoTM","ISO_{EMC} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC}",1000,0,100,1000,-10,190);
  fCeIsoVsEtPhoTM->Sumw2();
  fOutputList->Add(fCeIsoVsEtPhoTM);

  fTrIsoVsEtPhoTM = new TH2F("fTrIsoVsEtPhoTM","ISO_{TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{TRK}",1000,0,100,1000,-10,190);
  fTrIsoVsEtPhoTM->Sumw2();
  fOutputList->Add(fTrIsoVsEtPhoTM);

  fAllIsoVsEtPhoTM = new TH2F("fAllIsoVsEtPhoTM","ISO_{EMC+TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC+TRK}",1000,0,100,1000,-10,190);
  fAllIsoVsEtPhoTM->Sumw2();
  fOutputList->Add(fAllIsoVsEtPhoTM);

  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"))
    isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
  else
    isSelected =  (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
  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);

  // Track loop to fill a pT spectrum
  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
    AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(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;
    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;
    fM02Et->Fill(Et, c->GetM02());
    if(ceiso-cecore<1)
      fM02EtCeIso1->Fill(Et, c->GetM02());
    if(ceiso-cecore<2)
      fM02EtCeIso2->Fill(Et, c->GetM02());
    if(ceiso-cecore<5)
      fM02EtCeIso5->Fill(Et, c->GetM02());
    if(triso-trcore<1)
      fM02EtTrIso1->Fill(Et, c->GetM02());
    if(triso-trcore<2)
      fM02EtTrIso2->Fill(Et, c->GetM02());
    if(triso-trcore<5)
      fM02EtTrIso5->Fill(Et, c->GetM02());
    if(alliso-allcore<1)
      fM02EtAllIso1->Fill(Et, c->GetM02());
    if(alliso-allcore<2)
      fM02EtAllIso2->Fill(Et, c->GetM02());
    if(alliso-allcore<5)
      fM02EtAllIso5->Fill(Et, c->GetM02());
    if(c->GetM02()>0.1 && c->GetM02()<0.3){
      fCeIsoVsEtPho->Fill(Et, ceiso - cecore - ceisoue);
      fTrIsoVsEtPho->Fill(Et, triso - trcore - trisoue);
      fAllIsoVsEtPho->Fill(Et, alliso - allcore - allisoue);
      }
    Double_t dR = TMath::Sqrt(pow(c->GetTrackDx(),2)+pow(c->GetTrackDz(),2));
    if(dR<0.014){
      fM02EtTM->Fill(Et, c->GetM02());
      if(ceiso-cecore<1)
	fM02EtCeIso1TM->Fill(Et, c->GetM02());
      if(ceiso-cecore<2)
	fM02EtCeIso2TM->Fill(Et, c->GetM02());
      if(ceiso-cecore<5)
	fM02EtCeIso5TM->Fill(Et, c->GetM02());
      if(triso-trcore<1)
	fM02EtTrIso1TM->Fill(Et, c->GetM02());
      if(triso-trcore<2)
	fM02EtTrIso2TM->Fill(Et, c->GetM02());
      if(triso-trcore<5)
	fM02EtTrIso5TM->Fill(Et, c->GetM02());
      if(alliso-allcore<1)
	fM02EtAllIso1TM->Fill(Et, c->GetM02());
      if(alliso-allcore<2)
	fM02EtAllIso2TM->Fill(Et, c->GetM02());
      if(alliso-allcore<5)
	fM02EtAllIso5TM->Fill(Et, c->GetM02());
      if(c->GetM02()>0.1 && c->GetM02()<0.3){
	fCeIsoVsEtPhoTM->Fill(Et, ceiso);
	fTrIsoVsEtPhoTM->Fill(Et, triso);
	fAllIsoVsEtPhoTM->Fill(Et, alliso);
      }
    }
    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++){
    AliESDtrack *track = static_cast<AliESDtrack*> (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>
	10	#include <TLorentzVector.h>
	11	#include <TTree.h>
	12
	13	#include "AliAnalysisManager.h"
	14	#include "AliAnalysisTask.h"
	15	#include "AliEMCALGeometry.h"
	16	#include "AliEMCALRecoUtils.h"
	17	#include "AliESDCaloCells.h"
	18	#include "AliESDCaloCluster.h"
30159e6f	19	#include "AliESDEvent.h"
30159e6f	20	#include "AliESDHeader.h"
30159e6f	21	#include "AliESDInputHandler.h"
bd092f0f	22	#include "AliESDUtils.h"
30159e6f	23	#include "AliESDpid.h"
30159e6f	24	#include "AliESDtrack.h"
	25	#include "AliESDtrackCuts.h"
	26	#include "AliESDv0.h"
bd092f0f	27	#include "AliKFParticle.h"
	28	#include "AliMCEvent.h"
	29	#include "AliMCEventHandler.h"
	30	#include "AliStack.h"
30159e6f	31	#include "AliV0vertexer.h"
30159e6f	32	#include "AliVCluster.h"
30159e6f	33
30159e6f	34	ClassImp(AliAnalysisTaskEMCALIsoPhoton)
	35
	36	//________________________________________________________________________
bd092f0f	37	AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() :
	38	AliAnalysisTaskSE(),
	39	fCaloClusters(0),
	40	fSelPrimTracks(0),
	41	fEMCalCells(0),
	42	fPrTrCuts(0),
	43	fGeom(0x0),
	44	fESD(0),
	45	fOutputList(0),
	46	fEvtSel(0),
	47	fPVtxZ(0),
	48	fCellAbsIdVsAmpl(0),
	49	fNClusHighClusE(0),
	50	fM02Et(0),
	51	fM02EtTM(0),
	52	fM02EtCeIso1(0),
	53	fM02EtCeIso2(0),
	54	fM02EtCeIso5(0),
	55	fM02EtTrIso1(0),
	56	fM02EtTrIso2(0),
	57	fM02EtTrIso5(0),
	58	fM02EtAllIso1(0),
	59	fM02EtAllIso2(0),
	60	fM02EtAllIso5(0),
	61	fCeIsoVsEtPho(0),
	62	fTrIsoVsEtPho(0),
	63	fAllIsoVsEtPho(0),
	64	fM02EtCeIso1TM(0),
	65	fM02EtCeIso2TM(0),
	66	fM02EtCeIso5TM(0),
	67	fM02EtTrIso1TM(0),
	68	fM02EtTrIso2TM(0),
	69	fM02EtTrIso5TM(0),
	70	fM02EtAllIso1TM(0),
	71	fM02EtAllIso2TM(0),
	72	fM02EtAllIso5TM(0),
	73	fCeIsoVsEtPhoTM(0),
	74	fTrIsoVsEtPhoTM(0),
	75	fAllIsoVsEtPhoTM(0)
	76	{
	77	// Default constructor.
	78	}
	79
	80	//________________________________________________________________________
	81	AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
	82	AliAnalysisTaskSE(name),
30159e6f	83	fCaloClusters(0),
	84	fSelPrimTracks(0),
	85	fEMCalCells(0),
	86	fPrTrCuts(0),
	87	fGeom(0x0),
	88	fGeoName("EMCAL_COMPLETEV1"),
	89	fPeriod("LHC11c"),
	90	fIsTrain(0),
	91	fExoticCut(0.97),
	92	fIsoConeR(0.4),
30159e6f	93	fESD(0),
30159e6f	94	fOutputList(0),
30159e6f	95	fEvtSel(0),
bd092f0f	96	fPVtxZ(0),
	97	fCellAbsIdVsAmpl(0),
	98	fNClusHighClusE(0),
	99	fM02Et(0),
	100	fM02EtTM(0),
	101	fM02EtCeIso1(0),
	102	fM02EtCeIso2(0),
	103	fM02EtCeIso5(0),
	104	fM02EtTrIso1(0),
	105	fM02EtTrIso2(0),
	106	fM02EtTrIso5(0),
	107	fM02EtAllIso1(0),
	108	fM02EtAllIso2(0),
	109	fM02EtAllIso5(0),
	110	fCeIsoVsEtPho(0),
	111	fTrIsoVsEtPho(0),
	112	fAllIsoVsEtPho(0),
	113	fM02EtCeIso1TM(0),
	114	fM02EtCeIso2TM(0),
	115	fM02EtCeIso5TM(0),
	116	fM02EtTrIso1TM(0),
	117	fM02EtTrIso2TM(0),
	118	fM02EtTrIso5TM(0),
	119	fM02EtAllIso1TM(0),
	120	fM02EtAllIso2TM(0),
	121	fM02EtAllIso5TM(0),
	122	fCeIsoVsEtPhoTM(0),
	123	fTrIsoVsEtPhoTM(0),
	124	fAllIsoVsEtPhoTM(0)
30159e6f	125	{
	126	// Constructor
	127
	128	// Define input and output slots here
	129	// Input slot #0 works with a TChain
	130	DefineInput(0, TChain::Class());
	131	// Output slot #0 id reserved by the base class for AOD
	132	// Output slot #1 writes into a TH1 container
	133	DefineOutput(1, TList::Class());
	134	}
bd092f0f	135
30159e6f	136	//________________________________________________________________________
	137	void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
	138	{
bd092f0f	139	// Create histograms, called once.
30159e6f	140
	141	fCaloClusters = new TRefArray();
	142	fSelPrimTracks = new TObjArray();
	143
	144
	145	fOutputList = new TList();
	146	fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
	147
	148	fGeom = AliEMCALGeometry::GetInstance(fGeoName);
	149
	150	fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
	151	fOutputList->Add(fEvtSel);
	152
	153
	154	fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
	155	fOutputList->Add(fPVtxZ);
	156
	157	fCellAbsIdVsAmpl = new TH2F("hCellAbsIdVsAmpl","cell abs id vs cell amplitude (energy);E (GeV);ID",200,0,100,244810,-0.5,244810-0.5);
	158	fOutputList->Add(fPVtxZ);
	159
	160	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);
	161	fOutputList->Add(fNClusHighClusE);
	162
	163	fM02Et = new TH2F("fM02Et","M02 vs Et for all clusters;E_{T} (GeV);M02",1000,0,100,400,0,4);
	164	fM02Et->Sumw2();
	165	fOutputList->Add(fM02Et);
	166
	167	fM02EtTM = new TH2F("fM02EtTM","M02 vs Et for all track-matched clusters;E_{T} (GeV);M02",1000,0,100,400,0,4);
	168	fM02EtTM->Sumw2();
	169	fOutputList->Add(fM02EtTM);
	170
	171	fM02EtCeIso1 = new TH2F("fM02EtCeIso1","M02 vs Et for all clusters (ISO_{EMC}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
	172	fM02EtCeIso1->Sumw2();
	173	fOutputList->Add(fM02EtCeIso1);
	174
	175	fM02EtCeIso2 = new TH2F("fM02EtCeIso2","M02 vs Et for all clusters (ISO_{EMC}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
	176	fM02EtCeIso2->Sumw2();
	177	fOutputList->Add(fM02EtCeIso2);
	178
	179	fM02EtCeIso5 = new TH2F("fM02EtCeIso5","M02 vs Et for all clusters (ISO_{EMC}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
	180	fM02EtCeIso5->Sumw2();
	181	fOutputList->Add(fM02EtCeIso5);
	182
	183	fM02EtTrIso1 = new TH2F("fM02EtTrIso1","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
	184	fM02EtTrIso1->Sumw2();
	185	fOutputList->Add(fM02EtTrIso1);
	186
	187	fM02EtTrIso2 = new TH2F("fM02EtTrIso2","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
	188	fM02EtTrIso2->Sumw2();
	189	fOutputList->Add(fM02EtTrIso2);
	190
	191	fM02EtTrIso5 = new TH2F("fM02EtTrIso5","M02 vs Et for all clusters (ISO_{TRK}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
	192	fM02EtTrIso5->Sumw2();
	193	fOutputList->Add(fM02EtTrIso5);
	194
	195	fM02EtAllIso1 = new TH2F("fM02EtAllIso1","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
	196	fM02EtAllIso1->Sumw2();
	197	fOutputList->Add(fM02EtAllIso1);
	198
	199	fM02EtAllIso2 = new TH2F("fM02EtAllIso2","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
	200	fM02EtAllIso2->Sumw2();
	201	fOutputList->Add(fM02EtAllIso2);
	202
	203	fM02EtAllIso5 = new TH2F("fM02EtAllIso5","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
204	fM02EtAllIso5->Sumw2();
205	fOutputList->Add(fM02EtAllIso5);
206
207	fCeIsoVsEtPho = new TH2F("fCeIsoVsEtPho","ISO_{EMC} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC}",1000,0,100,1000,-10,190);
208	fCeIsoVsEtPho->Sumw2();
209	fOutputList->Add(fCeIsoVsEtPho);
210
211	fTrIsoVsEtPho = new TH2F("fTrIsoVsEtPho","ISO_{TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{TRK}",1000,0,100,1000,-10,190);
212	fTrIsoVsEtPho->Sumw2();
213	fOutputList->Add(fTrIsoVsEtPho);
214
215	fAllIsoVsEtPho = new TH2F("fAllIsoVsEtPho","ISO_{EMC+TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC+TRK}",1000,0,100,1000,-10,190);
216	fAllIsoVsEtPho->Sumw2();
217	fOutputList->Add(fAllIsoVsEtPho);
218
219	fM02EtCeIso1TM = new TH2F("fM02EtCeIso1TM","M02 vs Et for all clusters (ISO_{EMC}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
220	fM02EtCeIso1TM->Sumw2();
221	fOutputList->Add(fM02EtCeIso1TM);
222
223	fM02EtCeIso2TM = new TH2F("fM02EtCeIso2TM","M02 vs Et for all clusters (ISO_{EMC}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
224	fM02EtCeIso2TM->Sumw2();
225	fOutputList->Add(fM02EtCeIso2TM);
226
227	fM02EtCeIso5TM = new TH2F("fM02EtCeIso5TM","M02 vs Et for all clusters (ISO_{EMC}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
228	fM02EtCeIso5TM->Sumw2();
229	fOutputList->Add(fM02EtCeIso5TM);
230
231	fM02EtTrIso1TM = new TH2F("fM02EtTrIso1TM","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
232	fM02EtTrIso1TM->Sumw2();
233	fOutputList->Add(fM02EtTrIso1TM);
234
235	fM02EtTrIso2TM = new TH2F("fM02EtTrIso2TM","M02 vs Et for all clusters (ISO_{TRK}<1GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
236	fM02EtTrIso2TM->Sumw2();
237	fOutputList->Add(fM02EtTrIso2TM);
238
239	fM02EtTrIso5TM = new TH2F("fM02EtTrIso5TM","M02 vs Et for all clusters (ISO_{TRK}<5GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
240	fM02EtTrIso5TM->Sumw2();
241	fOutputList->Add(fM02EtTrIso5TM);
242
243	fM02EtAllIso1TM = new TH2F("fM02EtAllIso1TM","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
244	fM02EtAllIso1TM->Sumw2();
245	fOutputList->Add(fM02EtAllIso1TM);
246
247	fM02EtAllIso2TM = new TH2F("fM02EtAllIso2TM","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
248	fM02EtAllIso2TM->Sumw2();
249	fOutputList->Add(fM02EtAllIso2TM);
250
251	fM02EtAllIso5TM = new TH2F("fM02EtAllIso5TM","M02 vs Et for all clusters (ISO_{EMC+TRK}<2GeV);E_{T} (GeV);M02",1000,0,100,400,0,4);
252	fM02EtAllIso5TM->Sumw2();
253	fOutputList->Add(fM02EtAllIso5TM);
254
255	fCeIsoVsEtPhoTM = new TH2F("fCeIsoVsEtPhoTM","ISO_{EMC} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC}",1000,0,100,1000,-10,190);
256	fCeIsoVsEtPhoTM->Sumw2();
257	fOutputList->Add(fCeIsoVsEtPhoTM);
258
259	fTrIsoVsEtPhoTM = new TH2F("fTrIsoVsEtPhoTM","ISO_{TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{TRK}",1000,0,100,1000,-10,190);
260	fTrIsoVsEtPhoTM->Sumw2();
261	fOutputList->Add(fTrIsoVsEtPhoTM);
262
263	fAllIsoVsEtPhoTM = new TH2F("fAllIsoVsEtPhoTM","ISO_{EMC+TRK} vs. E_{T}^{clus} (0.1<#lambda_{0}^{2}<0.3);E_{T} (GeV);ISO_{EMC+TRK}",1000,0,100,1000,-10,190);
264	fAllIsoVsEtPhoTM->Sumw2();
265	fOutputList->Add(fAllIsoVsEtPhoTM);
266
267	PostData(1, fOutputList);
268	}
269
270	//________________________________________________________________________
271	void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
272	{
bd092f0f	273	// Main loop, called for each event.
	274
	275	// event trigger selection
30159e6f	276	Bool_t isSelected = 0;
	277	if(fPeriod.Contains("11a"))
	278	isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
	279	else
	280	isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
	281	if(!isSelected )
	282	return;
30159e6f	283
30159e6f	284	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
	285	if (!fESD) {
	286	printf("ERROR: fESD not available\n");
	287	return;
	288	}
	289
	290	fEvtSel->Fill(0);
	291
	292	AliESDVertex pv = (AliESDVertex)fESD->GetPrimaryVertex();
	293	if(!pv)
	294	return;
	295	fPVtxZ->Fill(pv->GetZ());
	296	if(TMath::Abs(pv->GetZ())>15)
	297	return;
	298
	299	fEvtSel->Fill(1);
bd092f0f	300
30159e6f	301	// Track loop to fill a pT spectrum
	302	for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
	303	AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
	304	if (!track)
	305	continue;
	306	if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
	307	fSelPrimTracks->Add(track);
	308	//printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi());
	309	}
bd092f0f	310	}
bd092f0f	311
30159e6f	312	if(!fIsTrain){
	313	for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
	314	if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
	315	break;
	316	fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
	317	}
	318	}
	319	fESD->GetEMCALClusters(fCaloClusters);
	320	fEMCalCells = fESD->GetEMCALCells();
	321
	322
	323	FillClusHists();
	324
	325	fCaloClusters->Clear();
	326	fSelPrimTracks->Clear();
	327	PostData(1, fOutputList);
	328	}
bd092f0f	329
30159e6f	330	//________________________________________________________________________
	331	void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
	332	{
bd092f0f	333	// Fill cluster histograms.
bd092f0f	334
30159e6f	335	if(!fCaloClusters)
	336	return;
	337	const Int_t nclus = fCaloClusters->GetEntries();
	338	if(nclus==0)
	339	return;
e681d9ce	340	Double_t maxE = 0;
30159e6f	341	for(Int_t ic=0;ic<nclus;ic++){
	342	maxE=0;
	343	AliESDCaloCluster c = static_cast<AliESDCaloCluster>(fCaloClusters->At(ic));
	344	if(!c)
	345	continue;
	346	if(!c->IsEMCAL())
	347	continue;
	348	Short_t id;
	349	Double_t Emax = GetMaxCellEnergy( c, id);
	350	Double_t Ecross = GetCrossEnergy( c, id);
	351	if((1-Ecross/Emax)>fExoticCut)
	352	continue;
	353	Float_t clsPos[3] = {0,0,0};
	354	c->GetPosition(clsPos);
	355	TVector3 clsVec(clsPos);
3dc2825e	356	Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
30159e6f	357	Float_t ceiso, cephiband, cecore;
	358	Float_t triso, trphiband, trcore;
	359	Float_t alliso, allphiband, allcore;
	360	Float_t phibandArea = (1.4 - 2fIsoConeR)2*fIsoConeR;
	361	Float_t netConeArea = TMath::Pi()(fIsoConeRfIsoConeR - 0.04*0.04);
	362	GetCeIso(clsVec, ceiso, cephiband, cecore);
	363	GetTrIso(clsVec, triso, trphiband, trcore);
	364	alliso = ceiso + triso;
	365	allphiband = cephiband + trphiband;
	366	allcore = cecore + trcore;
	367	Float_t ceisoue = cephiband/phibandArea*netConeArea;
	368	Float_t trisoue = trphiband/phibandArea*netConeArea;
	369	Float_t allisoue = allphiband/phibandArea*netConeArea;
30159e6f	370	fM02Et->Fill(Et, c->GetM02());
70484e97	371	if(ceiso-cecore<1)
30159e6f	372	fM02EtCeIso1->Fill(Et, c->GetM02());
70484e97	373	if(ceiso-cecore<2)
30159e6f	374	fM02EtCeIso2->Fill(Et, c->GetM02());
70484e97	375	if(ceiso-cecore<5)
30159e6f	376	fM02EtCeIso5->Fill(Et, c->GetM02());
70484e97	377	if(triso-trcore<1)
30159e6f	378	fM02EtTrIso1->Fill(Et, c->GetM02());
70484e97	379	if(triso-trcore<2)
30159e6f	380	fM02EtTrIso2->Fill(Et, c->GetM02());
70484e97	381	if(triso-trcore<5)
30159e6f	382	fM02EtTrIso5->Fill(Et, c->GetM02());
70484e97	383	if(alliso-allcore<1)
30159e6f	384	fM02EtAllIso1->Fill(Et, c->GetM02());
70484e97	385	if(alliso-allcore<2)
30159e6f	386	fM02EtAllIso2->Fill(Et, c->GetM02());
70484e97	387	if(alliso-allcore<5)
30159e6f	388	fM02EtAllIso5->Fill(Et, c->GetM02());
	389	if(c->GetM02()>0.1 && c->GetM02()<0.3){
	390	fCeIsoVsEtPho->Fill(Et, ceiso - cecore - ceisoue);
	391	fTrIsoVsEtPho->Fill(Et, triso - trcore - trisoue);
	392	fAllIsoVsEtPho->Fill(Et, alliso - allcore - allisoue);
	393	}
	394	Double_t dR = TMath::Sqrt(pow(c->GetTrackDx(),2)+pow(c->GetTrackDz(),2));
	395	if(dR<0.014){
	396	fM02EtTM->Fill(Et, c->GetM02());
70484e97	397	if(ceiso-cecore<1)
30159e6f	398	fM02EtCeIso1TM->Fill(Et, c->GetM02());
70484e97	399	if(ceiso-cecore<2)
30159e6f	400	fM02EtCeIso2TM->Fill(Et, c->GetM02());
70484e97	401	if(ceiso-cecore<5)
30159e6f	402	fM02EtCeIso5TM->Fill(Et, c->GetM02());
70484e97	403	if(triso-trcore<1)
30159e6f	404	fM02EtTrIso1TM->Fill(Et, c->GetM02());
70484e97	405	if(triso-trcore<2)
30159e6f	406	fM02EtTrIso2TM->Fill(Et, c->GetM02());
70484e97	407	if(triso-trcore<5)
30159e6f	408	fM02EtTrIso5TM->Fill(Et, c->GetM02());
70484e97	409	if(alliso-allcore<1)
30159e6f	410	fM02EtAllIso1TM->Fill(Et, c->GetM02());
70484e97	411	if(alliso-allcore<2)
30159e6f	412	fM02EtAllIso2TM->Fill(Et, c->GetM02());
70484e97	413	if(alliso-allcore<5)
30159e6f	414	fM02EtAllIso5TM->Fill(Et, c->GetM02());
	415	if(c->GetM02()>0.1 && c->GetM02()<0.3){
	416	fCeIsoVsEtPhoTM->Fill(Et, ceiso);
	417	fTrIsoVsEtPhoTM->Fill(Et, triso);
	418	fAllIsoVsEtPhoTM->Fill(Et, alliso);
	419	}
	420	}
	421	if(c->E()>maxE)
	422	maxE = c->E();
	423	}
	424	fNClusHighClusE->Fill(maxE,nclus);
	425	}
bd092f0f	426
30159e6f	427	//________________________________________________________________________
	428	void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
	429	{
bd092f0f	430	// Get cell isolation.
bd092f0f	431
30159e6f	432	if(!fEMCalCells)
	433	return;
	434	const Int_t ncells = fEMCalCells->GetNumberOfCells();
	435	Float_t totiso=0;
	436	Float_t totband=0;
	437	Float_t totcore=0;
	438	Float_t etacl = vec.Eta();
	439	Float_t phicl = vec.Phi();
	440	Float_t thetacl = vec.Theta();
	441	if(phicl<0)
	442	phicl+=TMath::TwoPi();
	443	Int_t absid = -1;
	444	Float_t eta=-1, phi=-1;
	445	for(int icell=0;icell<ncells;icell++){
	446	absid = TMath::Abs(fEMCalCells->GetCellNumber(icell));
	447	if(!fGeom)
	448	return;
	449	fGeom->EtaPhiFromIndex(absid,eta,phi);
	450	Float_t dphi = TMath::Abs(phi-phicl);
	451	Float_t deta = TMath::Abs(eta-etacl);
	452	Float_t R = TMath::Sqrt(detadeta + dphidphi);
	453	Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
	454	if(R<fIsoConeR){
	455	totiso += etcell;
	456	if(R<0.04)
	457	totcore += etcell;
	458	}
	459	else{
	460	if(dphi>fIsoConeR)
	461	continue;
	462	if(deta<fIsoConeR)
	463	continue;
	464	totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
	465	}
	466	}
	467	iso = totiso;
	468	phiband = totband;
	469	core = totcore;
	470	}
bd092f0f	471
30159e6f	472	//________________________________________________________________________
	473	void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
	474	{
bd092f0f	475	// Get track isolation.
bd092f0f	476
30159e6f	477	if(!fSelPrimTracks)
	478	return;
	479	const Int_t ntracks = fSelPrimTracks->GetEntries();
	480	Double_t totiso=0;
	481	Double_t totband=0;
	482	Double_t totcore=0;
	483	Double_t etacl = vec.Eta();
	484	Double_t phicl = vec.Phi();
	485	if(phicl<0)
	486	phicl+=TMath::TwoPi();
	487	for(int itrack=0;itrack<ntracks;itrack++){
	488	AliESDtrack track = static_cast<AliESDtrack> (fSelPrimTracks->At(itrack));
	489	if(!track)
	490	continue;
	491	Double_t dphi = TMath::Abs(track->Phi()-phicl);
	492	Double_t deta = TMath::Abs(track->Eta()-etacl);
	493	Double_t R = TMath::Sqrt(detadeta + dphidphi);
	494	Double_t pt = track->Pt();
	495	if(R<fIsoConeR){
	496	totiso += track->Pt();
	497	if(R<0.04)
	498	totcore += pt;
	499	}
	500	else{
	501	if(dphi>fIsoConeR)
	502	continue;
	503	if(deta<fIsoConeR)
	504	continue;
	505	totband += track->Pt();
	506	}
	507	}
	508	iso = totiso;
	509	phiband = totband;
	510	core = totcore;
	511	}
bd092f0f	512
30159e6f	513	//________________________________________________________________________
	514	Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
	515	{
	516	// Calculate the energy of cross cells around the leading cell.
	517
	518	AliVCaloCells *cells = 0;
	519	cells = fESD->GetEMCALCells();
	520	if (!cells)
	521	return 0;
	522
30159e6f	523	if (!fGeom)
	524	return 0;
	525
	526	Int_t iSupMod = -1;
	527	Int_t iTower = -1;
	528	Int_t iIphi = -1;
	529	Int_t iIeta = -1;
	530	Int_t iphi = -1;
	531	Int_t ieta = -1;
	532	Int_t iphis = -1;
	533	Int_t ietas = -1;
	534
	535	Double_t crossEnergy = 0;
	536
	537	fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
	538	fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
	539
	540	Int_t ncells = cluster->GetNCells();
	541	for (Int_t i=0; i<ncells; i++) {
	542	Int_t cellAbsId = cluster->GetCellAbsId(i);
	543	fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
	544	fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
	545	Int_t aphidiff = TMath::Abs(iphi-iphis);
	546	if (aphidiff>1)
	547	continue;
	548	Int_t aetadiff = TMath::Abs(ieta-ietas);
	549	if (aetadiff>1)
	550	continue;
	551	if ( (aphidiff==1 && aetadiff==0) \|\|
	552	(aphidiff==0 && aetadiff==1) ) {
	553	crossEnergy += cells->GetCellAmplitude(cellAbsId);
	554	}
	555	}
	556
	557	return crossEnergy;
	558	}
	559
30159e6f	560	//________________________________________________________________________
	561	Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
	562	{
	563	// Get maximum energy of attached cell.
	564
	565	id = -1;
	566
	567	AliVCaloCells *cells = 0;
	568	cells = fESD->GetEMCALCells();
	569	if (!cells)
	570	return 0;
	571
	572	Double_t maxe = 0;
	573	Int_t ncells = cluster->GetNCells();
	574	for (Int_t i=0; i<ncells; i++) {
	575	Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
	576	if (e>maxe) {
	577	maxe = e;
	578	id = cluster->GetCellAbsId(i);
	579	}
	580	}
	581	return maxe;
	582	}
bd092f0f	583
30159e6f	584	//________________________________________________________________________
	585	void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
	586	{
bd092f0f	587	// Called once at the end of the query.
30159e6f	588	}