From 1c662fe8739189662fd24f92e47bd456f3e1ae46 Mon Sep 17 00:00:00 2001 From: gconesab Date: Thu, 26 Jun 2014 21:49:21 +0200 Subject: [PATCH] add new reference task for photon isolation studies --- PWGGA/CMakelibPWGGAEMCALTasks.pkg | 1 + .../AliAnalysisTaskEMCALPhotonIsolation.cxx | 1342 +++++++++++++++++ .../AliAnalysisTaskEMCALPhotonIsolation.h | 187 +++ .../macros/AddTaskEMCALPhotonIsolation.C | 81 + .../macros/AddTaskPhotonPreparation.C | 99 ++ PWGGA/PWGGAEMCALTasksLinkDef.h | 1 + 6 files changed, 1711 insertions(+) create mode 100644 PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhotonIsolation.cxx create mode 100644 PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhotonIsolation.h create mode 100644 PWGGA/EMCALTasks/macros/AddTaskEMCALPhotonIsolation.C create mode 100644 PWGGA/EMCALTasks/macros/AddTaskPhotonPreparation.C diff --git a/PWGGA/CMakelibPWGGAEMCALTasks.pkg b/PWGGA/CMakelibPWGGAEMCALTasks.pkg index 5b7253481b7..92f6df84e1f 100644 --- a/PWGGA/CMakelibPWGGAEMCALTasks.pkg +++ b/PWGGA/CMakelibPWGGAEMCALTasks.pkg @@ -40,6 +40,7 @@ set ( SRCS EMCALTasks/AliAnalysisTaskEMCALPi0V2ShSh.cxx EMCALTasks/AliEMCalpi0ClusterEvaluationTask.cxx EMCALTasks/AliAnalysisTaskEMCALIsolation.cxx + EMCALTasks/AliAnalysisTaskEMCALPhotonIsolation.cxx EMCALTasks/AliAnalysisTaskEMCALMesonGGSDM.cxx EMCALTasks/AliAnalysisTaskEMCALMesonGGSDMpPb.cxx EMCALTasks/AliAnalysisTaskSDMGammaMC.cxx diff --git a/PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhotonIsolation.cxx b/PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhotonIsolation.cxx new file mode 100644 index 00000000000..7be3bb21c20 --- /dev/null +++ b/PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhotonIsolation.cxx @@ -0,0 +1,1342 @@ +// $Id$ +// +// Emcal Neutral Cluster analysis base task. +// +// Authors: D.Lodato,L.Ronflette, M.Marquard + + + +#include +#include +#include +#include +#include +#include +#include +#include +#include "AliAnalysisManager.h" +#include "AliCentrality.h" +#include "AliEMCALGeometry.h" +#include "AliESDEvent.h" +#include "AliAODEvent.h" +#include "AliLog.h" +#include "AliVCluster.h" +#include "AliVEventHandler.h" +#include "AliVParticle.h" +#include "AliClusterContainer.h" +#include "AliVTrack.h" +#include "AliEmcalParticle.h" +#include "AliParticleContainer.h" +#include "AliAODCaloCluster.h" +#include "AliESDCaloCluster.h" +#include "AliVCaloCells.h" +#include "AliPicoTrack.h" + +#include "AliAnalysisTaskEMCALPhotonIsolation.h" + + +ClassImp(AliAnalysisTaskEMCALPhotonIsolation) + +//________________________________________________________________________ +AliAnalysisTaskEMCALPhotonIsolation::AliAnalysisTaskEMCALPhotonIsolation() : +AliAnalysisTaskEmcal("AliAnalysisTaskEMCALPhotonIsolation",kTRUE), +//fParticleCollArray(), +fNCluster(0), +fWho(-1), +//fOutputList(0), +fTrackMult(0), +fTrackMultEMCAL(0), +fClustMult(0), +fPVZBefore(0), +fEtaPhiCell(0), +fEtaPhiClus(0), +fClusEvsClusT(0), +fGoodEventsOnPVZ(0), +fPT(0), +fM02(0), +fNLM(0), +fDeltaETAClusTrackVSpT(0), +fDeltaPHIClusTrackVSpT(0), +fEtIsoCells(0), +fEtIsoClust(0), +fPtIsoTrack(0), +fPtEtIsoTC(0), +fPhiBandUEClust(0), +fEtaBandUEClust(0), +fPhiBandUECells(0), +fEtaBandUECells(0), +fPhiBandUETracks(0), +fEtaBandUETracks(0), +fPerpConesUETracks(0), +fTPCWithoutIsoConeB2BbandUE(0), +fNTotClus10GeV(0), +fRecoPV(0), +fEtIsolatedCells(0), +fEtIsolatedClust(0), +fEtIsolatedTracks(0), +fTest(0), +fOutputTHnS(0), +fOutPTMAX(0), +fOutputTree(0), +fIsoConeRadius(0.4), +fEtIsoMethod(0), +fEtIsoThreshold(4), +fdetacut(0.025), +fdphicut(0.03), +fM02mincut(0.1), +fM02maxcut(0.3), +fQA(0), +fIsMC(0), +fTPC4Iso(0), +fIsoMethod(0), +fUEMethod(0), +fVertex(0), +fNDimensions(0), +fisLCAnalysis(0), +fevents(0), +flambda0T(0), +fEtT(0), +fPtT(0), +fEtisoT(0), +fPtisoT(0), +fetaT(0), +fphiT(0), +fsumEtisoconeT(0), +fsumEtUE(0) +//Tracks(0), +//clusters(0) + +{ + // Default constructor. + + //fParticleCollArray.SetOwner(kTRUE); + for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0; + + SetMakeGeneralHistograms(kTRUE); +} + +//________________________________________________________________________ +AliAnalysisTaskEMCALPhotonIsolation::AliAnalysisTaskEMCALPhotonIsolation(const char *name, Bool_t histo) : +AliAnalysisTaskEmcal(name, histo), +//fParticleCollArray(), +fNCluster(0), +fWho(-1), +//fOutputList(0), +fTrackMult(0), +fTrackMultEMCAL(0), +fClustMult(0), +fPVZBefore(0), +fEtaPhiCell(0), +fEtaPhiClus(0), +fClusEvsClusT(0), +fGoodEventsOnPVZ(0), +fPT(0), +fM02(0), +fNLM(0), +fDeltaETAClusTrackVSpT(0), +fDeltaPHIClusTrackVSpT(0), +fEtIsoCells(0), +fEtIsoClust(0), +fPtIsoTrack(0), +fPtEtIsoTC(0), +fPhiBandUEClust(0), +fEtaBandUEClust(0), +fPhiBandUECells(0), +fEtaBandUECells(0), +fPhiBandUETracks(0), +fEtaBandUETracks(0), +fPerpConesUETracks(0), +fTPCWithoutIsoConeB2BbandUE(0), +fNTotClus10GeV(0), +fRecoPV(0), +fEtIsolatedCells(0), +fEtIsolatedClust(0), +fEtIsolatedTracks(0), +fTest(0), +fOutputTHnS(0), +fOutPTMAX(0), +fOutputTree(0), +fIsoConeRadius(0.4), +fEtIsoMethod(0), +fEtIsoThreshold(4), +fdetacut(0.025), +fdphicut(0.03), +fM02mincut(0.1), +fM02maxcut(0.3), +fQA(0), +fIsMC(0), +fTPC4Iso(0), +fIsoMethod(0), +fUEMethod(0), +fVertex(0), +fNDimensions(0), +fisLCAnalysis(0), +fevents(0), +flambda0T(0), +fEtT(0), +fPtT(0), +fEtisoT(0), +fPtisoT(0), +fetaT(0), +fphiT(0), +fsumEtisoconeT(0), +fsumEtUE(0) +//Tracks(0), +//clusters(0) + +{ + // Standard constructor. + + //fParticleCollArray.SetOwner(kTRUE); + for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0; + + SetMakeGeneralHistograms(kTRUE); +} + +//________________________________________________________________________ +AliAnalysisTaskEMCALPhotonIsolation::~AliAnalysisTaskEMCALPhotonIsolation(){ + // Destructor +} + + +//________________________________________________________________________ +void AliAnalysisTaskEMCALPhotonIsolation::UserCreateOutputObjects(){ + + AliAnalysisTaskEmcal::UserCreateOutputObjects(); + + + if ((fIsoMethod == 0 || fIsoMethod == 1) && fTPC4Iso){ + cout<<"Error: Iso_Methods with CELLS and CLUSTERS work only within EMCAL "< 1) { + cout<<"Error: UE_Methods with CELLS and CLUSTERS work only within EMCAL"<1 || fWho==-1){ + cout<<"Error!!! OutputMode Can Only Be 0: TTree; 1: THnSparse"<SetOwner(); + //Initialize the common Output histograms + switch (fWho) + { + case 0: + //Initialization by Lucile and Marco + + fOutputTree = new TTree("OutTree","OutTree"); + fOutputTree->Branch("fevents",&fevents); + fOutputTree->Branch("flambda0T",&flambda0T); + fOutputTree->Branch("fEtT",&fEtT); + fOutputTree->Branch("fPtT",&fPtT); + fOutputTree->Branch("fEtisoT",&fEtisoT); + fOutputTree->Branch("fPtTiso",&fPtisoT); + fOutputTree->Branch("fetaT",&fetaT); + fOutputTree->Branch("fphiT",&fphiT); + fOutputTree->Branch("fsumEtisoconeT",&fsumEtisoconeT); + fOutputTree->Branch("fsumEtUE",&fsumEtUE); + + fOutput->Add(fOutputTree); + + break; + case 1: + //Initialization by Davide; + + TString sTitle; + Int_t binTrackMult=100, binPT=70, binM02=100, binETiso=100, binETUE=110, binETisoUE=110, binetacl=140,binphicl=100, binPTMC=70; + Int_t bins[] = {binTrackMult, binPT, binM02, binETiso, binETUE, binETisoUE, binetacl, binphicl, binPTMC, binPTMC}; + + fNDimensions = sizeof(bins)/sizeof(Int_t); + const Int_t ndims = fNDimensions; + + Double_t xmin[]= { 0., 0., 0., -10., -10., -10.,-1.0, 1. ,-10.,-10.}; + + Double_t xmax[]= {1000., 70., 2., 100., 100., 100., 1.0, 3.5, 60., 60.}; + + sTitle = Form("Direct Photons: Track Multiplicity, p_{T} , M02 , E_{T} Iso%s in %s, E_{T} UE %s in %s, E_{T} Iso_%s - E_{T} UE_%s in %s, #eta_{clus} distribution,#phi_{clus} distribution,MC_pT,MC_pT_incone; N_{ch}; p_{T} (GeV/c); M02; E_{T}^{iso%s} (GeV/c) ; E_{T}^{UE%s} (GeV/c); E_{T}^{iso%s}-E_{T}^{UE%s} (GeV/c); #eta_{cl}; #phi_{cl}; MC p_{T}; MC p_{T}^{incone}", sIsoMethod.Data(), sBoundaries.Data(), sUEMethod.Data(), sBoundaries.Data(), sIsoMethod.Data(), sUEMethod.Data(), sBoundaries.Data(), sIsoMethod.Data(), sUEMethod.Data(), sIsoMethod.Data(), sUEMethod.Data()); + + fOutputTHnS = new THnSparseF("fHnOutput",sTitle.Data(), ndims, bins, xmin, xmax); + + fOutput->Add(fOutputTHnS); + + Int_t binsbis[] = {binPT, binM02, binETisoUE}; + Double_t xminbis[] = { 0., 0., -10.}; + Double_t xmaxbis[] = {100., 2., 100.}; + + fOutPTMAX = new THnSparseF ("fOutPTMAX","3D matrix E_{#gamma} VS M02 VS pT_{max}^{cone}; E_{T}^{#gamma} (GeV/c); M02; p_{T}^{Iso}(GeV/c)",3,binsbis,xminbis,xmaxbis); + + fOutput->Add(fOutPTMAX); + break; + } + } + if(fQA){ + //Include QA plots to the OutputList //DEFINE BETTER THE BINNING AND THE AXES LIMITS + fTrackMult = new TH1D ("hTrackMult","Tracks multiplicity Distribution",250,0.,1000.); + fTrackMult->Sumw2(); + fOutput->Add(fTrackMult); + + fTrackMultEMCAL = new TH1D ("hTrackMultEMCAL","Tracks multiplicity Distribution inside EMCAL acceptance",250,0.,1000.); + fTrackMultEMCAL->Sumw2(); + fOutput->Add(fTrackMultEMCAL); + + fClustMult = new TH1D ("hClustMult","Clusters multiplicity Distribution",250,0.,1000.); + fClustMult->Sumw2(); + fOutput->Add(fClustMult); + + fRecoPV = new TH1D("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5); + fRecoPV->Sumw2(); + fOutput->Add(fRecoPV); + + fPVZBefore = new TH1D ("hPVZDistr","Z Distribution for the Reconstructed Vertex",200,0.,40.); + fPVZBefore->Sumw2(); + fOutput->Add(fPVZBefore); + + fEtaPhiCell = new TH2D ("hEtaPhiCellActivity","",250,0.,1000., 250, 0., 1000.); + fEtaPhiCell->Sumw2(); + fOutput->Add(fEtaPhiCell); + + fEtaPhiClus = new TH2D ("hEtaPhiClusActivity","",250,0.,1000., 250, 0., 1000.); + fEtaPhiClus->Sumw2(); + fOutput->Add(fEtaPhiClus); + + fClusEvsClusT = new TH2D ("hEnVSTime","",250,0.,1000., 250,0.,1000.); + fClusEvsClusT->Sumw2(); + fOutput->Add(fClusEvsClusT); + + fDeltaETAClusTrackVSpT = new TH2D("hTC_Dz","Track-Cluster Dz vs pT of Cluster",100,0.,100.,100,-0.05,0.05); + fDeltaETAClusTrackVSpT->Sumw2(); + fOutput->Add(fDeltaETAClusTrackVSpT); + + fDeltaPHIClusTrackVSpT = new TH2D("hTC_Dx","Track-Cluster Dx vs pT of Cluster",100,0.,100.,100,-0.05,0.05); + fDeltaPHIClusTrackVSpT->Sumw2(); + fOutput->Add(fDeltaPHIClusTrackVSpT); + + } + //Initialize only the Common THistos for the Three different output + + fGoodEventsOnPVZ = new TH1D ("hGOODwrtPVZ","Number of Selected Events wrt Cut on Primary Vertex Z (0=disregarded,1=selected)",2,0.,2.); + fGoodEventsOnPVZ->Sumw2(); + fOutput->Add(fGoodEventsOnPVZ); + + fPT = new TH1D("hPt_NC","P_{T} distribution for Neutral Clusters",100,0.,100.); + fPT->Sumw2(); + fOutput->Add(fPT); + + fM02 = new TH2D("hM02_NC","M02 distribution for Neutral Clusters vs E",100,0.,100.,500,0.,5.); + fM02->Sumw2(); + fOutput->Add(fM02); + + fNLM = new TH1D("hNLM_NC","NLM distribution for Neutral Clusters",200,0.,4.); + fNLM->Sumw2(); + fOutput->Add(fNLM); + + fEtIsoCells = new TH1D("hEtIsoCell_NC","E_{T}^{iso cone} in iso cone distribution for Neutral Clusters with EMCAL Cells",100,0.,100.); + fEtIsoCells->SetXTitle("#Sigma E_{T}^{iso cone} (GeV/c)"); + fEtIsoCells->Sumw2(); + fOutput->Add(fEtIsoCells); + + fEtIsoClust = new TH1D("hEtIsoClus_NC","#Sigma E_{T}^{iso cone} in iso cone distribution for Neutral Clusters with EMCAL Clusters",100,0.,100.); + fEtIsoClust->SetXTitle("#Sigma E_{T}^{iso cone} (GeV/c)"); + fEtIsoClust->Sumw2(); + fOutput->Add(fEtIsoClust); + + fPtIsoTrack = new TH1D("hPtIsoTrack_NC"," #Sigma P_{T}^{iso cone} in iso cone distribution for Neutral Clusters with Tracks",100,0.,100.); + fPtIsoTrack->SetXTitle("#Sigma P_{T}^{iso cone} (GeV/c)"); + fPtIsoTrack->Sumw2(); + fOutput->Add(fPtIsoTrack); + + fPtEtIsoTC = new TH1D("hPtEtIsoTrackClust_NC","#Sigma P_{T}^{iso cone} + #Sigma E_{T}^{iso cone} in iso cone distribution for Neutral Clusters with Tracks and Clusters",100,0.,100.); + fPtEtIsoTC->SetXTitle("#Sigma P_{T}^{iso cone} + #Sigma E_{T}^{iso cone} (GeV/c)"); + fPtEtIsoTC->Sumw2(); + fOutput->Add(fPtEtIsoTC); + + fPhiBandUEClust = new TH2D(Form("hPhiBandUE_Cluster"),Form("UE Estimation with Phi Band Clusters"),100,0.,100.,100,0.,100.); + fPhiBandUEClust->SetXTitle("E_{T}"); + fPhiBandUEClust->SetYTitle("#Sigma E_{T}^{UE}"); + fPhiBandUEClust->Sumw2(); + fOutput->Add(fPhiBandUEClust); + + fEtaBandUEClust = new TH2D(Form("hEtaBandUE_Cluster"),Form("UE Estimation with Phi Band Clusters"),100,0.,100.,100,0.,100.); + fEtaBandUEClust->SetXTitle("E_{T}"); + fEtaBandUEClust->SetYTitle("#Sigma E_{T}^{UE}"); + fEtaBandUEClust->Sumw2(); + fOutput->Add(fEtaBandUEClust); + + fPhiBandUECells = new TH2D(Form("hPhiBandUE_CELLS"),Form("UE Estimation with Phi Band CELLS"),100,0.,100.,100,0.,100.); + fPhiBandUECells->SetXTitle("E_{T}"); + fPhiBandUECells->SetYTitle("#Sigma E_{T}^{UE}"); + fPhiBandUECells->Sumw2(); + fOutput->Add(fPhiBandUECells); + + fEtaBandUECells = new TH2D(Form("hEtaBandUE_CELLS"),Form("UE Estimation with Phi Band and CELLS"),100,0.,100.,100,0.,100.); + fEtaBandUECells->SetXTitle("E_{T}"); + fEtaBandUECells->SetYTitle("#Sigma E_{T}^{UE}"); + fEtaBandUECells->Sumw2(); + fOutput->Add(fEtaBandUECells); + + fPhiBandUETracks = new TH2D(Form("hPhiBandUE_TPC"),Form("UE Estimation with Phi Band TPC "),100,0.,100.,100,0.,100.); + fPhiBandUETracks->SetXTitle("E_{T}"); + fPhiBandUETracks->SetYTitle("#Sigma P_{T}^{UE}"); + fPhiBandUETracks->Sumw2(); + fOutput->Add(fPhiBandUETracks); + + fEtaBandUETracks = new TH2D(Form("hEtaBandUE_TPC"),Form("UE Estimation with Phi Band and TPC"),100,0.,100.,100,0.,100.); + fEtaBandUETracks->SetXTitle("E_{T}"); + fEtaBandUETracks->SetYTitle("#Sigma P_{T}^{UE}"); + fEtaBandUETracks->Sumw2(); + fOutput->Add(fEtaBandUETracks); + + fPerpConesUETracks = new TH2D("hConesUE","UE Estimation with Perpendicular Cones in TPC",100,0.,100.,100,0.,100.); + fPerpConesUETracks->SetXTitle("E_{T}"); + fPerpConesUETracks->SetYTitle("#Sigma P_{T}^{UE}"); + fPerpConesUETracks->Sumw2(); + fOutput->Add(fPerpConesUETracks); + + fTPCWithoutIsoConeB2BbandUE = new TH2D("hFullTPCUE","UE Estimation with almost Full TPC",100,0.,100.,100,0.,100.); + fPhiBandUEClust->SetXTitle("E_{T}"); + fPhiBandUEClust->SetYTitle("#Sigma E_{T}^{UE}"); + fTPCWithoutIsoConeB2BbandUE->Sumw2(); + fOutput->Add(fTPCWithoutIsoConeB2BbandUE); + + fEtIsolatedClust = new TH1D("hEtIsolatedClust","E_{T} distribution for Isolated Photons with clusters; #Sigma E_{T}^{iso cone}SetXTitle("E_{T}^{iso}"); + fEtIsolatedClust->Sumw2(); + fOutput->Add(fEtIsolatedClust); + + fEtIsolatedCells = new TH1D("hEtIsolatedCells","E_{T} distribution for Isolated Photons with cells; #Sigma E_{T}^{iso cone}SetXTitle("E_{T}^{iso}"); + fEtIsolatedCells->Sumw2(); + fOutput->Add(fEtIsolatedCells); + + fEtIsolatedTracks = new TH1D("hEtIsolatedTracks","E_{T} distribution for Isolated Photons with tracks; #Sigma P_{T}^{iso cone}SetXTitle("E_{T}^{iso}"); + fEtIsolatedTracks->Sumw2(); + fOutput->Add(fEtIsolatedTracks); + + fTest = new TH1D ("hTest","test cluster collection",100,-2.,6.); + fTest->Sumw2(); + fOutput->Add(fTest); + + + + + PostData(1, fOutput); + // return; +} + +// //________________________________________________________________________ +Float_t* AliAnalysisTaskEMCALPhotonIsolation::GenerateFixedBinArray(Int_t n, Float_t min, Float_t max) const +{ + Float_t *bins = new Float_t[n+1]; + + Float_t binWidth = (max-min)/n; + bins[0] = min; + for (Int_t i = 1; i <= n; i++) { + bins[i] = bins[i-1]+binWidth; + } + + return bins; +} + +//________________________________________________________________________ +void AliAnalysisTaskEMCALPhotonIsolation::ExecOnce() +{ + // Init the analysis. + + + + if (fParticleCollArray.GetEntriesFast()<2) { + AliError(Form("Wrong number of particle collections (%d), required 2",fParticleCollArray.GetEntriesFast())); + return; + } + + +// for (Int_t i = 0; i < 2; i++) { +// AliParticleContainer *contain = static_cast(fParticleCollArray.At(i)); +// // contain->GetClassName("AliEmcalParticle"); +// } + + AliAnalysisTaskEmcal::ExecOnce(); + if (!fInitialized) { + + AliError(Form("AliAnalysisTask not initialized")); + return; + } + + fevents+=1; + +} + +//______________________________________________________________________________________ +Bool_t AliAnalysisTaskEMCALPhotonIsolation::Run() +{ + + + + AliParticleContainer *clusters = static_cast(fParticleCollArray.At(1)); + AliEmcalParticle *emccluster = 0; + + // delete output USEFUL LATER FOR CONTAINER CREATION !! + //fOutClusters->Delete(); + + //Int_t clusCount = 0; AliError(Form("Should be here each time")); + // loop over all clusters + clusters->ResetCurrentID(); + Int_t index=-1; + + //Double_t ETleadingclust = 0., M02leadingcluster = 0., lambda0cluster = 0., phileadingclust = 0., etaleadingclust = 0., ptmc = 0.,mcptsum = 0.; + //Int_t Ntracks; + //Definition of the Array for Davide's Output + //const Int_t ndims = fNDimensions; + //Double_t outputValues[ndims]; + + + + if (fisLCAnalysis) { + // AliError(Form("Check the loop")); + // get the leading particle + + + + emccluster = static_cast(clusters->GetLeadingParticle()); + + if(!emccluster){ + + AliError(Form("no leading one")); + return kTRUE; + } + + + + // emccluster = clusters->GetLeadingParticle(); + index = emccluster->IdInCollection(); + + //add a command to get the index of the leading cluster! + if (!emccluster->IsEMCAL()) return kTRUE; + + AliVCluster *COI = emccluster->GetCluster(); + if (!COI) return kTRUE; + + TLorentzVector VecCOI; + COI->GetMomentum(VecCOI,fVertex); + + + if(!CheckBoundaries(VecCOI)) + return kTRUE; + + if(ClustTrackMatching(COI)) + return kTRUE; + + else + FillGeneralHistograms(COI,VecCOI, index); + } + else{ + //get the entries of the Cluster Container + //whatever is a RETURN in LCAnalysis here is a CONTINUE, + //since there are more than 1 Cluster per Event + + while((emccluster = static_cast(clusters->GetNextAcceptParticle()))){ + index++; + if (!emccluster->IsEMCAL()) return kTRUE; + + + AliVCluster *COI = emccluster->GetCluster(); + if(!COI) return kTRUE; + + TLorentzVector VecCOI; + COI->GetMomentum(VecCOI,fVertex); + + if(!CheckBoundaries(VecCOI)) + return kTRUE; + + + if(ClustTrackMatching(COI)) + continue; + + else + FillGeneralHistograms(COI,VecCOI, index); + + } + } + // PostData(1, fOutput); + return kTRUE; +} + + +// //__________________________________________________________________________ +Bool_t AliAnalysisTaskEMCALPhotonIsolation::ClustTrackMatching(AliVCluster *cluster) { + + Double_t deta=999; + Double_t dphi=999; + + AliParticleContainer *Tracks = static_cast(fParticleCollArray.At(0)); + // for an incoming cluster reject it if there is a track corresponding with a deta and dphi lower than the cuts + TLorentzVector nPart; + cluster->GetMomentum(nPart, fVertex); + + AliVTrack *mt = NULL; + AliAODCaloCluster *acl = dynamic_cast(cluster); + if(acl) { + if(acl->GetNTracksMatched()>1) + mt = static_cast(acl->GetTrackMatched(0)); + } + else { + AliESDCaloCluster *ecl = dynamic_cast(cluster); + Int_t im = ecl->GetTrackMatchedIndex(); + if(Tracks && im>=0) { + mt = static_cast(Tracks->GetParticle(im)); + } + } + // if(mt && mt->TestFilterBit(768)) { //Hybrid Tracks with AOD + AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta); + fDeltaETAClusTrackVSpT->Fill(nPart.Pt(), deta); + fDeltaPHIClusTrackVSpT->Fill(nPart.Pt(), dphi); + + + + + if(mt) return kTRUE; + + else return kFALSE; + +} + + + +//__________________________________________________________________________ +void AliAnalysisTaskEMCALPhotonIsolation::EtIsoCellPhiBand(TLorentzVector c, Float_t &EtIso, Float_t &PhiBandcells){ + + + AliEMCALGeometry* EmcalGeom = AliEMCALGeometry::GetInstance(); + + Float_t sumEnergyPhiBandCells=0., sumEnergyConeCells=0.; + + + // check the cell corresponding to the leading cluster + Int_t absId = 999; + //maybe best to call it LeadingCellIdinClus or better maxId since it won't be used anywhere else ??? + Bool_t CellLeadingClustId = EmcalGeom->GetAbsCellIdFromEtaPhi(c.Eta(),c.Phi(),absId); + if(!CellLeadingClustId) return; + + else{ + Int_t iTower = -1; + Int_t iModule = -1; + Int_t imEta=-1, imPhi=-1; + Int_t ieta =-1, iphi =-1; + + EmcalGeom->GetCellIndex(absId,iModule,iTower,imPhi,imEta); // to get the module, the tower, eta and phi for the cell corresponding to the leading cluster + EmcalGeom->GetCellPhiEtaIndexInSModule(iModule,iTower,imPhi,imEta,iphi,ieta); // to get the cell eta and phi in the super module for the cell coming from the leading cluster + + // Get the row and the column of the cell corresponding to the leading cluster in EMCAL + Int_t colcellleadingclust = ieta; + if(iModule % 2) colcellleadingclust = AliEMCALGeoParams::fgkEMCALCols + ieta ; // if the SM number is even you need to translate to have the one corresponding in EMCAL + Int_t rowcellleadingclust = iphi + AliEMCALGeoParams::fgkEMCALRows*int(iModule/2); // to have the corresponding row in EMCAL + + // total number or rows and columns in EMCAL + Int_t ntotalrows = AliEMCALGeoParams::fgkEMCALRows*16/3 ; // 5 + 2/3 supermodules in a row + Int_t ntotalcols = 2*AliEMCALGeoParams::fgkEMCALCols; // 2 supermodules in a column + + Int_t nbconesize = int(fIsoConeRadius/0.0143); //0.0143 tower acceptance + + // Get the cells + AliVCaloCells * cells =InputEvent()->GetEMCALCells(); + + // define the max and min row and column corresponding to the isolation cone around the seed cell from the leading cluster + Int_t irowmincone = rowcellleadingclust-nbconesize; + if(irowmincone<0) irowmincone=0; + + Int_t irowmaxcone = rowcellleadingclust+nbconesize; + if(irowmaxcone>AliEMCALGeoParams::fgkEMCALRows) irowmaxcone=AliEMCALGeoParams::fgkEMCALRows; // AliEMCALGeoParams::fgkEMCALRows = 24 in a supermodule + + Int_t icolmincone = colcellleadingclust - nbconesize; + if(icolmincone<0) icolmincone = 0; + + Int_t icolmaxcone = colcellleadingclust+nbconesize; + if(icolmaxcone>AliEMCALGeoParams::fgkEMCALCols) icolmaxcone=AliEMCALGeoParams::fgkEMCALCols; // AliEMCALGeoParams::fgkEMCALCols = 48 in a supermodule + + // loop on all cells + for(Int_t icol=0; icol AliEMCALGeoParams::fgkEMCALCols - 1){ // if the SM number is even the column isn't the one corresponding in the supermodule + inModule = 2*iSector; + icolLoc = icol-AliEMCALGeoParams::fgkEMCALCols; + } + Int_t irowLoc = irow - AliEMCALGeoParams::fgkEMCALRows*iSector ; + + if(TMath::Abs(icol-colcellleadingclust)nbconesize){ + if(irowirowmincone){ + Int_t iabsId = EmcalGeom->GetAbsCellIdFromCellIndexes(iModule,irow,icol); // verifier les irow et icol + sumEnergyPhiBandCells+=cells->GetAmplitude(iabsId); //should be Et + } + } + else if (TMath::Abs(icol-colcellleadingclust)>nbconesize && TMath::Abs(icol+colcellleadingclust)GetAbsCellIdFromCellIndexes(iModule,irowLoc,icolLoc); // verifier les irow et icol + sumEnergyConeCells+=cells->GetAmplitude(iabsId); //should be Et + } + } + } + } + EtIso = sumEnergyConeCells; + PhiBandcells = sumEnergyPhiBandCells; +} + + +//__________________________________________________________________________ +void AliAnalysisTaskEMCALPhotonIsolation::EtIsoCellEtaBand(TLorentzVector c, Float_t &EtIso, Float_t &EtaBandcells){ + + + + AliEMCALGeometry* EmcalGeom = AliEMCALGeometry::GetInstance(); + + Float_t sumEnergyEtaBandCells=0., sumEnergyConeCells=0.; + + + + // check the cell corresponding to the leading cluster + Int_t absId = 999; + //maybe best to call it LeadingCellIdinClus or better maxId since it won't be used anywhere else ??? + Bool_t CellLeadingClustId = EmcalGeom->GetAbsCellIdFromEtaPhi(c.Eta(),c.Phi(),absId); + if(!CellLeadingClustId) return; + + else{ + + Int_t iTower = -1; + Int_t iModule = -1; + Int_t imEta=-1, imPhi=-1; + Int_t ieta =-1, iphi =-1; + + EmcalGeom->GetCellIndex(absId,iModule,iTower,imPhi,imEta); // to get the module, the tower, eta and phi for the cell corresponding to the leading cluster + EmcalGeom->GetCellPhiEtaIndexInSModule(iModule,iTower,imPhi,imEta,iphi,ieta); // to get the cell eta and phi in the super module for the cell coming from the leading cluster + + // Get the row and the column of the cell corresponding to the leading cluster in EMCAL + Int_t colcellleadingclust = ieta; + if(iModule % 2) colcellleadingclust = AliEMCALGeoParams::fgkEMCALCols + ieta ; // if the SM number is even you need to translate to have the one corresponding in EMCAL + Int_t rowcellleadingclust = iphi + AliEMCALGeoParams::fgkEMCALRows*int(iModule/2); // to have the corresponding row in EMCAL + + // total number or rows and columns in EMCAL + Int_t ntotalrows = AliEMCALGeoParams::fgkEMCALRows*16/3 ; // 5 + 2/3 supermodules in a row + Int_t ntotalcols = 2*AliEMCALGeoParams::fgkEMCALCols; // 2 supermodules in a column + + Int_t nbconesize = int(fIsoConeRadius/0.0143); //0.0143 tower acceptance + + // Get the cells + AliVCaloCells * cells =InputEvent()->GetEMCALCells(); + + // define the max and min row and column corresponding to the isolation cone around the seed cell from the leading cluster + Int_t irowmincone = rowcellleadingclust-nbconesize; + if(irowmincone<0) irowmincone=0; + + Int_t irowmaxcone = rowcellleadingclust+nbconesize; + if(irowmaxcone>AliEMCALGeoParams::fgkEMCALRows) irowmaxcone=AliEMCALGeoParams::fgkEMCALRows; // AliEMCALGeoParams::fgkEMCALRows = 24 in a supermodule + + Int_t icolmincone = colcellleadingclust-nbconesize; + if(icolmincone<0) icolmincone = 0; + + Int_t icolmaxcone = colcellleadingclust+nbconesize; + if(icolmaxcone>AliEMCALGeoParams::fgkEMCALCols) icolmaxcone=AliEMCALGeoParams::fgkEMCALCols; // AliEMCALGeoParams::fgkEMCALCols = 48 in a supermodule + + // loop on all cells + for(Int_t icol=0; icol AliEMCALGeoParams::fgkEMCALCols - 1){ // if the SM number is even the column isn't the one corresponding in the supermodule + inModule = 2*iSector; + icolLoc = icol-AliEMCALGeoParams::fgkEMCALCols; + } + + Int_t irowLoc = irow - AliEMCALGeoParams::fgkEMCALRows*iSector ; + + if(TMath::Abs(icol-colcellleadingclust)nbconesize){ + if(icolicolmincone){ + Int_t iabsId = EmcalGeom->GetAbsCellIdFromCellIndexes(iModule,irow,icol); // verifier les irow et icol + sumEnergyEtaBandCells+=cells->GetAmplitude(iabsId); //should be Et + } + } + else if (TMath::Abs(icol-colcellleadingclust)>nbconesize && TMath::Abs(icol+colcellleadingclust)GetAbsCellIdFromCellIndexes(iModule,irowLoc,icolLoc); // verifier les irow et icol + sumEnergyConeCells+=cells->GetAmplitude(iabsId); //should be Et + } + } + } + } + EtIso = sumEnergyConeCells; + EtaBandcells = sumEnergyEtaBandCells; +} + + +//__________________________________________________________________________ +void AliAnalysisTaskEMCALPhotonIsolation::EtIsoClusPhiBand(TLorentzVector c, Float_t &EtIso, Float_t &PhiBandclus, Int_t index){ + + Float_t sumEnergyPhiBandClus=0., sumEnergyConeClus=0.; + + //needs a check on the same cluster + AliParticleContainer *clusters = static_cast(fParticleCollArray.At(1)); + AliEmcalParticle *clust; + + Int_t localindex=-1; + while((clust = static_cast(clusters->GetNextAcceptParticle()))){ //check the position of other clusters in respect to the leading cluster + localindex++; + if(localindex==index) continue; + + AliVCluster *cluster= clust->GetCluster(); + + TLorentzVector nClust; //STILL NOT INITIALIZED + cluster->GetMomentum(nClust,fVertex); + Float_t phiclust =nClust.Phi(); + Float_t etaclust= nClust.Eta(); + //redefine phi/c.Eta() from the cluster we passed to the function + + Float_t radius = TMath::Sqrt(TMath::Power(phiclust- c.Phi(),2)+TMath::Power(etaclust-c.Eta(),2)); // define the radius between the leading cluster and the considered cluster + + if(radius>fIsoConeRadius){ // the cluster is outside the isolation cone in this case study UE + + // actually phi band here + if(TMath::Abs(etaclust - c.Eta()) < fIsoConeRadius){ + sumEnergyPhiBandClus += nClust.Pt(); + } + } + else // if the cluster is in the isolation cone, add the cluster pT + sumEnergyConeClus += nClust.Pt(); + + } + EtIso = sumEnergyConeClus; + PhiBandclus = sumEnergyPhiBandClus; +} + + +//__________________________________________________________________________ +void AliAnalysisTaskEMCALPhotonIsolation::EtIsoClusEtaBand(TLorentzVector c, Float_t &EtIso, Float_t &EtaBandclus, Int_t index){ + + Float_t sumEnergyEtaBandClus =0., sumEnergyConeClus=0.; + AliParticleContainer *clusters = static_cast(fParticleCollArray.At(1)); + + AliEmcalParticle *clust; + + Int_t localindex=-1; + while((clust = static_cast(clusters->GetNextAcceptParticle()))){ //check the position of other clusters in respect to the leading cluster + localindex++; + if(localindex==index) continue; + + AliVCluster *cluster= clust->GetCluster(); + + TLorentzVector nClust; //STILL NOT INITIALIZED + cluster->GetMomentum(nClust,fVertex); + + Float_t phiclust =nClust.Phi(); + Float_t etaclust= nClust.Eta(); + //redefine phi/c.Eta() from the cluster we passed to the function + + // define the radius between the leading cluster and the considered cluster + Float_t radius = TMath::Sqrt(TMath::Power(phiclust-c.Phi(),2)+TMath::Power(etaclust-c.Eta(),2)); + + if(radius>fIsoConeRadius){ // the cluster is outside the isolation cone in this case study UE + + // actually eta band here + if(TMath::Abs(etaclust - c.Eta()) < fIsoConeRadius){ + sumEnergyEtaBandClus += nClust.Pt(); + } + } + else // if the cluster is in the isolation cone, add the cluster pT + sumEnergyConeClus += nClust.Pt(); + + } + EtIso = sumEnergyConeClus; + EtaBandclus = sumEnergyEtaBandClus; +} + + +//__________________________________________________________________________ +void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackPhiBand(TLorentzVector c, Float_t &PtIso, Float_t &PhiBandtrack){ + + //INSERT BOUNDARIES ACCORDING TO THE FLAG WE WANT TO USE + Float_t sumpTConeCharged=0., sumpTPhiBandTrack=0.; + Float_t minphi= 0., maxphi= 2*TMath::Pi(), mineta = -0.9, maxeta= 0.9; + + AliParticleContainer *Tracks = static_cast(fParticleCollArray.At(0)); + + if(!fTPC4Iso){ + mineta = -0.7; + maxeta = 0.7; + minphi = 1.4; + maxphi = TMath::Pi(); + } + + AliVParticle *track = Tracks->GetNextAcceptParticle(0); + + while(track){ + + //CHECK IF TRACK IS IN BOUNDARIES + Float_t phitrack = track->Phi(); + Float_t etatrack = track->Eta(); + // define the radius between the leading cluster and the considered cluster + //redefine phi/c.Eta() from the cluster we passed to the function + if(phitrack < maxphi && phitrack > minphi && etatrack < maxeta && etatrack > mineta){ + Float_t radius = TMath::Sqrt(TMath::Power(phitrack - c.Phi(),2)+TMath::Power(etatrack - c.Eta(),2)); + + if(radius>fIsoConeRadius){ // if tracks are outside the isolation cone study + + // actually phi band here --- ADD Boundaries conditions + if(TMath::Abs(etatrack - c.Eta()) < fIsoConeRadius){ + sumpTPhiBandTrack += track->Pt(); + } + } + else + sumpTConeCharged+=track->Pt(); // should not double count if the track matching is already done + } + track=Tracks->GetNextAcceptParticle(); + } + PtIso = sumpTConeCharged; + PhiBandtrack = sumpTPhiBandTrack; +} + + +//__________________________________________________________________________ +void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackEtaBand(TLorentzVector c, Float_t &PtIso, Float_t &EtaBandtrack){ + + //INSERT BOUNDARIES ACCORDING TO THE FLAG WE WANT TO USE + Float_t sumpTConeCharged=0., sumpTEtaBandTrack=0.; + Float_t minphi= 0., maxphi= 2*TMath::Pi(), mineta = -0.9, maxeta= 0.9; + + AliParticleContainer *Tracks = static_cast(fParticleCollArray.At(0)); + + + if(!fTPC4Iso){ + mineta = -0.7; + maxeta = 0.7; + minphi = 1.4; + maxphi = TMath::Pi(); + } + + + AliVParticle *track = Tracks->GetNextAcceptParticle(0); + + while(track){ + + Float_t phitrack = track->Phi(); + Float_t etatrack = track->Eta(); + //redefine phi/c.Eta() from the cluster we passed to the function + if(phitrack < maxphi && phitrack > minphi && etatrack < maxeta && etatrack > mineta){ + Float_t radius = TMath::Sqrt(TMath::Power(phitrack-c.Phi(),2)+TMath::Power(etatrack-c.Eta(),2)); // define the radius between the leading cluster and the considered cluster + if(radius>fIsoConeRadius){ // if tracks are outside the isolation cone study UE + + // actually eta band here --- ADD Boundaries conditions + if(TMath::Abs(phitrack - c.Phi()) < fIsoConeRadius){ + sumpTEtaBandTrack += track->Pt(); + } + } + else + sumpTConeCharged += track->Pt(); // should not double count if the track matching is already done + } + track=Tracks->GetNextAcceptParticle(); + } + PtIso = sumpTConeCharged; + EtaBandtrack = sumpTEtaBandTrack; +} + + +//__________________________________________________________________________ +void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackOrthCones(TLorentzVector c, Float_t &PtIso, Float_t &Cones){ + + Float_t sumpTConeCharged=0., sumpTPerpConeTrack=0.; + + AliParticleContainer *Tracks = static_cast(fParticleCollArray.At(0)); + + Float_t EtaClus = c.Eta(); + Float_t PhiClus = c.Phi(); + Float_t PhiCone1 = PhiClus - TMath::PiOver2(); + Float_t PhiCone2 = PhiClus + TMath::PiOver2(); + + if (PhiCone1 < 0.) PhiCone1 += 2*TMath::Pi(); + + AliVParticle *track = Tracks->GetNextAcceptParticle(0); + + while(track){ + + Float_t phitrack = track->Phi(); + Float_t etatrack = track->Eta(); + + Float_t Dist2Clust = TMath::Sqrt(TMath::Power(etatrack-EtaClus, 2)+TMath::Power(phitrack-PhiClus, 2)); + if (Dist2ClustPt(); + + else{//tracks outside the IsoCone + //Distances from the centres of the two Orthogonal Cones + Float_t Dist2Cone1 = TMath::Sqrt(TMath::Power(etatrack-EtaClus, 2)+TMath::Power(phitrack-PhiCone1, 2)); + Float_t Dist2Cone2 = TMath::Sqrt(TMath::Power(etatrack-EtaClus, 2)+TMath::Power(phitrack-PhiCone2, 2)); + + //Is the Track Inside one of the two Cones ->Add to UE + if((Dist2Cone1 < fIsoConeRadius) || (Dist2Cone2 < fIsoConeRadius)) + sumpTPerpConeTrack += track->Pt(); + } + track=Tracks->GetNextAcceptParticle(); + } + PtIso = sumpTConeCharged; + Cones = sumpTPerpConeTrack; +} + +//__________________________________________________________________________ +void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackFullTPC(TLorentzVector c, Float_t &PtIso, Float_t &Full){ + + Float_t sumpTConeCharged=0., sumpTTPCexceptB2B=0.; + + AliParticleContainer *Tracks = static_cast(fParticleCollArray.At(0)); + AliVParticle *track = Tracks->GetNextAcceptParticle(0); + + while(track){ + + Float_t phitrack = track->Phi(); + Float_t etatrack = track->Eta(); + //redefine phi/c.Eta() from the cluster we passed to the function + + Float_t radius = TMath::Sqrt(TMath::Power(phitrack-c.Phi(),2)+TMath::Power(etatrack-c.Eta(),2)); // define the radius between the leading cluster and the considered cluster + if(radius>fIsoConeRadius){ // if tracks are outside the isolation cone study UE + Double_t dphiup = c.Phi() + TMath::Pi() - fIsoConeRadius; + Double_t dphidown = c.Phi() + TMath::Pi() + fIsoConeRadius; + // TPC except B2B + if(TMath::Power(phitrack-c.Phi(),2) +TMath::Power(etatrack-c.Eta(),2) > TMath::Power((fIsoConeRadius+0.1),2) && phitrack < dphidown && phitrack> dphiup) + sumpTTPCexceptB2B += track->Pt(); + } + else + sumpTConeCharged += track->Pt(); // should not double count if the track matching is already done + + track=Tracks->GetNextAcceptParticle(); + } + PtIso = sumpTConeCharged; + Full = sumpTTPCexceptB2B; +} + +//__________________________________________________________________________ +Bool_t AliAnalysisTaskEMCALPhotonIsolation::CheckBoundaries(TLorentzVector VecCOI){ + + Float_t minphibound= 1.5 , minetabound= -0.5, maxphibound= TMath::Pi()-0.1, maxetabound= 0.5; + Bool_t isINBoundaries; + + if(!fTPC4Iso){ + minphibound = 1.4+0.4; //to be changed with fIsoConeR + maxphibound = TMath::Pi()-0.4; + minetabound = -0.7+0.4; + maxetabound = 0.7-0.4; + } + if(VecCOI.Eta() > maxetabound || VecCOI.Eta() < minetabound || VecCOI.Phi() > maxphibound || VecCOI.Phi() >minphibound) + isINBoundaries=kFALSE; + else + isINBoundaries=kTRUE; + + return isINBoundaries; +} + +//__________________________________________________________________________ +Bool_t AliAnalysisTaskEMCALPhotonIsolation::FillGeneralHistograms(AliVCluster *COI, TLorentzVector VecCOI, Int_t index){ + + + AliParticleContainer *Tracks = static_cast(fParticleCollArray.At(0)); + + AliEmcalParticle *emctrack = 0; + + int NTracks=0; + Tracks->ResetCurrentID(); + while ((emctrack = static_cast(Tracks->GetNextAcceptParticle()))) { + AliVTrack *track = emctrack->GetTrack(); + if(!track) continue; + // if(!(track->TestFilterBit("kHybrid"))) continue; + + NTracks++; + } + + Double_t ETCOI = 0., M02COI = 0., lambda0cluster = 0., phiCOI = 0., etaCOI = 0., ptmc = 0., mcptsum = 0.; + //Int_t Ntracks; + //Definition of the Array for Davide's Output + const Int_t ndims = fNDimensions; + Double_t outputValues[ndims]; + + ETCOI = VecCOI.Et(); + M02COI = COI->GetM02(); + etaCOI = VecCOI.Eta(); + phiCOI = VecCOI.Phi(); + + // ******** Isolation and UE calculation with different methods ********* + + Double_t EtThreshold = 5; + + switch(fEtIsoMethod) + { + case 0: // SumEtFill(ETCOI); + } + break; + } + + //DO NOT CHANGE EVER AGAIN THE FOLLOWING DEFINITIONS + Float_t IsoConeArea = TMath::Pi()*fIsoConeRadius*fIsoConeRadius; + Float_t etabandArea = (1.4 - 2.*fIsoConeRadius)*2.*fIsoConeRadius; + Float_t phibandArea = ((5./9.)*TMath::Pi()- 2.*fIsoConeRadius)*2.*fIsoConeRadius; + Float_t etabandAreaTr = (1.8 - 2.*fIsoConeRadius)*2.*fIsoConeRadius; + Float_t phibandAreaTr = (2.*TMath::Pi() - 4.*fIsoConeRadius)*2.*fIsoConeRadius;//there is a 2 more because of the JET CONE B2B + Float_t PerpConesArea = 2.*IsoConeArea; + Float_t FullTPCArea = 1.8*2.*TMath::Pi()-fIsoConeRadius*(TMath::Pi()*fIsoConeRadius + 3.6); + + Float_t ISOLATION, UE; + + if(!fTPC4Iso){ //EMCAL Only for Acceptance of Clusters + switch(fIsoMethod) + { + case 0: //EMCAL CELLS + + switch (fUEMethod) + { + case 0: //phi band + EtIsoCellPhiBand(VecCOI, ISOLATION, UE); + //Normalisation UE wrt Area - TO DO- + UE = UE * (IsoConeArea / phibandArea); + fPhiBandUECells->Fill(VecCOI.Pt() , UE); + fEtIsoCells->Fill(ISOLATION); + if(ISOLATIONFill(ETCOI); + fEtisoT=ETCOI; + fPtisoT=VecCOI.Pt(); + } + break; + case 1: //eta band + EtIsoCellEtaBand(VecCOI, ISOLATION, UE); + //Normalisation UE wrt Area - TO DO- + UE = UE * (IsoConeArea / etabandArea); + fEtaBandUECells->Fill(VecCOI.Pt() , UE); + fEtIsoCells->Fill(ISOLATION); + if(ISOLATIONFill(ETCOI); + fEtisoT=ETCOI; + fPtisoT=VecCOI.Pt(); + } + break; + } + break; + + case 1: //EMCAL CLUSTERS + + switch (fUEMethod) + { + case 0: //phi band + EtIsoClusPhiBand(VecCOI, ISOLATION, UE,index); + //Normalisation UE wrt Area - TO DO- + UE = UE * (IsoConeArea / phibandArea); + fPhiBandUEClust->Fill(VecCOI.Pt() , UE); + fEtIsoClust->Fill(ISOLATION); + break; + case 1: //eta band + EtIsoClusEtaBand(VecCOI, ISOLATION, UE,index); + //Normalisation UE wrt Area - TO DO- + UE = UE * (IsoConeArea / etabandArea); + fEtaBandUEClust->Fill(VecCOI.Pt() , UE); + fEtIsoClust->Fill(ISOLATION); + if(ISOLATIONFill(ETCOI); + fEtisoT=ETCOI; + fPtisoT=VecCOI.Pt(); + } + break; + } + case 2: //TRACKS (TBD which tracks) //EMCAL Tracks + switch (fUEMethod) + { + case 0: //phi band + PtIsoTrackPhiBand(VecCOI, ISOLATION, UE); + //Normalisation UE wrt Area - TO DO- + UE = UE * (IsoConeArea / phibandAreaTr); + fPhiBandUETracks->Fill(VecCOI.Pt() , UE); + fPtIsoTrack->Fill(ISOLATION); + if(ISOLATIONFill(ETCOI); + fEtisoT=ETCOI; + fPtisoT=VecCOI.Pt(); + } + break; + case 1: //eta band + PtIsoTrackEtaBand(VecCOI, ISOLATION, UE); + //Normalisation UE wrt Area - TO DO- + UE = UE * (IsoConeArea / etabandAreaTr); + fEtaBandUETracks->Fill(VecCOI.Pt() , UE); + fPtIsoTrack->Fill(ISOLATION); + if(ISOLATIONFill(ETCOI); + fEtisoT=ETCOI; + fPtisoT=VecCOI.Pt(); + } + break; + // case 2: //Cones + // PtIsoTrackOrthCones(VecCOI, absId, ISOLATION, UE); + // break; + // case 3: //Full TPC + // PtIsoTrackFullTPC(VecCOI, absId, ISOLATION, UE); + // break; + } + } + } + else{ //EMCAL + TPC (Only Tracks for the Isolation since IsoCone Goes Out of EMCAL) + switch (fUEMethod) + { + case 0: //phi band + PtIsoTrackPhiBand(VecCOI, ISOLATION, UE); + //Normalisation UE wrt Area - TO DO- + UE = UE * (IsoConeArea / phibandAreaTr); + fPhiBandUETracks->Fill(VecCOI.Pt() , UE); + fPtIsoTrack->Fill(ISOLATION); + if(ISOLATIONFill(ETCOI); + fEtisoT=ETCOI; + fPtisoT=VecCOI.Pt(); + } + break; + case 1: //eta band + PtIsoTrackEtaBand(VecCOI, ISOLATION, UE); + //Normalisation UE wrt Area - TO DO- + UE = UE * (IsoConeArea / etabandAreaTr); + fEtaBandUETracks->Fill(VecCOI.Pt() , UE); + fPtIsoTrack->Fill(ISOLATION); + if(ISOLATIONFill(ETCOI); + fEtisoT=ETCOI; + fPtisoT=VecCOI.Pt(); + } + break; + case 2: //Cones + PtIsoTrackOrthCones(VecCOI, ISOLATION, UE); + //Normalisation UE wrt Area - TO DO- + UE = UE * (IsoConeArea / PerpConesArea); + fPerpConesUETracks ->Fill(VecCOI.Pt() , UE); + fPtIsoTrack->Fill(ISOLATION); + if(ISOLATIONFill(ETCOI); + fEtisoT=ETCOI; + fPtisoT=VecCOI.Pt(); + } + break; + case 3: //Full TPC + PtIsoTrackFullTPC(VecCOI, ISOLATION, UE); + //Normalisation UE wrt Area - TO DO- + UE = UE * (IsoConeArea / FullTPCArea); + fTPCWithoutIsoConeB2BbandUE->Fill(VecCOI.Pt() , UE); + fPtIsoTrack->Fill(ISOLATION); + if(ISOLATIONFill(ETCOI); + fEtisoT=ETCOI; + fPtisoT=VecCOI.Pt(); + } + break; + } + } + + + //Here we should call something to know the number of tracks... + //Soon I'll put in this version the "old way"; please let me know if + //any of you could do the same with the JET framework + + switch(fWho) { + case 0: + flambda0T=M02COI; + fEtT=VecCOI.Et(); + fPtT=VecCOI.Pt(); + fetaT=VecCOI.Eta(); + fphiT=VecCOI.Phi(); + fsumEtisoconeT=ISOLATION; + // AliError(Form("lambda 0 %f",flambda0T)); + fsumEtUE=UE; + + fOutputTree->Fill(); + break; + + case 1: + outputValues[0] = NTracks; + outputValues[1] = ETCOI; + outputValues[2] = lambda0cluster; + outputValues[3] = ISOLATION; + outputValues[4] = UE; + outputValues[5] = ISOLATION-UE; + outputValues[6] = VecCOI.Eta(); + outputValues[7] = VecCOI.Phi(); + if (fIsMC) { + outputValues[8] = ptmc; + outputValues[9] = mcptsum; + } + fOutputTHnS -> Fill(outputValues); + break; + // // fOutPTMAX -> Fill(outputValues[1],outputValues[2],); + } + return kTRUE; +} + + + diff --git a/PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhotonIsolation.h b/PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhotonIsolation.h new file mode 100644 index 00000000000..f5b60fd3a9d --- /dev/null +++ b/PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhotonIsolation.h @@ -0,0 +1,187 @@ +#ifndef ALIANALYSISTASKEMCALPHOTONISOLATION_H +#define ALIANALYSISTASKEMCALPHOTONISOLATION_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + + ////////////////////////////////////////////////////////////////////////// + // // + // Task for Isolated Gamma in p-p,p-Pb and eventually g-h Correlation // + // // + // Author: Davide Francesco Lodato (Utrecht University) // + // Lucile Ronflette (Subatech, Nantes) // + // Marco Marquard (University Frankfurt am Main) // + ////////////////////////////////////////////////////////////////////////// + + //ROOT System + +class TH1; +class TH2; +class TH3; +class THnSparse; +class TList; +class TObjArray; +class AliEMCALGeometry; +class AliESDCaloCells; +class AliESDEvent; +class AliESDtrack; +class TClonesArray; +class TList; +class TString; +class AliVCluster; +class AliVParticle; +class AliESDtrackCuts; +class AliAODEvent; +class AliAODCaloCells; +class AliVCluster; +class AliMCEvent; +class AliStack; +class TParticle; +class AliClusterContainer; +class AliParticleContainer; +class AliEmcalParticle; + //AliRoot System +class AliEMCALTrack; + //class AliMagF; +class AliEMCALRecoUtils; + //class AliAnalysisFilter; +class AliAODTrack; +class AliAODCaloCluster; +class AliESDCaloCluster; +class AliVCaloCells; + //class AliEventPoolManager; + +#include "AliAnalysisTaskEmcal.h" + +class AliAnalysisTaskEMCALPhotonIsolation : public AliAnalysisTaskEmcal { +public: + AliAnalysisTaskEMCALPhotonIsolation(); + AliAnalysisTaskEMCALPhotonIsolation(const char *name, Bool_t histo=kFALSE); + virtual ~AliAnalysisTaskEMCALPhotonIsolation(); + + void UserCreateOutputObjects(); + + void SetIsoConeRadius(Float_t r) { fIsoConeRadius = r ;} + void SetCTMdeltaEta (Float_t r) { fdetacut = r ;} + void SetCTMdeltaPhi (Float_t r) { fdphicut = r ;} + void SetIsoMethod (Int_t r ) { fIsoMethod = r ;} + void SetUEMethod (Int_t UE) { fUEMethod = UE;} + void SetOutputFormat (Int_t iOut) { fWho = iOut;} + void SetQA (Bool_t QA) { fQA = QA;} + void SetMC (Bool_t MC) { fIsMC = MC;} + void SetUSEofTPC (Bool_t TPC) { fTPC4Iso = TPC;} + void SetLCAnalysis (Bool_t LC) { fisLCAnalysis = LC;} + +protected: + + void EtIsoCellPhiBand(TLorentzVector c, Float_t &EtIso, Float_t &PhiBand); //EIsoCone via Cells UE via PhiBand EMCAL + void EtIsoCellEtaBand(TLorentzVector c, Float_t &EtIso, Float_t &EtaBand); //EIsoCone via Cells UE via EtaBand EMCAL + void EtIsoClusPhiBand(TLorentzVector c, Float_t &EtIso, Float_t &EtaBand, Int_t index); //EIsoCone via Clusters UE via EtaBand EMCAL + void EtIsoClusEtaBand(TLorentzVector c, Float_t &EtIso, Float_t &EtaBand, Int_t index); //EIsoCone via Clusters UE via EtaBand EMCAL + void PtIsoTrackPhiBand(TLorentzVector c, Float_t &PtIso, Float_t &PhiBand); //PIsoCone via Track UE via PhiBand TPC + void PtIsoTrackEtaBand(TLorentzVector c, Float_t &PtIso, Float_t &EtaBand); //PIsoCone via Track UE via EtaBand TPC + // void PtIsoTraClusPhiBand(TLorentzVector c, Float_t &PtIso, Float_t &PhiBand); //(P+E)IsoCone via Track/Clus UE via PhiBand TPC+EMCAL + // void PtIsoTraClusEtaBand(TLorentzVector c, Float_t &PtIso, Float_t &EtaBand); //(P+E)IsoCone via Track/Clus UE via EtaBand TPC+EMCAL + void PtIsoTrackOrthCones(TLorentzVector c, Float_t &PtIso, Float_t &Cones); //PIsoCone via Tracks UE via Orthogonal Cones in Phi + void PtIsoTrackFullTPC(TLorentzVector c, Float_t &PtIso, Float_t &Full); //PIsoCone via Tracks UE via FullTPC - IsoCone - B2BEtaBand + Bool_t ClustTrackMatching(AliVCluster *cluster); + Bool_t CheckBoundaries(TLorentzVector VecCOI); + // void FillNCOutput(AliVCluster *COI, TLorentzVector VecCOI, Int_t index); + + Float_t* GenerateFixedBinArray(Int_t n, Float_t min, Float_t max) const; + void ExecOnce(); + Bool_t Run(); + + using AliAnalysisTaskEmcal::FillGeneralHistograms; + Bool_t FillGeneralHistograms(AliVCluster *COI, TLorentzVector VecCOI, Int_t index); + //Bool_t FillGeneralHistograms(AliVCluster *COI, TLorentzVector VecCOI, Int_t index); + + + // TObjArray fParticleCollArray; // Neutral Clusters collection array + TClonesArray *fNCluster; // Neutral clusters + + Int_t fWho; // MODE for the Output Object (TTree or THnSparse) + + + + + + //IMPLEMENT ALL THE HISTOGRAMS AND ALL THE OUTPUT OBJECTS WE WANT!!! + // TList *fOutputList; //! Output list + TGeoHMatrix *fGeomMatrix[12];//! Geometry misalignment matrices for EMCal + + + TH1 *fTrackMult; //!Track Multiplicity ---QA + TH1 *fTrackMultEMCAL; //!Track Multiplicity EMCAL ---QA + TH1 *fClustMult; //!Cluster Multiplicity EMCAL ---QA + TH1 *fPVZBefore; //!Z Vertex distribution before cuts. ---QA + TH2 *fEtaPhiCell; //!EMCAL Active Cells Distribution EtaPhi ---QA + TH2 *fEtaPhiClus; //!EMCAL Cluster Distribution EtaPhi ---QA + TH2 *fClusEvsClusT; //!Cluster Energy vs Cluster Time ---QA + TH1 *fGoodEventsOnPVZ; //!Number of selected events After the Cut on Primary Vertex + TH1 *fPT; //!Pt distribution + TH2 *fM02; //!Squared_Lambda0 distribution + TH1 *fNLM; //!NLM distribution + TH2 *fDeltaETAClusTrackVSpT; //!dEta Cluster-Track VS pT! + TH2 *fDeltaPHIClusTrackVSpT; //!dPhi Cluster-Track VS pT! + TH1 *fEtIsoCells; //!Isolation Energy with EMCAL Cells + TH1 *fEtIsoClust; //!Isolation Energy with EMCAL Clusters + TH1 *fPtIsoTrack; //!Isolation Pt with Tracks + TH1 *fPtEtIsoTC; //!Isolation with Pt from Tracks and Et from NON-Matched Clusters + TH2 *fPhiBandUEClust; //!UE with Phi Band (Clusters) + TH2 *fEtaBandUEClust; //!UE with Eta Band (Clusters) + TH2 *fPhiBandUECells; //!UE with Phi Band (Cells) + TH2 *fEtaBandUECells; //!UE with Eta Band (Cells) + TH2 *fPhiBandUETracks; //!UE with Phi Band (Tracks) + TH2 *fEtaBandUETracks; //!UE with Eta Band (Tracks) + TH2 *fPerpConesUETracks; //!UE with Cones (Tracks ONLY) + TH2 *fTPCWithoutIsoConeB2BbandUE; //!UE with Full TPC except IsoCone and EtaBand in Back2Back + TH1 *fNTotClus10GeV; //!number of TOTAL clusters with Energy bigger than 10 GeV + TH1 *fRecoPV; //! primary vertex reconstruction + TH1 *fEtIsolatedCells; //! Isolated photons, isolation with cells + TH1 *fEtIsolatedClust; //! Isolated photons, isolation with clusters + TH1 *fEtIsolatedTracks; //! Isolated photons, isolation with tracks + TH1 *fTest; //! + + THnSparse *fOutputTHnS; //! 1st Method 4 Output + THnSparse *fOutPTMAX; //! 1st Method 4 Isolation on pTMax + + TTree *fOutputTree; //! 2nd Method 4 Output + + Float_t fIsoConeRadius; // Radius for the Isolation Cont + Int_t fEtIsoMethod; // Isolation definition 0=SumEtSetOutputFormat(iOutput); + task->SetLCAnalysis(kFALSE); + task->SetQA(kTRUE); + task->SetIsoMethod(1); + task->SetUEMethod(1); + task->SetUSEofTPC(kFALSE); + + + AliParticleContainer *trackCont = task->AddParticleContainer(ntracks); + AliParticleContainer *clusterCont = task->AddParticleContainer(nclusters); + + printf("Task for neutral cluster analysis created and configured, pass it to AnalysisManager\n"); + // #### Add analysis task + manager->AddTask(task); + + // AliAnalysisDataContainer *contHistos = manager->CreateContainer(myContName.Data(), TList::Class(), AliAnalysisManager::kOutputContainer, Form("%s:NeutralCluster", AliAnalysisManager::GetCommonFileName())); + AliAnalysisDataContainer *contHistos = manager->CreateContainer(myContName.Data(), TList::Class(), AliAnalysisManager::kOutputContainer,Form("%s:NeutralCluster",AliAnalysisManager::GetCommonFileName())); + AliAnalysisDataContainer *cinput = manager->GetCommonInputContainer(); + manager->ConnectInput(task, 0, cinput); + manager->ConnectOutput(task, 1, contHistos); + + /*if(isEMCalTrain) + RequestMemory(task,200*1024); + + // #### Do some nasty piggybacking on demand + if (externalMacro) + gROOT->Macro(externalMacro); +*/ + return task; +} diff --git a/PWGGA/EMCALTasks/macros/AddTaskPhotonPreparation.C b/PWGGA/EMCALTasks/macros/AddTaskPhotonPreparation.C new file mode 100644 index 00000000000..a526852d9eb --- /dev/null +++ b/PWGGA/EMCALTasks/macros/AddTaskPhotonPreparation.C @@ -0,0 +1,99 @@ +AliAnalysisTaskSE* AddTaskPhotonPreparation( + const char* periodstr = "LHC11h", + const char* pTracksName = "PicoTracks", + const char* usedMCParticles = "MCParticlesSelected", + const char* usedClusters = "CaloClusters", + const UInt_t pSel = AliVEvent::kAny, + const Bool_t doHistos = kTRUE, + const Bool_t makePicoTracks = kTRUE, + const Bool_t makeTrigger = kTRUE, + const Bool_t isEmcalTrain = kFALSE, + const Double_t trackeff = 1.0, + const Bool_t doAODTrackProp = kTRUE, + const Bool_t modifyMatchObjs = kTRUE, + const Int_t iOutput = 1 +) +{ + + printf("Preparing neutral cluster analysis\n"); + + // #### Define manager and data container names + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error("AddTaskNeutralCluster", "No analysis manager found."); + return NULL; + } + + AliVEventHandler *evhand = mgr->GetInputEventHandler(); + if (!evhand) { + Error("AddTasNeutralCluster", "This task requires an input event handler"); + return NULL; + } + TString period(periodstr); + TString clusterColName(usedClusters); + TString particleColName(usedMCParticles); + TString picoTracksName(pTracksName); + + TString dType("ESD"); + if (!evhand->InheritsFrom("AliESDInputHandler")) + dType = "AOD"; + if ((dType == "AOD") && (clusterColName == "CaloClusters")) + clusterColName = "caloClusters"; + if ((dType == "ESD") && (clusterColName == "caloClusters")) + clusterColName = "CaloClusters"; + + if (0) { + gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEmcalTrackPropagator.C"); + + cout<<"AddTaskEmcalTrackPropagator"<SelectCollisionCandidates(pSel); + } + + //----------------------- Trigger Maker ----------------------------------------------------- + if (makeTrigger) { + gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEmcalTriggerMaker.C"); + AliEmcalTriggerMaker *emcalTriggers = AddTaskEmcalTriggerMaker("EmcalTriggers"); + + cout<<"AddTaskEmcalTriggerMaker"<SelectCollisionCandidates(pSel); + } + + //----------------------- Track Matching tasks ----------------------------------------------------- + gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskMatchingChain.C"); + AliEmcalClusTrackMatcherTask *emcalClus = AddTaskMatchingChain(periodstr,pSel, + clusterColName, + trackeff,doAODTrackProp, + 0.1,modifyMatchObjs,doHistos); + + cout<<"AddTaskMatchingChain"<LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEmcalPicoTrackMaker.C"); + AliEmcalPicoTrackMaker *pTrackTask = AddTaskEmcalPicoTrackMaker(picoTracksName, inputTracks); + pTrackTask->SelectCollisionCandidates(pSel); + } + + + printf("Creating container names for cluster analysis\n"); + TString myContName(""); + + myContName = Form("Photon_Preperation"); + + + gROOT->LoadMacro("$ALICE_ROOT/PWGGA/EMCALTasks/macros/AddTaskEMCALPhotonIsolation.C"); + + AliAnalysisTaskEMCALPhotonIsolation *task =AddTaskEMCALPhotonIsolation(emctracks,emcclusters,kTRUE, iOutput,kFALSE); + task->SelectCollisionCandidates(pSel); + if(isEmcalTrain) + RequestMemory(task,500*1024); + + + return task; +} diff --git a/PWGGA/PWGGAEMCALTasksLinkDef.h b/PWGGA/PWGGAEMCALTasksLinkDef.h index 515f5168190..faff47f9c9b 100644 --- a/PWGGA/PWGGAEMCALTasksLinkDef.h +++ b/PWGGA/PWGGAEMCALTasksLinkDef.h @@ -30,6 +30,7 @@ #pragma link C++ class AliAnalysisTaskEMCALPi0V2ShSh+; #pragma link C++ class AliEMCalpi0ClusterEvaluationTask+; #pragma link C++ class AliAnalysisTaskEMCALIsolation+; +#pragma link C++ class AliAnalysisTaskEMCALPhotonIsolation+; #pragma link C++ class AliAnalysisTaskEMCALMesonGGSDM+; #pragma link C++ class AliAnalysisTaskEMCALMesonGGSDMpPb+; #pragma link C++ class AliAnalysisTaskSDMGammaMC+; -- 2.43.0