[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 <TH3F.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 "AliESDtrack.h"
#include "AliESDtrackCuts.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),
  fIsMc(0),
  fDebug(0),
  fExoticCut(0.97),
  fIsoConeR(0.4),
  fNDimensions(7),
  fECut(3.),
  fTrackMult(0),        
  fESD(0),
  fMCEvent(0),
  fStack(0),
  fOutputList(0),
  fEvtSel(0),
  fNClusEt10(0),
  fPVtxZ(0),  
  fMCDirPhotonPtEtaPhi(0),
  fDecayPhotonPtMC(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),
  fIsMc(0),
  fDebug(0),
  fExoticCut(0.97),
  fIsoConeR(0.4),
  fNDimensions(7),
  fECut(3.),
  fTrackMult(0),        
  fESD(0),
  fMCEvent(0),
  fStack(0),
  fOutputList(0),
  fEvtSel(0),
  fNClusEt10(0),
  fPVtxZ(0),            
  fMCDirPhotonPtEtaPhi(0),
  fDecayPhotonPtMC(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);
  
  fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5);
  fOutputList->Add(fNClusEt10);
  
  fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
  fOutputList->Add(fPVtxZ);

  fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",1000,0,100,154,-0.77,0.77,130,1.38,3.20);
  fMCDirPhotonPtEtaPhi->Sumw2();
  fOutputList->Add(fMCDirPhotonPtEtaPhi);

  fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",1000,0,100);
  fDecayPhotonPtMC->Sumw2();
  fOutputList->Add(fDecayPhotonPtMC);

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

  Int_t nEt=1000, nM02=400, nCeIso=1000, nTrIso=1000,  nAllIso=1000,  nCeIsoNoUE=1000,  nAllIsoNoUE=1000, nTrClDphi=200, nTrClDeta=100, nClEta=140, nClPhi=128, nTime=60, nMult=20;
  Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult};
  fNDimensions = sizeof(bins)/sizeof(Int_t);
  const Int_t ndims =   fNDimensions;
  Double_t xmin[] = { 0.,   0.,  -10.,   -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5};
  Double_t xmax[] = { 100., 4., 190., 190., 190.,  190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5};
  fHnOutput =  new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, CeIsoNoUESub, AllIsoNoUESub, d#phi_{trk},d#eta_{trk},#eta_{clus},#phi_{clus},T_{max},mult", 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(fIsMc)
    isSelected = kTRUE;
  if(fDebug)
    printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc);
  if(!isSelected )
        return; 

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

  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);
    }
  }
  AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
  fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);  

  fESD->GetEMCALClusters(fCaloClusters);
  fEMCalCells = fESD->GetEMCALCells();
  
    
  FillClusHists(); 
  if(fDebug)
    printf("passed calling of FillClusHists\n");

  fCaloClusters->Clear();
  fSelPrimTracks->Clear();

  fMCEvent = MCEvent();
  if(fMCEvent){
    if(fDebug)
      std::cout<<"MCevent exists\n";
    fStack = (AliStack*)fMCEvent->Stack();
  }
  else{
    if(fDebug)
      std::cout<<"ERROR: NO MC EVENT!!!!!!\n";
  }
  FillMcHists();
  if(fDebug)
    printf("passed calling of FillMcHists\n");


  PostData(1, fOutputList);
}      

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

  if(!fCaloClusters)
    return;
  const Int_t nclus = fCaloClusters->GetEntries();
  if(nclus==0)
    return;
  if(fDebug)
    printf("Inside FillClusHists and there are %d clusters\n",nclus);
  Double_t maxE = 0;
  Int_t nclus10 = 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());
    if(Et>10)
      nclus10++;
    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] = ceiso-Et;
    outputValues[6] = alliso-Et;
    outputValues[7] = c->GetTrackDx();
    outputValues[8] = c->GetTrackDz();
    outputValues[9] = clsVec.Eta();
    outputValues[10] = clsVec.Phi();
    outputValues[11] = fEMCalCells->GetCellTime(id);
    outputValues[12] = fTrackMult;
    fHnOutput->Fill(outputValues);
    if(c->E()>maxE)
      maxE = c->E();
  }
  fNClusHighClusE->Fill(maxE,nclus);
  fNClusEt10->Fill(nclus10);
} 

