including task to produce final histograms on isolated photon analysis
authormcosenti <mcosenti@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Jun 2012 14:02:40 +0000 (14:02 +0000)
committermcosenti <mcosenti@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Jun 2012 14:02:40 +0000 (14:02 +0000)
PWGGA/CMakelibPWGGAEMCALTasks.pkg
PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx [new file with mode: 0644]
PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.h [new file with mode: 0644]
PWGGA/EMCALTasks/macros/AddTaskEMCALIsoPhoton.C [new file with mode: 0644]
PWGGA/PWGGAEMCALTasksLinkDef.h

index 001188b..eeb920e 100644 (file)
@@ -52,6 +52,7 @@ set ( SRCS
  EMCALTasks/AliAnalysisTaskEmcal.cxx
  EMCALTasks/AliEmcalClusTrackMatcherTask.cxx
  EMCALTasks/AliAnalysisTaskPi0V2.cxx
+ EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx
 )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
diff --git a/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx b/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx
new file mode 100644 (file)
index 0000000..518b006
--- /dev/null
@@ -0,0 +1,550 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliESDEvent.h"
+#include "AliESDHeader.h"
+#include "AliESDUtils.h"
+#include "AliESDInputHandler.h"
+#include "AliESDpid.h"
+#include "AliKFParticle.h"
+
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+
+#include "AliESDtrack.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDv0.h"
+#include "AliV0vertexer.h"
+#include "AliESDCaloCluster.h"
+#include "AliESDCaloCells.h"
+#include "AliEMCALGeometry.h"
+#include "AliEMCALRecoUtils.h"
+#include "TLorentzVector.h"
+#include "AliVCluster.h"
+
+
+#include "AliAnalysisTaskEMCALIsoPhoton.h"
+#include "TFile.h"
+
+
+ClassImp(AliAnalysisTaskEMCALIsoPhoton)
+
+//________________________________________________________________________
+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),                   //!primary vertex Z before cut
+    fCellAbsIdVsAmpl(0),         //!cell abs id vs cell amplitude (energy)
+    fNClusHighClusE(0),          //!total number of clusters vs. highest clus energy in the event
+    fM02Et(0),                   //!M02 vs Et for all clusters
+    fM02EtTM(0),                 //!M02 vs Et for clusters with track-match (dEta=0.01 && dPhi=0.025)
+    fM02EtCeIso1(0),             //!M02 vs Et for clusters with isolation neutral Et<1GeV
+    fM02EtCeIso2(0),             //!M02 vs Et for clusters with isolation neutral Et<2GeV
+    fM02EtCeIso5(0),             //!M02 vs Et for clusters with isolation neutral Et<5GeV
+    fM02EtTrIso1(0),             //!M02 vs Et for clusters with isolation charged Et<1GeV
+    fM02EtTrIso2(0),             //!M02 vs Et for clusters with isolation charged Et<2GeV
+    fM02EtTrIso5(0),             //!M02 vs Et for clusters with isolation charged Et<5GeV
+    fM02EtAllIso1(0),            //!M02 vs Et for clusters with isolation total Et<1GeV
+    fM02EtAllIso2(0),            //!M02 vs Et for clusters with isolation total Et<2GeV
+    fM02EtAllIso5(0),            //!M02 vs Et for clusters with isolation total Et<5GeV
+    fCeIsoVsEtPho(0),            //!Neutral isolation Et vs. cluster Et, 0.05<M02<0.30
+    fTrIsoVsEtPho(0),            //!Charged isolation Et vs. cluster Et, 0.05<M02<0.30
+    fAllIsoVsEtPho(0),           //!Total isolation Et vs. cluster Et, 0.05<M02<0.30
+    //track matched stuff
+    fM02EtCeIso1TM(0),           //!Track-matched M02 vs Et for clusters with isolation neutral Et<1GeV
+    fM02EtCeIso2TM(0),           //!Track-matched M02 vs Et for clusters with isolation neutral Et<2GeV
+    fM02EtCeIso5TM(0),           //!Track-matched M02 vs Et for clusters with isolation neutral Et<5GeV
+    fM02EtTrIso1TM(0),           //!Track-matched M02 vs Et for clusters with isolation charged Et<1GeV
+    fM02EtTrIso2TM(0),           //!Track-matched M02 vs Et for clusters with isolation charged Et<2GeV
+    fM02EtTrIso5TM(0),           //!Track-matched M02 vs Et for clusters with isolation charged Et<5GeV
+    fM02EtAllIso1TM(0),          //!Track-matched M02 vs Et for clusters with isolation total Et<1GeV
+    fM02EtAllIso2TM(0),          //!Track-matched M02 vs Et for clusters with isolation total Et<2GeV
+    fM02EtAllIso5TM(0),          //!Track-matched M02 vs Et for clusters with isolation total Et<5GeV
+    fCeIsoVsEtPhoTM(0),          //!Track-matched Neutral isolation Et vs. cluster Et, 0.05<M02<0.30
+    fTrIsoVsEtPhoTM(0),          //!Track-matched Charged isolation Et vs. cluster Et, 0.05<M02<0.30
+    fAllIsoVsEtPhoTM(0)         //!Track-matched Total isolation Et vs. cluster Et, 0.05<M02<0.30
+
+
+  
+{
+  // 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 *) 
+{
+  //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; 
+  // Main loop
+  // Called for each event
+
+  // Post output data.
+  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());
+    }
+  } //track loop 
+  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()
+{
+  if(!fCaloClusters)
+    return;
+  const Int_t nclus = fCaloClusters->GetEntries();
+  if(nclus==0)
+    return;
+  Double_t maxE;
+  Int_t nthresholds = 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);
+    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;
+    Double_t Et = c->E()*TMath::Sin(clsVec.Theta());
+    fM02Et->Fill(Et, c->GetM02());
+    if(ceiso>1)
+      fM02EtCeIso1->Fill(Et, c->GetM02());
+    if(ceiso>2)
+      fM02EtCeIso2->Fill(Et, c->GetM02());
+    if(ceiso>5)
+      fM02EtCeIso5->Fill(Et, c->GetM02());
+    if(triso>1)
+      fM02EtTrIso1->Fill(Et, c->GetM02());
+    if(triso>2)
+      fM02EtTrIso2->Fill(Et, c->GetM02());
+    if(triso>5)
+      fM02EtTrIso5->Fill(Et, c->GetM02());
+    if(alliso>1)
+      fM02EtAllIso1->Fill(Et, c->GetM02());
+    if(alliso>2)
+      fM02EtAllIso2->Fill(Et, c->GetM02());
+    if(alliso>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>1)
+       fM02EtCeIso1TM->Fill(Et, c->GetM02());
+      if(ceiso>2)
+       fM02EtCeIso2TM->Fill(Et, c->GetM02());
+      if(ceiso>5)
+       fM02EtCeIso5TM->Fill(Et, c->GetM02());
+      if(triso>1)
+       fM02EtTrIso1TM->Fill(Et, c->GetM02());
+      if(triso>2)
+       fM02EtTrIso2TM->Fill(Et, c->GetM02());
+      if(triso>5)
+       fM02EtTrIso5TM->Fill(Et, c->GetM02());
+      if(alliso>1)
+       fM02EtAllIso1TM->Fill(Et, c->GetM02());
+      if(alliso>2)
+       fM02EtAllIso2TM->Fill(Et, c->GetM02());
+      if(alliso>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)
+{
+  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)
+{
+  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 *) 
+{
+  // Draw result to the screen
+  // Called once at the end of the query
+
+}
diff --git a/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.h b/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.h
new file mode 100644 (file)
index 0000000..2461a95
--- /dev/null
@@ -0,0 +1,152 @@
+#ifndef AliAnalysisTaskEMCALIsoPhoton_cxx\r
+#define AliAnalysisTaskEMCALIsoPhoton_cxx\r
+\r
+class TH1;\r
+class TH2;\r
+class TObjArray;\r
+class AliESDEvent;\r
+class AliESDtrack;\r
+class AliESDtrackCuts;\r
+class AliVCluster;\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class AliAnalysisTaskEMCALIsoPhoton : public AliAnalysisTaskSE {\r
+ public:\r
+  AliAnalysisTaskEMCALIsoPhoton() : \r
+  AliAnalysisTaskSE(), \r
+  \r
+    fCaloClusters(0),\r
+    fSelPrimTracks(0),\r
+    fEMCalCells(0),\r
+    fPrTrCuts(0),\r
+\r
+    fGeom(0x0),\r
+  \r
+    fESD(0),\r
+  \r
+    fOutputList(0),\r
+\r
+    fEvtSel(0),\r
+\r
+    fPVtxZ(0),                   //!primary vertex Z before cut\r
+    fCellAbsIdVsAmpl(0),         //!cell abs id vs cell amplitude (energy)\r
+    fNClusHighClusE(0),          //!total number of clusters vs. highest clus energy in the event\r
+    fM02Et(0),                   //!M02 vs Et for all clusters\r
+    fM02EtTM(0),                 //!M02 vs Et for clusters with track-match (dEta=0.01 && dPhi=0.025)\r
+    fM02EtCeIso1(0),             //!M02 vs Et for clusters with isolation neutral Et<1GeV\r
+    fM02EtCeIso2(0),             //!M02 vs Et for clusters with isolation neutral Et<2GeV\r
+    fM02EtCeIso5(0),             //!M02 vs Et for clusters with isolation neutral Et<5GeV\r
+    fM02EtTrIso1(0),             //!M02 vs Et for clusters with isolation charged Et<1GeV\r
+    fM02EtTrIso2(0),             //!M02 vs Et for clusters with isolation charged Et<2GeV\r
+    fM02EtTrIso5(0),             //!M02 vs Et for clusters with isolation charged Et<5GeV\r
+    fM02EtAllIso1(0),            //!M02 vs Et for clusters with isolation total Et<1GeV\r
+    fM02EtAllIso2(0),            //!M02 vs Et for clusters with isolation total Et<2GeV\r
+    fM02EtAllIso5(0),            //!M02 vs Et for clusters with isolation total Et<5GeV\r
+    fCeIsoVsEtPho(0),            //!Neutral isolation Et vs. cluster Et, 0.05<M02<0.30\r
+    fTrIsoVsEtPho(0),            //!Charged isolation Et vs. cluster Et, 0.05<M02<0.30\r
+    fAllIsoVsEtPho(0),           //!Total isolation Et vs. cluster Et, 0.05<M02<0.30\r
+    //track matched stuff\r
+    fM02EtCeIso1TM(0),           //!Track-matched M02 vs Et for clusters with isolation neutral Et<1GeV\r
+    fM02EtCeIso2TM(0),           //!Track-matched M02 vs Et for clusters with isolation neutral Et<2GeV\r
+    fM02EtCeIso5TM(0),           //!Track-matched M02 vs Et for clusters with isolation neutral Et<5GeV\r
+    fM02EtTrIso1TM(0),           //!Track-matched M02 vs Et for clusters with isolation charged Et<1GeV\r
+    fM02EtTrIso2TM(0),           //!Track-matched M02 vs Et for clusters with isolation charged Et<2GeV\r
+    fM02EtTrIso5TM(0),           //!Track-matched M02 vs Et for clusters with isolation charged Et<5GeV\r
+    fM02EtAllIso1TM(0),          //!Track-matched M02 vs Et for clusters with isolation total Et<1GeV\r
+    fM02EtAllIso2TM(0),          //!Track-matched M02 vs Et for clusters with isolation total Et<2GeV\r
+    fM02EtAllIso5TM(0),          //!Track-matched M02 vs Et for clusters with isolation total Et<5GeV\r
+    fCeIsoVsEtPhoTM(0),          //!Track-matched Neutral isolation Et vs. cluster Et, 0.05<M02<0.30\r
+    fTrIsoVsEtPhoTM(0),          //!Track-matched Charged isolation Et vs. cluster Et, 0.05<M02<0.30\r
+    fAllIsoVsEtPhoTM(0)          //!Track-matched Total isolation Et vs. cluster Et, 0.05<M02<0.30\r
+\r
+\r
+    \r
+  \r
+  \r
+  {}\r
+  AliAnalysisTaskEMCALIsoPhoton(const char *name);\r
+  virtual ~AliAnalysisTaskEMCALIsoPhoton() {}\r
+\r
+  void   UserCreateOutputObjects();\r
+  void   UserExec(Option_t *option);\r
+  void   Terminate(Option_t *);\r
+\r
+  void         SetGeoName(const char *n)              { fGeoName            = n;       }\r
+  void         SetPeriod(const char *n)               { fPeriod             = n;       }\r
+  void         SetTrainMode(Bool_t t)                 { fIsTrain            = t;       }\r
+  void         SetExotCut(Double_t c)                 { fExoticCut          = c;       }\r
+  void         SetIsoConeR(Double_t r)                { fIsoConeR           = r;       }\r
+  void         SetPrimTrackCuts(AliESDtrackCuts *c)   { fPrTrCuts           = c;       }\r
+  void         FillClusHists();\r
+  void         GetCeIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core);\r
+  void         GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core);\r
+  Double_t     GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax);\r
+  Double_t     GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const; \r
+\r
+\r
+\r
+\r
+  \r
+ protected:\r
+  TRefArray             *fCaloClusters;          //!pointer to EMCal clusters\r
+  TObjArray             *fSelPrimTracks;         //!pointer to ESD primary tracks\r
+  AliESDCaloCells       *fEMCalCells;            //!pointer to EMCal cells\r
+  AliESDtrackCuts       *fPrTrCuts;              //!pointer to hold the prim track cuts\r
+  AliEMCALGeometry      *fGeom;                   // geometry utils\r
+  TString               fGeoName;                // geometry name (def = EMCAL_FIRSTYEARV1)\r
+  TString               fPeriod;                 // string to the LHC period\r
+  Bool_t                fIsTrain;                //variable to set train mode\r
+  Double_t              fExoticCut;              //variable to set the cut on exotic clusters\r
+  Double_t              fIsoConeR;               //variable to set the isolation cone radius\r
+  \r
+ private:\r
+  AliESDEvent *fESD;      //! ESD object\r
+\r
+\r
+  \r
+  TList       *fOutputList; //! Output list\r
+  //histograms for events with 1+ track pt>1\r
+  \r
+  TH1F        *fEvtSel;                  //!evt selection counter: 0=all trg, 1=pv cut \r
+  TH1F        *fPVtxZ;                   //!primary vertex Z before cut\r
+  TH2F        *fCellAbsIdVsAmpl;         //!cell abs id vs cell amplitude (energy)\r
+  TH2F        *fNClusHighClusE;          //!total number of clusters vs. highest clus energy in the event\r
+  TH2F        *fM02Et;                   //!M02 vs Et for all clusters\r
+  TH2F        *fM02EtTM;                 //!M02 vs Et for clusters with track-match (dEta=0.01 && dPhi=0.025)\r
+  TH2F        *fM02EtCeIso1;             //!M02 vs Et for clusters with isolation neutral Et<1GeV\r
+  TH2F        *fM02EtCeIso2;             //!M02 vs Et for clusters with isolation neutral Et<2GeV\r
+  TH2F        *fM02EtCeIso5;             //!M02 vs Et for clusters with isolation neutral Et<5GeV\r
+  TH2F        *fM02EtTrIso1;             //!M02 vs Et for clusters with isolation charged Et<1GeV\r
+  TH2F        *fM02EtTrIso2;             //!M02 vs Et for clusters with isolation charged Et<2GeV\r
+  TH2F        *fM02EtTrIso5;             //!M02 vs Et for clusters with isolation charged Et<5GeV\r
+  TH2F        *fM02EtAllIso1;            //!M02 vs Et for clusters with isolation total Et<1GeV\r
+  TH2F        *fM02EtAllIso2;            //!M02 vs Et for clusters with isolation total Et<2GeV\r
+  TH2F        *fM02EtAllIso5;            //!M02 vs Et for clusters with isolation total Et<5GeV\r
+  TH2F        *fCeIsoVsEtPho;            //!Neutral isolation Et vs. cluster Et, 0.05<M02<0.30\r
+  TH2F        *fTrIsoVsEtPho;            //!Charged isolation Et vs. cluster Et, 0.05<M02<0.30\r
+  TH2F        *fAllIsoVsEtPho;           //!Total isolation Et vs. cluster Et, 0.05<M02<0.30\r
+  //track matched stuff\r
+  TH2F        *fM02EtCeIso1TM;           //!Track-matched M02 vs Et for clusters with isolation neutral Et<1GeV\r
+  TH2F        *fM02EtCeIso2TM;           //!Track-matched M02 vs Et for clusters with isolation neutral Et<2GeV\r
+  TH2F        *fM02EtCeIso5TM;           //!Track-matched M02 vs Et for clusters with isolation neutral Et<5GeV\r
+  TH2F        *fM02EtTrIso1TM;           //!Track-matched M02 vs Et for clusters with isolation charged Et<1GeV\r
+  TH2F        *fM02EtTrIso2TM;           //!Track-matched M02 vs Et for clusters with isolation charged Et<2GeV\r
+  TH2F        *fM02EtTrIso5TM;           //!Track-matched M02 vs Et for clusters with isolation charged Et<5GeV\r
+  TH2F        *fM02EtAllIso1TM;          //!Track-matched M02 vs Et for clusters with isolation total Et<1GeV\r
+  TH2F        *fM02EtAllIso2TM;          //!Track-matched M02 vs Et for clusters with isolation total Et<2GeV\r
+  TH2F        *fM02EtAllIso5TM;          //!Track-matched M02 vs Et for clusters with isolation total Et<5GeV\r
+  TH2F        *fCeIsoVsEtPhoTM;          //!Track-matched Neutral isolation Et vs. cluster Et, 0.05<M02<0.30\r
+  TH2F        *fTrIsoVsEtPhoTM;          //!Track-matched Charged isolation Et vs. cluster Et, 0.05<M02<0.30\r
+  TH2F        *fAllIsoVsEtPhoTM;         //!Track-matched Total isolation Et vs. cluster Et, 0.05<M02<0.30\r
+\r
+\r
+\r
+   \r
+  AliAnalysisTaskEMCALIsoPhoton(const AliAnalysisTaskEMCALIsoPhoton&); // not implemented\r
+  AliAnalysisTaskEMCALIsoPhoton& operator=(const AliAnalysisTaskEMCALIsoPhoton&); // not implemented\r
+  \r
+  ClassDef(AliAnalysisTaskEMCALIsoPhoton, 1); // example of analysis\r
+};\r
+\r
+#endif\r
diff --git a/PWGGA/EMCALTasks/macros/AddTaskEMCALIsoPhoton.C b/PWGGA/EMCALTasks/macros/AddTaskEMCALIsoPhoton.C
new file mode 100644 (file)
index 0000000..834c247
--- /dev/null
@@ -0,0 +1,58 @@
+
+
+AliAnalysisTaskEMCALIsoPhoton *AddTaskEMCALIsoPhoton()
+{
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTaskEMCALIsoPhoton", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // Create the task and configure it.
+  //===========================================================================
+  AliAnalysisTaskEMCALIsoPhoton* ana = new  AliAnalysisTaskEMCALIsoPhoton("");
+  
+  ana->SelectCollisionCandidates( AliVEvent::kEMC1 | AliVEvent::kMB | AliVEvent::kEMC7 | AliVEvent::kINT7);
+  
+  Bool_t isMC = (mgr->GetMCtruthEventHandler() != NULL);
+
+  //ana->SetClusThreshold(clusTh);
+  ana->SetTrainMode(kTRUE);
+  //ana->SetGridMode(kTRUE);
+  // ana->SetMcMode(isMC);
+  
+  AliESDtrackCuts *cutsp = new AliESDtrackCuts;
+  cutsp->SetMinNClustersTPC(70);
+  cutsp->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
+  cutsp->SetMaxChi2PerClusterTPC(4);
+  cutsp->SetRequireTPCRefit(kTRUE);
+  cutsp->SetAcceptKinkDaughters(kFALSE);
+  cutsp->SetMaxDCAToVertexZ(3.2);
+  cutsp->SetMaxDCAToVertexXY(2.4);
+  cutsp->SetDCAToVertex2D(kTRUE);
+  cutsp->SetPtRange(0.2);
+  cutsp->SetEtaRange(-1.0,1.0);
+  ana->SetPrimTrackCuts(cutsp);
+  ana->SetPeriod(period.Data());
+  if(period.Contains("11"))
+    ana->SetGeoName("EMCAL_COMPLETEV1");
+  else
+    ana->SetGeoName("EMCAL_FIRSTYEARV1");
+
+  
+  mgr->AddTask(ana);
+  
+  // Create ONLY the output containers for the data produced by the task.
+  // Get and connect other common input/output containers via the manager as below
+  //==============================================================================
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("histosEMCALIsoPhoton", 
+                                                           TList::Class(),AliAnalysisManager::kOutputContainer,
+                                                           Form("%s", AliAnalysisManager::GetCommonFileName()));
+  
+  mgr->ConnectInput  (ana, 0, mgr->GetCommonInputContainer());
+  mgr->ConnectOutput (ana, 1, coutput1 );
+   
+  return ana;
+}
index 67b26e0..8f3cd93 100644 (file)
@@ -41,5 +41,6 @@
 #pragma link C++ class AliEmcalParticleMaker+;
 #pragma link C++ class AliEmcalClusTrackMatcherTask+;
 #pragma link C++ class AliAnalysisTaskPi0V2+;
+#pragma link C++ class AliAnalysisTaskEMCALIsoPhoton+;
 
 #endif