PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx

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