//________________________________________________________________________
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 ::FillMcHists()
{
  if(!fStack)
    return;
  Int_t nTracks = fStack->GetNtrack();
  if(fDebug)
    printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
  for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
    TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack));  
    if(!mcp)
      continue;  
    Int_t pdg = mcp->GetPdgCode();
    if(pdg!=22)
      continue;
    if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2)
      continue;
    Int_t imom = mcp->GetMother(0);
    if(imom<0 || imom>nTracks)
      continue;
    TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom));  
    if(!mcmom)
      continue;
    Int_t pdgMom = mcmom->GetPdgCode();
    if((imom==6 || imom==7) && pdgMom==22) {
      fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
      if(fDebug){
	printf("Found \"photonic\" parton at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d,\n",imom,mcmom->Pt(), mcmom->Eta(), mcmom->Phi(), mcmom->GetStatusCode());
	printf("with a final photon at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d\n",iTrack,mcp->Pt(), mcp->Eta(), mcp->Phi(),mcp->GetStatusCode());
      }
    }
    else{
      if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
	fDecayPhotonPtMC->Fill(mcp->Pt());
    }
  }
}
//________________________________________________________________________
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>
caaf99d3	10	#include <TH3F.h>
965c985f	11	#include <THnSparse.h>
bd092f0f	12	#include <TLorentzVector.h>
a62631a9	13	#include <TList.h>
bd092f0f	14
	15	#include "AliAnalysisManager.h"
	16	#include "AliAnalysisTask.h"
	17	#include "AliEMCALGeometry.h"
	18	#include "AliEMCALRecoUtils.h"
	19	#include "AliESDCaloCells.h"
	20	#include "AliESDCaloCluster.h"
30159e6f	21	#include "AliESDEvent.h"
30159e6f	22	#include "AliESDHeader.h"
30159e6f	23	#include "AliESDInputHandler.h"
bd092f0f	24	#include "AliESDUtils.h"
30159e6f	25	#include "AliESDtrack.h"
30159e6f	26	#include "AliESDtrackCuts.h"
bd092f0f	27	#include "AliMCEvent.h"
	28	#include "AliMCEventHandler.h"
	29	#include "AliStack.h"
30159e6f	30	#include "AliV0vertexer.h"
30159e6f	31	#include "AliVCluster.h"
30159e6f	32
30159e6f	33	ClassImp(AliAnalysisTaskEMCALIsoPhoton)
	34
	35	//________________________________________________________________________
bd092f0f	36	AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() :
	37	AliAnalysisTaskSE(),
	38	fCaloClusters(0),
	39	fSelPrimTracks(0),
3f4073ba	40	fTracks(0),
bd092f0f	41	fEMCalCells(0),
	42	fPrTrCuts(0),
	43	fGeom(0x0),
9148fa89	44	fGeoName("EMCAL_COMPLETEV1"),
	45	fPeriod("LHC11c"),
	46	fTrigBit("kEMC7"),
	47	fIsTrain(0),
c7bb0b43	48	fIsMc(0),
ecd47673	49	fDebug(0),
9148fa89	50	fExoticCut(0.97),
9148fa89	51	fIsoConeR(0.4),
965c985f	52	fNDimensions(7),
965c985f	53	fECut(3.),
16a4050e	54	fTrackMult(0),
bd092f0f	55	fESD(0),
c7bb0b43	56	fMCEvent(0),
c7bb0b43	57	fStack(0),
bd092f0f	58	fOutputList(0),
bd092f0f	59	fEvtSel(0),
0b21f686	60	fNClusEt10(0),
c7bb0b43	61	fPVtxZ(0),
caaf99d3	62	fMCDirPhotonPtEtaPhi(0),
c7bb0b43	63	fDecayPhotonPtMC(0),
bd092f0f	64	fCellAbsIdVsAmpl(0),
16a4050e	65	fNClusHighClusE(0),
965c985f	66	fHnOutput(0)
bd092f0f	67	{
	68	// Default constructor.
	69	}
	70
	71	//________________________________________________________________________
	72	AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) :
	73	AliAnalysisTaskSE(name),
30159e6f	74	fCaloClusters(0),
30159e6f	75	fSelPrimTracks(0),
3f4073ba	76	fTracks(0),
30159e6f	77	fEMCalCells(0),
	78	fPrTrCuts(0),
	79	fGeom(0x0),
	80	fGeoName("EMCAL_COMPLETEV1"),
	81	fPeriod("LHC11c"),
751194e8	82	fTrigBit("kEMC7"),
30159e6f	83	fIsTrain(0),
c7bb0b43	84	fIsMc(0),
ecd47673	85	fDebug(0),
30159e6f	86	fExoticCut(0.97),
30159e6f	87	fIsoConeR(0.4),
965c985f	88	fNDimensions(7),
965c985f	89	fECut(3.),
16a4050e	90	fTrackMult(0),
30159e6f	91	fESD(0),
c7bb0b43	92	fMCEvent(0),
c7bb0b43	93	fStack(0),
30159e6f	94	fOutputList(0),
30159e6f	95	fEvtSel(0),
0b21f686	96	fNClusEt10(0),
bd092f0f	97	fPVtxZ(0),
caaf99d3	98	fMCDirPhotonPtEtaPhi(0),
c7bb0b43	99	fDecayPhotonPtMC(0),
c7bb0b43	100	fCellAbsIdVsAmpl(0),
bd092f0f	101	fNClusHighClusE(0),
965c985f	102	fHnOutput(0)
30159e6f	103	{
	104	// Constructor
	105
	106	// Define input and output slots here
	107	// Input slot #0 works with a TChain
	108	DefineInput(0, TChain::Class());
	109	// Output slot #0 id reserved by the base class for AOD
	110	// Output slot #1 writes into a TH1 container
	111	DefineOutput(1, TList::Class());
	112	}
bd092f0f	113
30159e6f	114	//________________________________________________________________________
	115	void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects()
	116	{
bd092f0f	117	// Create histograms, called once.
30159e6f	118
	119	fCaloClusters = new TRefArray();
	120	fSelPrimTracks = new TObjArray();
	121
	122
	123	fOutputList = new TList();
a62631a9	124	fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
30159e6f	125
	126	fGeom = AliEMCALGeometry::GetInstance(fGeoName);
	127
	128	fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2);
	129	fOutputList->Add(fEvtSel);
0b21f686	130
	131	fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5);
	132	fOutputList->Add(fNClusEt10);
30159e6f	133
	134	fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100);
	135	fOutputList->Add(fPVtxZ);
c7bb0b43	136
caaf99d3	137	fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",1000,0,100,154,-0.77,0.77,130,1.38,3.20);
	138	fMCDirPhotonPtEtaPhi->Sumw2();
	139	fOutputList->Add(fMCDirPhotonPtEtaPhi);
c7bb0b43	140
	141	fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",1000,0,100);
	142	fDecayPhotonPtMC->Sumw2();
	143	fOutputList->Add(fDecayPhotonPtMC);
	144
30159e6f	145	fCellAbsIdVsAmpl = new TH2F("hCellAbsIdVsAmpl","cell abs id vs cell amplitude (energy);E (GeV);ID",200,0,100,244810,-0.5,244810-0.5);
a62631a9	146	fOutputList->Add(fCellAbsIdVsAmpl);
30159e6f	147
	148	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);
	149	fOutputList->Add(fNClusHighClusE);
	150
302b398c	151	Int_t nEt=1000, nM02=400, nCeIso=1000, nTrIso=1000, nAllIso=1000, nCeIsoNoUE=1000, nAllIsoNoUE=1000, nTrClDphi=200, nTrClDeta=100, nClEta=140, nClPhi=128, nTime=60, nMult=20;
302b398c	152	Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult};
c1f18270	153	fNDimensions = sizeof(bins)/sizeof(Int_t);
c1f18270	154	const Int_t ndims = fNDimensions;
302b398c	155	Double_t xmin[] = { 0., 0., -10., -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5};
	156	Double_t xmax[] = { 100., 4., 190., 190., 190., 190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5};
	157	fHnOutput = new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, CeIsoNoUESub, AllIsoNoUESub, d#phi_{trk},d#eta_{trk},#eta_{clus},#phi_{clus},T_{max},mult", ndims, bins, xmin, xmax);
965c985f	158	fOutputList->Add(fHnOutput);
	159
	160
	161
30159e6f	162	PostData(1, fOutputList);
	163	}
	164
	165	//________________________________________________________________________
	166	void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *)
	167	{
bd092f0f	168	// Main loop, called for each event.
	169
	170	// event trigger selection
30159e6f	171	Bool_t isSelected = 0;
751194e8	172	if(fPeriod.Contains("11a")){
	173	if(fTrigBit.Contains("kEMC"))
	174	isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
	175	else
	176	isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
	177	}
	178	else{
	179	if(fTrigBit.Contains("kEMC"))
	180	isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
	181	else
	182	isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7);
	183	}
c7bb0b43	184	if(fIsMc)
c7bb0b43	185	isSelected = kTRUE;
ecd47673	186	if(fDebug)
ecd47673	187	printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc);
30159e6f	188	if(!isSelected )
30159e6f	189	return;
30159e6f	190
30159e6f	191	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
	192	if (!fESD) {
	193	printf("ERROR: fESD not available\n");
	194	return;
	195	}
	196
	197	fEvtSel->Fill(0);
ecd47673	198	if(fDebug)
ecd47673	199	printf("fESD is ok\n");
30159e6f	200
	201	AliESDVertex pv = (AliESDVertex)fESD->GetPrimaryVertex();
	202	if(!pv)
	203	return;
	204	fPVtxZ->Fill(pv->GetZ());
	205	if(TMath::Abs(pv->GetZ())>15)
	206	return;
ecd47673	207	if(fDebug)
ecd47673	208	printf("passed vertex cut\n");
30159e6f	209
30159e6f	210	fEvtSel->Fill(1);
bd092f0f	211
3f4073ba	212	if (!fTracks)
3f4073ba	213	fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks"));
30159e6f	214	// Track loop to fill a pT spectrum
3f4073ba	215	const Int_t Ntracks = fTracks->GetEntriesFast();
	216	for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
	217	// for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
	218	//AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks);
	219	AliVTrack track = static_cast<AliVTrack>(fTracks->At(iTracks));
30159e6f	220	if (!track)
	221	continue;
	222	if (fPrTrCuts && fPrTrCuts->IsSelected(track)){
	223	fSelPrimTracks->Add(track);
	224	//printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi());
	225	}
bd092f0f	226	}
bd092f0f	227
30159e6f	228	if(!fIsTrain){
	229	for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
	230	if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
	231	break;
	232	fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
	233	}
	234	}
16a4050e	235	AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
	236	fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
	237
30159e6f	238	fESD->GetEMCALClusters(fCaloClusters);
	239	fEMCalCells = fESD->GetEMCALCells();
	240
	241
	242	FillClusHists();
ecd47673	243	if(fDebug)
ecd47673	244	printf("passed calling of FillClusHists\n");
30159e6f	245
	246	fCaloClusters->Clear();
	247	fSelPrimTracks->Clear();
c7bb0b43	248
c7bb0b43	249	fMCEvent = MCEvent();
ecd47673	250	if(fMCEvent){
ecd47673	251	if(fDebug)
ec9ab86b	252	std::cout<<"MCevent exists\n";
c7bb0b43	253	fStack = (AliStack*)fMCEvent->Stack();
ecd47673	254	}
	255	else{
	256	if(fDebug)
ec9ab86b	257	std::cout<<"ERROR: NO MC EVENT!!!!!!\n";
ecd47673	258	}
63e40cb8	259	FillMcHists();
ecd47673	260	if(fDebug)
	261	printf("passed calling of FillMcHists\n");
	262
c7bb0b43	263
30159e6f	264	PostData(1, fOutputList);
30159e6f	265	}
bd092f0f	266
30159e6f	267	//________________________________________________________________________
	268	void AliAnalysisTaskEMCALIsoPhoton::FillClusHists()
	269	{
bd092f0f	270	// Fill cluster histograms.
bd092f0f	271
30159e6f	272	if(!fCaloClusters)
	273	return;
	274	const Int_t nclus = fCaloClusters->GetEntries();
	275	if(nclus==0)
	276	return;
ecd47673	277	if(fDebug)
ecd47673	278	printf("Inside FillClusHists and there are %d clusters\n",nclus);
e681d9ce	279	Double_t maxE = 0;
0b21f686	280	Int_t nclus10 = 0;
30159e6f	281	for(Int_t ic=0;ic<nclus;ic++){
	282	maxE=0;
	283	AliESDCaloCluster c = static_cast<AliESDCaloCluster>(fCaloClusters->At(ic));
	284	if(!c)
	285	continue;
	286	if(!c->IsEMCAL())
	287	continue;
965c985f	288	if(c->E()<fECut)
965c985f	289	continue;
30159e6f	290	Short_t id;
	291	Double_t Emax = GetMaxCellEnergy( c, id);
	292	Double_t Ecross = GetCrossEnergy( c, id);
	293	if((1-Ecross/Emax)>fExoticCut)
	294	continue;
	295	Float_t clsPos[3] = {0,0,0};
	296	c->GetPosition(clsPos);
	297	TVector3 clsVec(clsPos);
3dc2825e	298	Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
0b21f686	299	if(Et>10)
0b21f686	300	nclus10++;
30159e6f	301	Float_t ceiso, cephiband, cecore;
	302	Float_t triso, trphiband, trcore;
	303	Float_t alliso, allphiband, allcore;
	304	Float_t phibandArea = (1.4 - 2fIsoConeR)2*fIsoConeR;
	305	Float_t netConeArea = TMath::Pi()(fIsoConeRfIsoConeR - 0.04*0.04);
	306	GetCeIso(clsVec, ceiso, cephiband, cecore);
	307	GetTrIso(clsVec, triso, trphiband, trcore);
	308	alliso = ceiso + triso;
	309	allphiband = cephiband + trphiband;
	310	allcore = cecore + trcore;
	311	Float_t ceisoue = cephiband/phibandArea*netConeArea;
	312	Float_t trisoue = trphiband/phibandArea*netConeArea;
	313	Float_t allisoue = allphiband/phibandArea*netConeArea;
965c985f	314	const Int_t ndims = fNDimensions;
	315	Double_t outputValues[ndims];
	316	outputValues[0] = Et;
	317	outputValues[1] = c->GetM02();
	318	outputValues[2] = ceiso-cecore-ceisoue;
	319	outputValues[3] = triso-trisoue;
	320	outputValues[4] = alliso-cecore-allisoue;
302b398c	321	outputValues[5] = ceiso-Et;
302b398c	322	outputValues[6] = alliso-Et;
717bb45b	323	outputValues[7] = c->GetTrackDx();
	324	outputValues[8] = c->GetTrackDz();
	325	outputValues[9] = clsVec.Eta();
	326	outputValues[10] = clsVec.Phi();
16a4050e	327	outputValues[11] = fEMCalCells->GetCellTime(id);
16a4050e	328	outputValues[12] = fTrackMult;
965c985f	329	fHnOutput->Fill(outputValues);
30159e6f	330	if(c->E()>maxE)
	331	maxE = c->E();
	332	}
	333	fNClusHighClusE->Fill(maxE,nclus);
0b21f686	334	fNClusEt10->Fill(nclus10);
30159e6f	335	}
bd092f0f	336
30159e6f	337	//________________________________________________________________________
	338	void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
	339	{
bd092f0f	340	// Get cell isolation.
bd092f0f	341
30159e6f	342	if(!fEMCalCells)
	343	return;
	344	const Int_t ncells = fEMCalCells->GetNumberOfCells();
	345	Float_t totiso=0;
	346	Float_t totband=0;
	347	Float_t totcore=0;
	348	Float_t etacl = vec.Eta();
	349	Float_t phicl = vec.Phi();
	350	Float_t thetacl = vec.Theta();
	351	if(phicl<0)
	352	phicl+=TMath::TwoPi();
	353	Int_t absid = -1;
	354	Float_t eta=-1, phi=-1;
	355	for(int icell=0;icell<ncells;icell++){
	356	absid = TMath::Abs(fEMCalCells->GetCellNumber(icell));
	357	if(!fGeom)
	358	return;
	359	fGeom->EtaPhiFromIndex(absid,eta,phi);
	360	Float_t dphi = TMath::Abs(phi-phicl);
	361	Float_t deta = TMath::Abs(eta-etacl);
	362	Float_t R = TMath::Sqrt(detadeta + dphidphi);
	363	Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
	364	if(R<fIsoConeR){
	365	totiso += etcell;
	366	if(R<0.04)
	367	totcore += etcell;
	368	}
	369	else{
	370	if(dphi>fIsoConeR)
	371	continue;
	372	if(deta<fIsoConeR)
	373	continue;
	374	totband += fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl);
	375	}
	376	}
	377	iso = totiso;
	378	phiband = totband;
	379	core = totcore;
	380	}
bd092f0f	381
30159e6f	382	//________________________________________________________________________
	383	void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core)
	384	{
bd092f0f	385	// Get track isolation.
bd092f0f	386
30159e6f	387	if(!fSelPrimTracks)
	388	return;
	389	const Int_t ntracks = fSelPrimTracks->GetEntries();
	390	Double_t totiso=0;
	391	Double_t totband=0;
	392	Double_t totcore=0;
	393	Double_t etacl = vec.Eta();
	394	Double_t phicl = vec.Phi();
	395	if(phicl<0)
	396	phicl+=TMath::TwoPi();
	397	for(int itrack=0;itrack<ntracks;itrack++){
3f4073ba	398	AliVTrack track = static_cast<AliVTrack> (fSelPrimTracks->At(itrack));
30159e6f	399	if(!track)
	400	continue;
	401	Double_t dphi = TMath::Abs(track->Phi()-phicl);
	402	Double_t deta = TMath::Abs(track->Eta()-etacl);
	403	Double_t R = TMath::Sqrt(detadeta + dphidphi);
	404	Double_t pt = track->Pt();
	405	if(R<fIsoConeR){
	406	totiso += track->Pt();
	407	if(R<0.04)
	408	totcore += pt;
	409	}
	410	else{
	411	if(dphi>fIsoConeR)
	412	continue;
	413	if(deta<fIsoConeR)
	414	continue;
	415	totband += track->Pt();
	416	}
	417	}
	418	iso = totiso;
	419	phiband = totband;
	420	core = totcore;
	421	}
bd092f0f	422
30159e6f	423	//________________________________________________________________________
	424	Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
	425	{
	426	// Calculate the energy of cross cells around the leading cell.
	427
	428	AliVCaloCells *cells = 0;
	429	cells = fESD->GetEMCALCells();
	430	if (!cells)
	431	return 0;
	432
30159e6f	433	if (!fGeom)
	434	return 0;
	435
	436	Int_t iSupMod = -1;
	437	Int_t iTower = -1;
	438	Int_t iIphi = -1;
	439	Int_t iIeta = -1;
	440	Int_t iphi = -1;
	441	Int_t ieta = -1;
	442	Int_t iphis = -1;
	443	Int_t ietas = -1;
	444
	445	Double_t crossEnergy = 0;
	446
	447	fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
	448	fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
	449
	450	Int_t ncells = cluster->GetNCells();
	451	for (Int_t i=0; i<ncells; i++) {
	452	Int_t cellAbsId = cluster->GetCellAbsId(i);
	453	fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
	454	fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
	455	Int_t aphidiff = TMath::Abs(iphi-iphis);
	456	if (aphidiff>1)
	457	continue;
	458	Int_t aetadiff = TMath::Abs(ieta-ietas);
	459	if (aetadiff>1)
	460	continue;
	461	if ( (aphidiff==1 && aetadiff==0) \|\|
	462	(aphidiff==0 && aetadiff==1) ) {
	463	crossEnergy += cells->GetCellAmplitude(cellAbsId);
	464	}
	465	}
	466
	467	return crossEnergy;
	468	}
	469
30159e6f	470	//________________________________________________________________________
	471	Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
	472	{
	473	// Get maximum energy of attached cell.
	474
	475	id = -1;
	476
	477	AliVCaloCells *cells = 0;
	478	cells = fESD->GetEMCALCells();
	479	if (!cells)
	480	return 0;
	481
	482	Double_t maxe = 0;
	483	Int_t ncells = cluster->GetNCells();
	484	for (Int_t i=0; i<ncells; i++) {
	485	Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
	486	if (e>maxe) {
	487	maxe = e;
	488	id = cluster->GetCellAbsId(i);
	489	}
	490	}
	491	return maxe;
	492	}
bd092f0f	493
c7bb0b43	494	//________________________________________________________________________
	495	void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists()
	496	{
	497	if(!fStack)
	498	return;
	499	Int_t nTracks = fStack->GetNtrack();
ecd47673	500	if(fDebug)
ecd47673	501	printf("Inside FillMcHists and there are %d mcparts\n",nTracks);
c7bb0b43	502	for(Int_t iTrack=0;iTrack<nTracks;iTrack++){
	503	TParticle mcp = static_cast<TParticle>(fStack->Particle(iTrack));
	504	if(!mcp)
	505	continue;
	506	Int_t pdg = mcp->GetPdgCode();
	507	if(pdg!=22)
	508	continue;
f2acdbe9	509	if(TMath::Abs(mcp->Eta())>0.7 \|\|mcp->Phi()<1.4 \|\| mcp->Phi()>3.2)
f2acdbe9	510	continue;
c7bb0b43	511	Int_t imom = mcp->GetMother(0);
	512	if(imom<0 \|\| imom>nTracks)
	513	continue;
	514	TParticle mcmom = static_cast<TParticle>(fStack->Particle(imom));
	515	if(!mcmom)
	516	continue;
	517	Int_t pdgMom = mcmom->GetPdgCode();
ec9ab86b	518	if((imom==6 \|\| imom==7) && pdgMom==22) {
caaf99d3	519	fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi());
40c4f1bf	520	if(fDebug){
f2acdbe9	521	printf("Found \"photonic\" parton at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d,\n",imom,mcmom->Pt(), mcmom->Eta(), mcmom->Phi(), mcmom->GetStatusCode());
40c4f1bf	522	printf("with a final photon at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d\n",iTrack,mcp->Pt(), mcp->Eta(), mcp->Phi(),mcp->GetStatusCode());
40c4f1bf	523	}
f2acdbe9	524	}
c7bb0b43	525	else{
	526	if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000)
	527	fDecayPhotonPtMC->Fill(mcp->Pt());
	528	}
	529	}
	530	}
30159e6f	531	//________________________________________________________________________
	532	void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *)
	533	{
bd092f0f	534	// Called once at the end of the query.
30159e6f	535	}