PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhotonIsolation.cxx

   1   // $Id$
   2   //
   3   // Emcal Neutral Cluster analysis base task.
   4   //
   5   // Authors: D.Lodato,L.Ronflette, M.Marquard
   6
   7
   8
   9 #include <TClonesArray.h>
  10 #include <TChain.h>
  11 #include <TList.h>
  12 #include <TVector3.h>
  13 #include <TLorentzVector.h>
  14 #include <TH1D.h>
  15 #include <TH2D.h>
  16 #include <TH3D.h>
  17 #include <THnSparse.h>
  18 #include "AliAnalysisManager.h"
  19 #include "AliCentrality.h"
  20 #include "AliEMCALGeometry.h"
  21 #include "AliESDEvent.h"
  22 #include "AliAODEvent.h"
  23 #include "AliLog.h"
  24 #include "AliVCluster.h"
  25 #include "AliVEventHandler.h"
  26 #include "AliVParticle.h"
  27 #include "AliClusterContainer.h"
  28 #include "AliVTrack.h"
  29 #include "AliEmcalParticle.h"
  30 #include "AliParticleContainer.h"
  31 #include "AliAODCaloCluster.h"
  32 #include "AliESDCaloCluster.h"
  33 #include "AliVCaloCells.h"
  34 #include "AliPicoTrack.h"
  35 #include "AliAODMCParticle.h"
  36 #include "AliAODMCHeader.h"
  37 #include "AliEMCALRecoUtils.h"
  38
  39
  40
  41 #include "AliAnalysisTaskEMCALPhotonIsolation.h"
  42
  43
  44 ClassImp(AliAnalysisTaskEMCALPhotonIsolation)
  45
  46   //________________________________________________________________________
  47 AliAnalysisTaskEMCALPhotonIsolation::AliAnalysisTaskEMCALPhotonIsolation() :
  48 AliAnalysisTaskEmcal("AliAnalysisTaskEMCALPhotonIsolation",kTRUE),
  49   //fParticleCollArray(),
  50 fAOD(0),
  51 fVevent(0),
  52 fNCluster(0),
  53 fAODMCParticles(0),
  54 fTracksAna(0),
  55 fStack(0),
  56 fWho(-1),
  57   //fOutputList(0),
  58 fTrackMult(0),
  59 fTrackMultEMCAL(0),
  60 fClustMult(0),
  61 fPVZBefore(0),
  62 fEtaPhiCell(0),
  63 fEtaPhiClus(0),
  64 fClusEvsClusT(0),
  65 fVz(0),
  66 fEvents(0),
  67 fPT(0),
  68 fE(0),
  69 fPtaftTime(0),
  70 fPtaftTM(0),
  71 fPtaftFC(0),
  72 fPtaftM02C(0),
  73 fClusTime(0),
  74 fM02(0),
  75 fNLM(0),
  76 fDeltaETAClusTrack(0),
  77 fDeltaPHIClusTrack(0),
  78 fDeltaETAClusTrackMatch(0),
  79 fDeltaPHIClusTrackMatch(0),
  80 fDeltaETAClusTrackVSpT(0),
  81 fDeltaPHIClusTrackVSpT(0),
  82 fEtIsoCells(0),
  83 fEtIsoClust(0),
  84 fPtIsoTrack(0),
  85 fPtEtIsoTC(0),
  86 fPhiBandUEClust(0),
  87 fEtaBandUEClust(0),
  88 fPhiBandUECells(0),
  89 fEtaBandUECells(0),
  90 fPhiBandUETracks(0),
  91 fEtaBandUETracks(0),
  92 fPerpConesUETracks(0),
  93 fTPCWithoutIsoConeB2BbandUE(0),
  94 fNTotClus10GeV(0),
  95 fRecoPV(0),
  96 fEtIsolatedCells(0),
  97 fEtIsolatedClust(0),
  98 fPtIsolatedNClust(0),
  99 fPtIsolatedNTracks(0),
 100 fEtIsolatedTracks(0),
 101 fPtvsM02iso(0),
 102 fPtvsM02noiso(0),
 103 fTestIndex(0),
 104 fTestIndexE(0),
 105 fTestLocalIndexE(0),
 106 fTestEnergyCone(0),
 107 fTestEtaPhiCone(0),
 108 fOutputTHnS(0),
 109 fOutMCTruth(0),
 110 fOutClustMC(0),
 111 fOutputQATree(0),
 112 fOutputTree(0),
 113 fphietaPhotons(0),
 114 fphietaOthers(0),
 115 fphietaOthersBis(0),
 116 fIsoConeRadius(0.4),
 117 fEtIsoMethod(0),
 118 fEtIsoThreshold(2),
 119 fdetacut(0.025),
 120 fdphicut(0.03),
 121 fM02mincut(0.1),
 122 fM02maxcut(0.3),
 123 fQA(0),
 124 fIsMC(0),
 125 fTPC4Iso(0),
 126 fIsoMethod(0),
 127 fUEMethod(0),
 128 fNDimensions(0),
 129 fMCDimensions(0),
 130 fMCQAdim(0),
 131 fisLCAnalysis(0),
 132 fTest1(0),
 133 fTest2(0),
 134 fEClustersT(0),
 135 fPtClustersT(0),
 136 fEtClustersT(0),
 137 fEtaClustersT(0),
 138 fPhiClustersT(0),
 139 fM02ClustersT(0),
 140 fevents(0),
 141 fNClustersT(0),
 142 flambda0T(0),
 143 fM02isoT(0),
 144 fM02noisoT(0),
 145 fPtnoisoT(0),
 146 fEtT(0),
 147 fPtT(0),
 148 fPtisoT(0),
 149 fEtisolatedT(0),
 150 fPtisolatedT(0),
 151 fetaT(0),
 152 fphiT(0),
 153 fsumEtisoconeT(0),
 154 fsumEtUE(0)
 155   //tracks(0),
 156   //clusters(0)
 157
 158 {
 159     // Default constructor.
 160
 161     //fParticleCollArray.SetOwner(kTRUE);
 162     // for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
 163
 164   SetMakeGeneralHistograms(kTRUE);
 165 }
 166
 167   //________________________________________________________________________
 168 AliAnalysisTaskEMCALPhotonIsolation::AliAnalysisTaskEMCALPhotonIsolation(const char *name, Bool_t histo) :
 169 AliAnalysisTaskEmcal(name, histo),
 170   //fParticleCollArray(),
 171 fAOD(0),
 172 fVevent(0),
 173 fNCluster(0),
 174 fAODMCParticles(0),
 175 fTracksAna(0),
 176 fStack(0),
 177 fWho(-1),
 178   //fOutputList(0),
 179 fTrackMult(0),
 180 fTrackMultEMCAL(0),
 181 fClustMult(0),
 182 fPVZBefore(0),
 183 fEtaPhiCell(0),
 184 fEtaPhiClus(0),
 185 fClusEvsClusT(0),
 186 fVz(0),
 187 fEvents(0),
 188 fPT(0),
 189 fE(0),
 190 fPtaftTime(0),
 191 fPtaftTM(0),
 192 fPtaftFC(0),
 193 fPtaftM02C(0),
 194 fClusTime(0),
 195 fM02(0),
 196 fNLM(0),
 197 fDeltaETAClusTrack(0),
 198 fDeltaPHIClusTrack(0),
 199 fDeltaETAClusTrackMatch(0),
 200 fDeltaPHIClusTrackMatch(0),
 201 fDeltaETAClusTrackVSpT(0),
 202 fDeltaPHIClusTrackVSpT(0),
 203 fEtIsoCells(0),
 204 fEtIsoClust(0),
 205 fPtIsoTrack(0),
 206 fPtEtIsoTC(0),
 207 fPhiBandUEClust(0),
 208 fEtaBandUEClust(0),
 209 fPhiBandUECells(0),
 210 fEtaBandUECells(0),
 211 fPhiBandUETracks(0),
 212 fEtaBandUETracks(0),
 213 fPerpConesUETracks(0),
 214 fTPCWithoutIsoConeB2BbandUE(0),
 215 fNTotClus10GeV(0),
 216 fRecoPV(0),
 217 fEtIsolatedCells(0),
 218 fEtIsolatedClust(0),
 219 fPtIsolatedNClust(0),
 220 fPtIsolatedNTracks(0),
 221 fEtIsolatedTracks(0),
 222 fPtvsM02iso(0),
 223 fPtvsM02noiso(0),
 224 fTestIndex(0),
 225 fTestIndexE(0),
 226 fTestLocalIndexE(0),
 227 fTestEnergyCone(0),
 228 fTestEtaPhiCone(0),
 229 fOutputTHnS(0),
 230 fOutMCTruth(0),
 231 fOutClustMC(0),
 232 fOutputQATree(0),
 233 fOutputTree(0),
 234 fphietaPhotons(0),
 235 fphietaOthers(0),
 236 fphietaOthersBis(0),
 237 fIsoConeRadius(0.4),
 238 fEtIsoMethod(0),
 239 fEtIsoThreshold(2),
 240 fdetacut(0.025),
 241 fdphicut(0.03),
 242 fM02mincut(0.1),
 243 fM02maxcut(0.3),
 244 fQA(0),
 245 fIsMC(0),
 246 fTPC4Iso(0),
 247 fIsoMethod(0),
 248 fUEMethod(0),
 249 fNDimensions(0),
 250 fMCDimensions(0),
 251 fMCQAdim(0),
 252 fisLCAnalysis(0),
 253 fTest1(0),
 254 fTest2(0),
 255 fEClustersT(0),
 256 fPtClustersT(0),
 257 fEtClustersT(0),
 258 fEtaClustersT(0),
 259 fPhiClustersT(0),
 260 fM02ClustersT(0),
 261 fevents(0),
 262 fNClustersT(0),
 263 flambda0T(0),
 264 fM02isoT(0),
 265 fM02noisoT(0),
 266 fPtnoisoT(0),
 267 fEtT(0),
 268 fPtT(0),
 269 fPtisoT(0),
 270 fEtisolatedT(0),
 271 fPtisolatedT(0),
 272 fetaT(0),
 273 fphiT(0),
 274 fsumEtisoconeT(0),
 275 fsumEtUE(0)
 276   //tracks(0),
 277   //clusters(0)
 278
 279 {
 280     // Standard constructor.
 281
 282     //fParticleCollArray.SetOwner(kTRUE);
 283     //    for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
 284
 285   SetMakeGeneralHistograms(kTRUE);
 286 }
 287
 288   //________________________________________________________________________
 289 AliAnalysisTaskEMCALPhotonIsolation::~AliAnalysisTaskEMCALPhotonIsolation(){
 290     // Destructor
 291 }
 292
 293
 294   //________________________________________________________________________
 295 void AliAnalysisTaskEMCALPhotonIsolation::UserCreateOutputObjects(){
 296     // Create ouput histograms and THnSparse and TTree
 297
 298   AliAnalysisTaskEmcal::UserCreateOutputObjects();
 299
 300
 301   if ((fIsoMethod == 0 || fIsoMethod == 1) && fTPC4Iso){
 302     cout<<"Error: Iso_Methods with CELLS and CLUSTERS work only within EMCAL "<<endl;
 303     cout<<"Please Set Iso_Method and TPC4Iso Accordingly!!"<<endl;
 304     return;
 305   }
 306   if ((fIsoMethod == 0 || fIsoMethod == 1) && fUEMethod> 1) {
 307     cout<<"Error: UE_Methods with CELLS and CLUSTERS work only within EMCAL"<<endl;
 308     cout<<"Please Set Iso_Method and UE_Method Accordingly!!"<<endl;
 309     return;
 310   }
 311
 312   if(fUEMethod>1 && !fTPC4Iso){
 313     cout<<"Please set UE Method Accordingly to the Use of the TPC for the Analysis"<<endl;
 314     return;
 315   }
 316
 317   TString sIsoMethod="\0",sUEMethod="\0",sBoundaries="\0";
 318
 319   if(fIsoMethod==0)
 320     sIsoMethod = "Cells";
 321   else if(fIsoMethod==1)
 322     sIsoMethod = "Clust";
 323   else if(fIsoMethod==2)
 324     sIsoMethod = "Tracks";
 325
 326   if(fUEMethod==0)
 327     sUEMethod = "PhiBand";
 328   else if(fUEMethod==1)
 329     sUEMethod = "EtaBand";
 330   else if(fUEMethod==2)
 331     sUEMethod = "PerpCones";
 332   else if(fUEMethod==3)
 333     sUEMethod = "FullTPC";
 334
 335   if(fTPC4Iso)
 336     sBoundaries = "TPC Acceptance";
 337   else
 338     sBoundaries = "EMCAL Acceptance";
 339
 340   if(fWho>1 || fWho==-1){
 341     cout<<"Error!!! OutputMode Can Only Be 0: TTree; 1: THnSparse"<<endl;
 342     return;
 343   }
 344   else{
 345     fOutput = new TList();
 346     fOutput->SetOwner();
 347       //Initialize the common Output histograms
 348     switch (fWho)
 349     {
 350       case 0:
 351
 352 // Tree for QA after cluster selection
 353         fOutputQATree = new TTree("OutQATree","OutQATree");
 354         fOutputQATree->Branch("fevents",&fevents);
 355         fOutputQATree->Branch("fNClustersT",&fNClustersT);
 356         fOutputQATree->Branch("fEClustersT",&fEClustersT);
 357         fOutputQATree->Branch("fPtClustersT",&fPtClustersT);
 358         fOutputQATree->Branch("fEtClustersT",&fEtClustersT);
 359         fOutputQATree->Branch("fEtaClustersT",&fEtaClustersT);
 360         fOutputQATree->Branch("fPhiClustersT",&fPhiClustersT);
 361         fOutputQATree->Branch("fM02ClustersT",&fM02ClustersT);
 362
 363         fOutput->Add(fOutputQATree);
 364
 365         fOutputTree = new TTree("OutTree",Form("OutTree; Iso Method %d, UE Method %d, TPC %d, LC %d, Iso Cone %f, CPV #eta %f #phi %f",fIsoMethod,fUEMethod,fTPC4Iso,fisLCAnalysis,fIsoConeRadius,fdetacut,fdphicut));
 366         fOutputTree->Branch("flambda0T",&flambda0T);
 367         fOutputTree->Branch("fEtT",&fEtT);
 368         fOutputTree->Branch("fPtT",&fPtT);
 369         fOutputTree->Branch("fEtisolatedT",&fEtisolatedT);
 370         fOutputTree->Branch("fPtTiso",&fPtisolatedT);
 371         fOutputTree->Branch("fetaT",&fetaT);
 372         fOutputTree->Branch("fphiT",&fphiT);
 373         fOutputTree->Branch("fsumEtisoconeT",&fsumEtisoconeT);
 374         fOutputTree->Branch("fsumEtUE",&fsumEtUE);
 375
 376         fOutput->Add(fOutputTree);
 377
 378
 379         break;
 380       case 1:
 381           //Initialization by Davide;
 382
 383         TString sTitle;
 384         Int_t binTrackMult=100, binPT=70, binM02=100, binETiso=100, binETUE=110, binETisoUE=110, binetacl=140,binphicl=100;
 385
 386         Int_t binMCMotherPDG=50;
 387
 388         Int_t bins[] = {binTrackMult, binPT, binM02, binETiso, binETUE, binETisoUE, binetacl, binphicl};
 389
 390         fNDimensions = sizeof(bins)/sizeof(Int_t);
 391         const Int_t ndims =   fNDimensions;
 392
 393         Double_t xmin[]= {  0.,  0., 0., -10., -10., -10.,-1.0, 1. };
 394
 395         Double_t xmax[]= {1000., 70., 2., 100., 100., 100., 1.0, 3.5,};
 396
 397         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());
 398
 399         fOutputTHnS =  new THnSparseF("fHnOutput",sTitle.Data(), ndims, bins, xmin, xmax);
 400
 401         fOutput->Add(fOutputTHnS);
 402
 403         Int_t binsMC[] = {binTrackMult, binPT, binETiso, binETUE, binMCMotherPDG };
 404
 405         if(fIsMC){
 406
 407           fMCDimensions = sizeof(binsMC)/sizeof(Int_t);
 408             //const Int_t nDimMC = fMCDimensions;
 409
 410           Double_t xminbis[] = {   0.,  0., -10., -10.,    0.};
 411           Double_t xmaxbis[] = {1000., 70., 100., 100., 1000.};
 412
 413           fOutMCTruth = new THnSparseF ("fOutMCTruth","Multiplicity, E_{#gamma}, UE, TruthET,TruthUETE, MomPDG; N_{Tracks}; E_{T}^{#gamma} (GeV/c); p_{T}^{Iso}(GeV/c); E_{T}^{gamma} True; E_{T} ^{UE} True",3,binsMC,xminbis,xmaxbis);
 414
 415           fOutput->Add(fOutMCTruth);
 416           break;
 417         }
 418     }
 419   }
 420
 421   if(fQA){
 422       //Include QA plots to the OutputList //DEFINE BETTER THE BINNING AND THE AXES LIMITS
 423     fTrackMult = new TH1D ("hTrackMult","Tracks multiplicity Distribution",250,0.,1000.);
 424     fTrackMult->Sumw2();
 425     fOutput->Add(fTrackMult);
 426
 427     fTrackMultEMCAL = new TH1D ("hTrackMultEMCAL","Tracks multiplicity Distribution inside EMCAL acceptance",250,0.,1000.);
 428     fTrackMultEMCAL->Sumw2();
 429     fOutput->Add(fTrackMultEMCAL);
 430
 431     fClustMult = new TH1D ("hClustMult","Clusters multiplicity Distribution",1000,0.,1000.);
 432     fClustMult->Sumw2();
 433     fOutput->Add(fClustMult);
 434
 435     fRecoPV = new TH1D("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5);
 436     fRecoPV->Sumw2();
 437     fOutput->Add(fRecoPV);
 438
 439     fPVZBefore = new TH1D ("hPVZDistr","Z Distribution for the Reconstructed Vertex",200,0.,40.);
 440     fPVZBefore->Sumw2();
 441     fOutput->Add(fPVZBefore);
 442
 443     fEtaPhiCell = new TH2D ("hEtaPhiCellActivity","",250,0.,1000., 250, 0., 1000.);
 444     fEtaPhiCell->Sumw2();
 445     fOutput->Add(fEtaPhiCell);
 446
 447     fEtaPhiClus = new TH2D ("hEtaPhiClusActivity","",250,-0.8,0.8, 250, 1.2, 3.4);
 448     fEtaPhiClus->Sumw2();
 449     fOutput->Add(fEtaPhiClus);
 450
 451     fClusEvsClusT = new TH2D ("hEnVSTime","",250,0.,1000., 250,0.,1000.);
 452     fClusEvsClusT->Sumw2();
 453     fOutput->Add(fClusEvsClusT);
 454
 455     fDeltaETAClusTrack = new TH1D("h_Dz","Track-Cluster Dz ",100,-0.05,0.05);
 456     fDeltaETAClusTrack ->Sumw2();
 457     fOutput->Add(fDeltaETAClusTrack);
 458
 459     fDeltaPHIClusTrack = new TH1D("h_Dx","Track-Cluster Dx",100,-0.05,0.05);
 460     fDeltaPHIClusTrack->Sumw2();
 461     fOutput->Add(fDeltaPHIClusTrack);
 462
 463     fDeltaETAClusTrackMatch = new TH1D("h_DzMatch","Track-Cluster Dz matching ",100,-0.05,0.05);
 464     fDeltaETAClusTrackMatch ->Sumw2();
 465     fOutput->Add(fDeltaETAClusTrackMatch);
 466
 467     fDeltaPHIClusTrackMatch = new TH1D("h_DxMatch","Track-Cluster Dx matching",100,-0.05,0.05);
 468     fDeltaPHIClusTrackMatch->Sumw2();
 469     fOutput->Add(fDeltaPHIClusTrackMatch);
 470
 471     fDeltaETAClusTrackVSpT = new TH2D("hTC_Dz","Track-Cluster Dz vs pT of Cluster",100,0.,100.,100,-0.05,0.05);
 472     fDeltaETAClusTrackVSpT->Sumw2();
 473     fOutput->Add(fDeltaETAClusTrackVSpT);
 474
 475     fDeltaPHIClusTrackVSpT = new TH2D("hTC_Dx","Track-Cluster Dx vs pT of Cluster",100,0.,100.,100,-0.05,0.05);
 476     fDeltaPHIClusTrackVSpT->Sumw2();
 477     fOutput->Add(fDeltaPHIClusTrackVSpT);
 478
 479   }
 480     //Initialize only the Common THistos for the Three different output
 481
 482   fVz = new TH1D("hVz_NC","Vertex Z distribution",100,-50.,50.);
 483   fVz->Sumw2();
 484   fOutput->Add(fVz);
 485
 486   fEvents = new TH1D("hEvents_NC","Events",100,0.,100.);
 487   fEvents->Sumw2();
 488   fOutput->Add(fEvents);
 489
 490   fPT = new TH1D("hPt_NC","P_{T} distribution for Neutral Clusters",100,0.,100.);
 491   fPT->Sumw2();
 492   fOutput->Add(fPT);
 493
 494   fE = new TH1D("hE_NC","E distribution for Clusters",200,0.,100.);
 495   fE->Sumw2();
 496   fOutput->Add(fE);
 497
 498   fPtaftTime = new TH1D("hPtaftTime_NC","p_{T} distribution for Clusters after cluster time cut",200,0.,100.);
 499   fPtaftTime->Sumw2();
 500   fOutput->Add(fPtaftTime);
 501
 502   fPtaftTM = new TH1D("hPtaftTM_NC","p_{T} distribution for Neutral Clusters",200,0.,100.);
 503   fPtaftTM->Sumw2();
 504   fOutput->Add(fPtaftTM);
 505
 506   fPtaftFC = new TH1D("hPtaftFC_NC","p_{T} distribution for Clusters after fiducial cut",200,0.,100.);
 507   fPtaftFC->Sumw2();
 508   fOutput->Add(fPtaftFC);
 509
 510   fPtaftM02C = new TH1D("hPtaftM02C_NC","p_{T} distribution for Clusters after shower shape cut",200,0.,100.);
 511   fPtaftM02C->Sumw2();
 512   fOutput->Add(fPtaftM02C);
 513
 514   fClusTime = new TH1D("hClusTime_NC","Time distribution for Clusters",800,-50.,50.);
 515   fClusTime->Sumw2();
 516   fOutput->Add(fClusTime);
 517
 518   fM02 = new TH2D("hM02_NC","M02 distribution for Neutral Clusters vs E",100,0.,100.,500,0.,5.);
 519   fM02->Sumw2();
 520   fOutput->Add(fM02);
 521
 522   fNLM = new TH1D("hNLM_NC","NLM distribution for Neutral Clusters",200,0.,4.);
 523   fNLM->Sumw2();
 524   fOutput->Add(fNLM);
 525
 526   fEtIsoCells = new TH1D("hEtIsoCell_NC","E_{T}^{iso cone} in iso cone distribution for Neutral Clusters with EMCAL Cells",200,-0.25,99.75);
 527   fEtIsoCells->SetXTitle("#Sigma E_{T}^{iso cone} (GeV/c)");
 528   fEtIsoCells->Sumw2();
 529   fOutput->Add(fEtIsoCells);
 530
 531   fEtIsoClust = new TH2D("hEtIsoClus_NC","#Sigma p_{T}^{iso cone} in iso cone distribution for Neutral Clusters with EMCAL Clusters",200,0.,100.,200,0.,100.);
 532   fEtIsoClust->SetYTitle("#Sigma P_{T}^{iso cone} (GeV/c)");
 533   fEtIsoClust->SetXTitle("p_{T}^{clust}");
 534   fEtIsoClust->Sumw2();
 535   fOutput->Add(fEtIsoClust);
 536
 537   fPtIsoTrack = new TH2D("hPtIsoTrack_NC"," #Sigma p_{T}^{iso cone} in iso cone distribution for Neutral Clusters with Tracks",200,0.,100.,200,0.,100.);
 538   fPtIsoTrack->SetYTitle("#Sigma p_{T}^{iso cone} (GeV/c)");
 539   fPtIsoTrack->SetXTitle("p_{T}^{clust}");
 540   fPtIsoTrack->Sumw2();
 541   fOutput->Add(fPtIsoTrack);
 542
 543   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",200,-0.25,99.75);
 544   fPtEtIsoTC->SetXTitle("#Sigma P_{T}^{iso cone} + #Sigma E_{T}^{iso cone} (GeV/c)");
 545   fPtEtIsoTC->Sumw2();
 546   fOutput->Add(fPtEtIsoTC);
 547
 548   fPhiBandUEClust = new TH2D(Form("hPhiBandUE_Cluster"),Form("UE Estimation with Phi Band Clusters"),100,0.,100.,100,0.,100.);
 549   fPhiBandUEClust->SetXTitle("E_{T}");
 550   fPhiBandUEClust->SetYTitle("#Sigma E_{T}^{UE}");
 551   fPhiBandUEClust->Sumw2();
 552   fOutput->Add(fPhiBandUEClust);
 553
 554   fEtaBandUEClust = new TH2D(Form("hEtaBandUE_Cluster"),Form("UE Estimation with Phi Band Clusters"),100,0.,100.,100,0.,100.);
 555   fEtaBandUEClust->SetXTitle("E_{T}");
 556   fEtaBandUEClust->SetYTitle("#Sigma E_{T}^{UE}");
 557   fEtaBandUEClust->Sumw2();
 558   fOutput->Add(fEtaBandUEClust);
 559
 560   fPhiBandUECells = new TH2D(Form("hPhiBandUE_CELLS"),Form("UE Estimation with Phi Band CELLS"),100,0.,100.,100,0.,100.);
 561   fPhiBandUECells->SetXTitle("E_{T}");
 562   fPhiBandUECells->SetYTitle("#Sigma E_{T}^{UE}");
 563   fPhiBandUECells->Sumw2();
 564   fOutput->Add(fPhiBandUECells);
 565
 566   fEtaBandUECells = new TH2D(Form("hEtaBandUE_CELLS"),Form("UE Estimation with Phi Band and CELLS"),100,0.,100.,100,0.,100.);
 567   fEtaBandUECells->SetXTitle("E_{T}");
 568   fEtaBandUECells->SetYTitle("#Sigma E_{T}^{UE}");
 569   fEtaBandUECells->Sumw2();
 570   fOutput->Add(fEtaBandUECells);
 571
 572   fPhiBandUETracks = new TH2D(Form("hPhiBandUE_TPC"),Form("UE Estimation with Phi Band TPC "),100,0.,100.,100,0.,100.);
 573   fPhiBandUETracks->SetXTitle("E_{T}");
 574   fPhiBandUETracks->SetYTitle("#Sigma P_{T}^{UE}");
 575   fPhiBandUETracks->Sumw2();
 576   fOutput->Add(fPhiBandUETracks);
 577
 578   fEtaBandUETracks = new TH2D(Form("hEtaBandUE_TPC"),Form("UE Estimation with Phi Band and TPC"),100,0.,100.,100,0.,100.);
 579   fEtaBandUETracks->SetXTitle("E_{T}");
 580   fEtaBandUETracks->SetYTitle("#Sigma P_{T}^{UE}");
 581   fEtaBandUETracks->Sumw2();
 582   fOutput->Add(fEtaBandUETracks);
 583
 584   fPerpConesUETracks = new TH2D("hConesUE","UE Estimation with Perpendicular Cones in TPC",100,0.,100.,100,0.,100.);
 585   fPerpConesUETracks->SetXTitle("E_{T}");
 586   fPerpConesUETracks->SetYTitle("#Sigma P_{T}^{UE}");
 587   fPerpConesUETracks->Sumw2();
 588   fOutput->Add(fPerpConesUETracks);
 589
 590   fTPCWithoutIsoConeB2BbandUE = new TH2D("hFullTPCUE","UE Estimation with almost Full TPC",100,0.,100.,100,0.,100.);
 591   fPhiBandUEClust->SetXTitle("E_{T}");
 592   fPhiBandUEClust->SetYTitle("#Sigma E_{T}^{UE}");
 593   fTPCWithoutIsoConeB2BbandUE->Sumw2();
 594   fOutput->Add(fTPCWithoutIsoConeB2BbandUE);
 595
 596   fEtIsolatedClust = new TH1D("hEtIsolatedClust","E_{T} distribution for Isolated Photons with clusters; #Sigma E_{T}^{iso cone}<Ethres",200,0.,100.);
 597   fEtIsolatedClust->SetXTitle("E_{T}^{iso}");
 598   fEtIsolatedClust->Sumw2();
 599   fOutput->Add(fEtIsolatedClust);
 600
 601   fPtIsolatedNClust = new TH1D("hEtIsolatedNClust","p_{T} distribution for neutral clusters; #Sigma p_{T}^{iso cone}<Pthres",200,0.,100.);
 602   fPtIsolatedNClust->SetXTitle("p_{T}^{iso}");
 603   fPtIsolatedNClust->Sumw2();
 604   fOutput->Add(fPtIsolatedNClust);
 605
 606   fPtIsolatedNTracks = new TH1D("hEtIsolatedNTracks","p_{T} distribution for neutral clusters; #Sigma p_{T}^{iso cone}<Pthres",200,0.,100.);
 607   fPtIsolatedNTracks->SetXTitle("p_{T}^{iso}");
 608   fPtIsolatedNTracks->Sumw2();
 609   fOutput->Add(fPtIsolatedNTracks);
 610
 611   fEtIsolatedCells = new TH1D("hEtIsolatedCells","E_{T} distribution for Isolated Photons with cells; #Sigma E_{T}^{iso cone}<Ethres",100,0.,100.);
 612   fEtIsolatedCells->SetXTitle("E_{T}^{iso}");
 613   fEtIsolatedCells->Sumw2();
 614   fOutput->Add(fEtIsolatedCells);
 615
 616   fEtIsolatedTracks = new TH1D("hEtIsolatedTracks","E_{T} distribution for Isolated Photons with tracks; #Sigma P_{T}^{iso cone}<Pthres",100,0.,100.);
 617   fEtIsolatedTracks->SetXTitle("E_{T}^{iso}");
 618   fEtIsolatedTracks->Sumw2();
 619   fOutput->Add(fEtIsolatedTracks);
 620
 621   fPtvsM02iso = new TH2D("hPtvsM02iso","p_{T} vs #lambda_{0}^{2} distribution for isolated clusters",200,0.,100.,500,0.,5.);
 622   fPtvsM02iso->SetXTitle("p_{T}^{iso}");
 623   fPtvsM02iso->SetYTitle("#lambda_{0}^{2}");
 624   fOutput->Add(fPtvsM02iso);
 625
 626   fPtvsM02noiso = new TH2D("hPtvsM02noiso","p_{T} vs #lambda_{0}^{2} distribution for non isolated clusters",200,0.,100.,500,0.,5.);
 627   fPtvsM02noiso->SetXTitle("p_{T}^{iso}");
 628   fPtvsM02noiso->SetYTitle("#lambda_{0}^{2}");
 629   fOutput->Add(fPtvsM02noiso);
 630
 631   fTestIndex= new TH2D("hTestIndex","Test index pour cluster",100,0.,100.,100,0.,100.);
 632   fTestIndex->SetXTitle("index");
 633   fTestIndex->SetYTitle("local index");
 634   fTestIndex->Sumw2();
 635   fOutput->Add(fTestIndex);
 636
 637     fTestIndexE= new TH2D("hTestIndexE","Test index vs energy pour cluster",200,0.,100.,100,0.,100.);
 638   fTestIndexE->SetXTitle("cluster energy");
 639   fTestIndexE->SetYTitle("index");
 640   fTestIndexE->Sumw2();
 641   fOutput->Add(fTestIndexE);
 642
 643     fTestLocalIndexE= new TH2D("hTestLocalIndexE","Test local index vs energy for cluster",200,0.,100.,100,0.,100.);
 644   fTestLocalIndexE->SetXTitle("cluster energy");
 645   fTestLocalIndexE->SetYTitle("local index");
 646   fTestLocalIndexE->Sumw2();
 647   fOutput->Add(fTestLocalIndexE);
 648
 649   fTestEnergyCone= new TH2D("hTestEnergyCone","Test energy clusters and tracks in cone",200,0.,100.,200,0.,100.);
 650   fTestEnergyCone->SetXTitle("#sum^{cone} p_{T}^{cluster}");
 651   fTestEnergyCone->SetYTitle("#sum^{cone} p_{T}^{track}");
 652   fTestEnergyCone->Sumw2();
 653   fOutput->Add(fTestEnergyCone);
 654
 655   fTestEtaPhiCone= new TH2D("hTestEtatPhiCone","Test eta phi neutral clusters candidates",200,0.,100.,200,0.,100.);
 656   fTestEtaPhiCone->SetXTitle("phi");
 657   fTestEtaPhiCone->SetYTitle("eta");
 658   fTestEtaPhiCone->Sumw2();
 659   fOutput->Add(fTestEtaPhiCone);
 660
 661   if(fIsMC){
 662       //CREATE THE TH2 specific for the MC Analysis Maybe to be added in the THNSparse, or cloning two or three axes and add the specific MC Truth info
 663   }
 664
 665   PostData(1, fOutput);
 666     //   return;
 667 }
 668
 669   //________________________________________________________________________
 670 Float_t* AliAnalysisTaskEMCALPhotonIsolation::GenerateFixedBinArray(Int_t n, Float_t min, Float_t max) const
 671 {
 672     // Generate the bin array for the ThnSparse
 673
 674   Float_t *bins = new Float_t[n+1];
 675
 676   Float_t binWidth = (max-min)/n;
 677   bins[0] = min;
 678   for (Int_t i = 1; i <= n; i++) {
 679     bins[i] = bins[i-1]+binWidth;
 680   }
 681
 682   return bins;
 683 }
 684
 685   //________________________________________________________________________
 686 void AliAnalysisTaskEMCALPhotonIsolation::ExecOnce()
 687 {
 688     //   Init the analysis.
 689
 690
 691
 692   if (fParticleCollArray.GetEntriesFast()<2) {
 693     AliError(Form("Wrong number of particle collections (%d), required 2",fParticleCollArray.GetEntriesFast()));
 694     return;
 695   }
 696
 697
 698     //    for (Int_t i = 0; i < 2; i++) {
 699     //       AliParticleContainer *contain = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
 700     //       contain->SetClassName("AliEmcalParticle");
 701     //    }
 702
 703
 704
 705   AliAnalysisTaskEmcal::ExecOnce();
 706   if (!fInitialized) {
 707
 708     AliError(Form("AliAnalysisTask not initialized"));
 709     return;
 710   }
 711
 712
 713
 714 }
 715
 716   //______________________________________________________________________________________
 717 Bool_t AliAnalysisTaskEMCALPhotonIsolation::Run()
 718 {
 719     // Run the analysis
 720 fTest1+=1;
 721     //vertex cuts
 722   if (fVertex[2]>10 || fVertex[2]<-10) return kFALSE;
 723
 724   AliParticleContainer *clusters = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
 725
 726 Int_t nbTracksEvent;
 727 nbTracksEvent =InputEvent()->GetNumberOfTracks();
 728 if(nbTracksEvent==0) return kFALSE;
 729
 730     // Fill events number histogram
 731   fEvents->Fill(0);
 732
 733     //Fill Vertex Z histogram
 734   fVz->Fill(fVertex[2]);
 735
 736     // delete output USEFUL LATER FOR CONTAINER CREATION !!
 737     //fOutClusters->Delete();
 738
 739
 740 Int_t index;
 741
 742
 743     //Double_t ETleadingclust = 0., M02leadingcluster = 0., lambda0cluster = 0., phileadingclust = 0., etaleadingclust = 0., ptmc = 0.,mcptsum = 0.;
 744     //Int_t Ntracks;
 745   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 746     //fVevent = dynamic_cast<AliVEvent*>(InputEvent());
 747
 748   if(fIsMC){
 749       //event=fAOD;
 750       //    AliMCEvent *mcEvent;
 751       //    AliMCEventHandler *eventHandler;
 752     AliAODMCHeader *mcHeader;
 753
 754     fAODMCParticles = dynamic_cast <TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 755     mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 756
 757     AnalyzeMC();
 758   }
 759
 760   if (fisLCAnalysis) {
 761
 762       // get the leading particle
 763    AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(clusters->GetLeadingParticle());
 764
 765     if(!emccluster){
 766
 767       AliError(Form("No leading cluster"));
 768       return kFALSE;
 769     }
 770
 771
 772     index = emccluster->IdInCollection();
 773       //add a command to get the index of the leading cluster!
 774     if (!emccluster->IsEMCAL()) return kFALSE;
 775
 776     AliVCluster *coi = emccluster->GetCluster();
 777     if (!coi) return kFALSE;
 778
 779     TLorentzVector vecCOI;
 780     coi->GetMomentum(vecCOI,fVertex);
 781
 782
 783     Double_t coiTOF = coi->GetTOF()*1e9;
 784     if(coiTOF<-30 || coiTOF>30)
 785       return kFALSE;
 786
 787
 788     if(ClustTrackMatching(emccluster))
 789       return kFALSE;
 790
 791     if(!CheckBoundaries(vecCOI))
 792       return kFALSE;
 793
 794     else
 795       FillGeneralHistograms(coi,vecCOI, index);
 796   }
 797   else{
 798       //get the entries of the Cluster Container
 799       //whatever is a RETURN in LCAnalysis here is a CONTINUE,
 800       //since there are more than 1 Cluster per Event
 801        // clusters->ResetCurrentID();
 802
 803    //    AliError("fonctionne bien");
 804         AliEmcalParticle *emccluster=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(0));
 805         index=0;
 806    //    AliError("fonctionne bien");
 807
 808     while(emccluster){
 809
 810       AliVCluster *coi = emccluster->GetCluster();
 811       if(!coi) {
 812        index++;
 813         emccluster = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
 814         continue;
 815       }
 816         //
 817       TLorentzVector vecCOI;
 818       coi->GetMomentum(vecCOI,fVertex);
 819       Double_t coiTOF = coi->GetTOF()*1e9;
 820    //   Double_t coiM02 = coi->GetM02();
 821
 822       FillQAHistograms(coi,vecCOI);
 823         //AliInfo(Form("Cluster number: %d; \t Cluster ToF: %lf ;\tCluster M02:%lf\n",index,coiTOF,coiM02));
 824
 825     //    if(vecCOI.E()<0.3){ // normally already done
 826     //              index++;
 827     //      emccluster=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
 828     //      continue;
 829     //    }
 830
 831       if(!fIsMC){
 832         if(coiTOF<-30 || coiTOF>30){
 833           index++;
 834           emccluster=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
 835           continue;
 836         }
 837       }
 838       else{//different timing cuts for NON CALIBRATED ToF Signal...
 839         if(coiTOF<-570 || coiTOF>630){
 840           index++;
 841           emccluster=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
 842           continue;
 843         }
 844       }
 845       fPtaftTime->Fill(vecCOI.Pt());
 846       if(ClustTrackMatching(emccluster)){
 847         index++;
 848         emccluster = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
 849         continue;
 850       }
 851      fPtaftTM->Fill(vecCOI.Pt());
 852
 853
 854       if(!CheckBoundaries(vecCOI)){
 855         index++;
 856         emccluster = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
 857         continue;
 858       }
 859
 860           fPtaftFC->Fill(vecCOI.Pt());
 861
 862       if(vecCOI.Et()<5.){
 863            index++;
 864         emccluster = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
 865           continue;
 866
 867       }
 868
 869 fTestIndexE->Fill(vecCOI.Pt(),index);
 870
 871         //AliInfo("Passed the CheckBoundaries conditions");
 872
 873  FillGeneralHistograms(coi,vecCOI, index);
 874       index++;
 875       emccluster = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
 876
 877     }
 878
 879   }
 880
 881     //  PostData(1, fOutput);
 882   return kTRUE;
 883 }
 884
 885
 886   //_________________________________________________________________________________
 887 void AliAnalysisTaskEMCALPhotonIsolation::FillQAHistograms(AliVCluster *coi, TLorentzVector vecCOI){
 888
 889         switch(fWho){
 890
 891       case 0:
 892         fevents=0;
 893         fEClustersT=vecCOI.E();
 894         fPtClustersT=vecCOI.Pt();
 895         fEtClustersT=vecCOI.Et();
 896         fEtaClustersT=vecCOI.Eta();
 897         fPhiClustersT=vecCOI.Phi();
 898         fM02ClustersT=coi->GetM02();
 899
 900         fOutputQATree->Fill();
 901
 902         break;
 903
 904       case 1:
 905
 906         break;
 907
 908     }
 909
 910     fPT->Fill(vecCOI.Pt());
 911     fE->Fill(vecCOI.E());
 912     fM02->Fill(vecCOI.E(),coi->GetM02());
 913     fEtaPhiClus->Fill(vecCOI.Eta(),vecCOI.Phi());
 914
 915     Double_t checktof = coi->GetTOF()*1e9;
 916     if(checktof>-30 && checktof<30){
 917     fClusTime->Fill(checktof);
 918    // fPtaftTime->Fill(vecCOI.Pt());
 919
 920   //  if(!ClustTrackMatching(coi)){
 921   //  fPtaftTM->Fill(vecCOI.Pt());
 922
 923     if(CheckBoundaries(vecCOI)){
 924   //  fPtaftFC->Fill(vecCOI.Pt());
 925
 926     Double_t checkM02=coi->GetM02();
 927     if(fM02mincut < checkM02 && checkM02 < fM02maxcut){
 928     fPtaftM02C->Fill(vecCOI.Pt());
 929     }
 930     }
 931    // }
 932 }
 933 }
 934
 935
 936   //__________________________________________________________________________
 937 Bool_t AliAnalysisTaskEMCALPhotonIsolation::ClustTrackMatching(AliEmcalParticle *partC) {
 938     // Check if the cluster match to a track
 939
 940
 941  AliVCluster *cluster = partC->GetCluster();
 942     TLorentzVector nPart;
 943   cluster->GetMomentum(nPart, fVertex);
 944
 945
 946 AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
 947
 948 AliVCluster *clust = partC -> GetCluster();
 949
 950 Int_t nbMObj = partC -> GetNumberOfMatchedObj();
 951
 952 if (nbMObj == 0) return kFALSE;
 953
 954 for(Int_t i=0;i<nbMObj;i++){
 955 Int_t imt = partC->GetMatchedObjId(i);
 956
 957 if(imt<0) continue;
 958
 959 AliEmcalParticle *partT = static_cast<AliEmcalParticle*>(tracks->GetAcceptParticle(imt));
 960 AliVTrack *mt = partT ->GetTrack();
 961
 962 if(!mt) continue;
 963
 964 Double_t deta = 999;
 965 Double_t dphi = 999;
 966
 967  Double_t veta = mt->GetTrackEtaOnEMCal();
 968   Double_t vphi = mt->GetTrackPhiOnEMCal();
 969
 970   Float_t pos[3] = {0};
 971   clust->GetPosition(pos);
 972   TVector3 cpos(pos);
 973   Double_t ceta     = cpos.Eta();
 974   Double_t cphi     = cpos.Phi();
 975   deta=veta-ceta;
 976   dphi=TVector2::Phi_mpi_pi(vphi-cphi);
 977
 978    fDeltaETAClusTrack->Fill(deta);
 979   fDeltaPHIClusTrack->Fill(dphi);
 980 //
 981
 982  if(TMath::Abs(dphi)<fdphicut && TMath::Abs(deta)<fdetacut){
 983     fDeltaETAClusTrackMatch->Fill(deta);
 984     fDeltaPHIClusTrackMatch->Fill(dphi);
 985     return kTRUE;
 986   }
 987
 988 }
 989
 990 return kFALSE;
 991 }
 992
 993
 994
 995   //__________________________________________________________________________
 996 void AliAnalysisTaskEMCALPhotonIsolation::EtIsoCellPhiBand(TLorentzVector c, Float_t &etIso, Float_t &phiBandcells){
 997     // Underlying events study with EMCAL cells in phi band // have to be tested
 998
 999   AliEMCALGeometry* emcalGeom = AliEMCALGeometry::GetInstance();
1000
1001   Float_t sumEnergyPhiBandCells=0., sumEnergyConeCells=0.;
1002
1003
1004     // check the cell corresponding to the leading cluster
1005   Int_t absId = 999;
1006     //maybe best to call it LeadingCellIdinClus or better maxId since it won't be used anywhere else ???
1007   Bool_t cellLeadingClustID = emcalGeom->GetAbsCellIdFromEtaPhi(c.Eta(),c.Phi(),absId);
1008   if(!cellLeadingClustID) return;
1009
1010   else{
1011     Int_t iTower = -1;
1012     Int_t iModule = -1;
1013     Int_t imEta=-1, imPhi=-1;
1014     Int_t iEta =-1, iPhi =-1;
1015
1016     emcalGeom->GetCellIndex(absId,iModule,iTower,imPhi,imEta); // to get the module, the tower, eta and phi for the cell corresponding to the leading cluster
1017     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
1018
1019       // Get the row and the column of the cell corresponding to the leading cluster in EMCAL
1020     Int_t colCellLeadingClust = iEta;
1021     if(iModule % 2) colCellLeadingClust = AliEMCALGeoParams::fgkEMCALCols + iEta ; // if the SM number is even you need to translate to have the one corresponding in EMCAL
1022     Int_t rowCellLeadingClust = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(iModule/2); // to have the corresponding row in EMCAL
1023
1024       // total number or rows and columns in EMCAL
1025     Int_t nTotalRows = AliEMCALGeoParams::fgkEMCALRows*16/3 ;  // 5 + 2/3 supermodules in a row
1026     Int_t nTotalCols = 2*AliEMCALGeoParams::fgkEMCALCols; // 2 supermodules in a column
1027
1028     Int_t nbConeSize = int(fIsoConeRadius/0.0143); //0.0143 tower acceptance
1029
1030       // Get the cells
1031     AliVCaloCells * cells =InputEvent()->GetEMCALCells();
1032
1033       // define the max and min row and column corresponding to the isolation cone around the seed cell from the leading cluster
1034     Int_t iRowMinCone = rowCellLeadingClust-nbConeSize;
1035     if(iRowMinCone<0) iRowMinCone=0;
1036
1037     Int_t iRowMaxCone = rowCellLeadingClust+nbConeSize;
1038     if(iRowMaxCone>AliEMCALGeoParams::fgkEMCALRows) iRowMaxCone=AliEMCALGeoParams::fgkEMCALRows;  // AliEMCALGeoParams::fgkEMCALRows = 24 in a supermodule
1039
1040     Int_t iColMinCone = colCellLeadingClust - nbConeSize;
1041     if(iColMinCone<0) iColMinCone = 0;
1042
1043     Int_t iColMaxCone = colCellLeadingClust+nbConeSize;
1044     if(iColMaxCone>AliEMCALGeoParams::fgkEMCALCols) iColMaxCone=AliEMCALGeoParams::fgkEMCALCols;  // AliEMCALGeoParams::fgkEMCALCols = 48 in a supermodule
1045
1046       // loop on all cells
1047     for(Int_t iCol=0; iCol<nTotalCols; iCol++){
1048       for(Int_t iRow=0; iRow<nTotalRows; iRow++){
1049           // now recover the cell indexes in a supermodule
1050         Int_t iSector = int(iRow/AliEMCALGeoParams::fgkEMCALRows); // check in which SM is the cell
1051         if(iSector==5) continue;  //
1052         Int_t inModule = -1;
1053         Int_t iColLoc  = -1;
1054         if(iCol < AliEMCALGeoParams::fgkEMCALCols){ // if the SM number is odd the column is the one corresponding in the supermodule
1055           inModule = 2*iSector + 1;
1056           iColLoc  = iCol;
1057         }
1058         else if(iCol > AliEMCALGeoParams::fgkEMCALCols - 1){ // if the SM number is even the column isn't the one corresponding in the supermodule
1059           inModule = 2*iSector;
1060           iColLoc  = iCol-AliEMCALGeoParams::fgkEMCALCols;
1061         }
1062
1063         Int_t iRowLoc  = iRow - AliEMCALGeoParams::fgkEMCALRows*iSector ;
1064
1065         if(TMath::Abs(iCol-colCellLeadingClust)<nbConeSize && TMath::Abs(iCol+colCellLeadingClust)>nbConeSize){
1066           if(iRow<iRowMaxCone && iRow>iRowMinCone){
1067             Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(iModule,iRow,iCol);  // verifier les iRow et iCol
1068             sumEnergyPhiBandCells+=cells->GetAmplitude(iabsId); //should be Et
1069           }
1070         }
1071         else if (TMath::Abs(iCol-colCellLeadingClust)>nbConeSize && TMath::Abs(iCol+colCellLeadingClust)<nbConeSize){
1072           Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(inModule,iRowLoc,iColLoc);  // verifier les iRow et iCol
1073           sumEnergyConeCells+=cells->GetAmplitude(iabsId); //should be Et
1074         }
1075       }
1076     }
1077   }
1078   etIso = sumEnergyConeCells;
1079   phiBandcells = sumEnergyPhiBandCells;
1080 }
1081
1082
1083   //__________________________________________________________________________
1084 void AliAnalysisTaskEMCALPhotonIsolation::EtIsoCellEtaBand(TLorentzVector c, Float_t &etIso, Float_t &etaBandcells){
1085     // Underlying events study with EMCAL cell in eta band // have to be tested
1086
1087
1088   AliEMCALGeometry* emcalGeom = AliEMCALGeometry::GetInstance();
1089
1090   Float_t sumEnergyEtaBandCells=0., sumEnergyConeCells=0.;
1091
1092
1093
1094     // check the cell corresponding to the leading cluster
1095   Int_t absId = 999;
1096     //maybe best to call it LeadingCellIdinClus or better maxId since it won't be used anywhere else ???
1097   Bool_t cellLeadingClustID = emcalGeom->GetAbsCellIdFromEtaPhi(c.Eta(),c.Phi(),absId);
1098   if(!cellLeadingClustID) return;
1099
1100   else{
1101
1102     Int_t iTower = -1;
1103     Int_t iModule = -1;
1104     Int_t imEta=-1, imPhi=-1;
1105     Int_t iEta =-1, iPhi =-1;
1106
1107     emcalGeom->GetCellIndex(absId,iModule,iTower,imPhi,imEta); // to get the module, the tower, eta and phi for the cell corresponding to the leading cluster
1108     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
1109
1110       // Get the row and the column of the cell corresponding to the leading cluster in EMCAL
1111     Int_t colCellLeadingClust = iEta;
1112     if(iModule % 2) colCellLeadingClust = AliEMCALGeoParams::fgkEMCALCols + iEta ; // if the SM number is even you need to translate to have the one corresponding in EMCAL
1113     Int_t rowCellLeadingClust = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(iModule/2); // to have the corresponding row in EMCAL
1114
1115       // total number or rows and columns in EMCAL
1116     Int_t nTotalRows = AliEMCALGeoParams::fgkEMCALRows*16/3 ;  // 5 + 2/3 supermodules in a row
1117     Int_t nTotalCols = 2*AliEMCALGeoParams::fgkEMCALCols; // 2 supermodules in a column
1118
1119     Int_t nbConeSize = int(fIsoConeRadius/0.0143); //0.0143 tower acceptance
1120
1121       // Get the cells
1122     AliVCaloCells * cells =InputEvent()->GetEMCALCells();
1123
1124       // define the max and min row and column corresponding to the isolation cone around the seed cell from the leading cluster
1125     Int_t iRowMinCone = rowCellLeadingClust-nbConeSize;
1126     if(iRowMinCone<0) iRowMinCone=0;
1127
1128     Int_t iRowMaxCone = rowCellLeadingClust+nbConeSize;
1129     if(iRowMaxCone>AliEMCALGeoParams::fgkEMCALRows) iRowMaxCone=AliEMCALGeoParams::fgkEMCALRows;  // AliEMCALGeoParams::fgkEMCALRows = 24 in a supermodule
1130
1131     Int_t iColMinCone = colCellLeadingClust-nbConeSize;
1132     if(iColMinCone<0) iColMinCone = 0;
1133
1134     Int_t iColMaxCone = colCellLeadingClust+nbConeSize;
1135     if(iColMaxCone>AliEMCALGeoParams::fgkEMCALCols) iColMaxCone=AliEMCALGeoParams::fgkEMCALCols;  // AliEMCALGeoParams::fgkEMCALCols = 48 in a supermodule
1136
1137       // loop on all cells
1138     for(Int_t iCol=0; iCol<nTotalCols; iCol++)
1139     {
1140       for(Int_t iRow=0; iRow<nTotalRows; iRow++)
1141       {
1142           // now recover the cell indexes in a supermodule
1143         Int_t iSector = int(iRow/AliEMCALGeoParams::fgkEMCALRows); // check in which SM is the cell
1144         if(iSector==5) continue;  //
1145         Int_t inModule = -1;
1146         Int_t iColLoc  = -1;
1147         if(iCol < AliEMCALGeoParams::fgkEMCALCols){ // if the SM number is odd the column is the one corresponding in the supermodule
1148           inModule = 2*iSector + 1;
1149           iColLoc  = iCol;
1150         }
1151         else if(iCol > AliEMCALGeoParams::fgkEMCALCols - 1){ // if the SM number is even the column isn't the one corresponding in the supermodule
1152           inModule = 2*iSector;
1153           iColLoc  = iCol-AliEMCALGeoParams::fgkEMCALCols;
1154         }
1155
1156         Int_t iRowLoc  = iRow - AliEMCALGeoParams::fgkEMCALRows*iSector ;
1157
1158         if(TMath::Abs(iCol-colCellLeadingClust)<nbConeSize && TMath::Abs(iCol+colCellLeadingClust)>nbConeSize){
1159           if(iCol<iColMaxCone && iCol>iColMinCone){
1160             Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(iModule,iRow,iCol);  // verifier les iRow et iCol
1161             sumEnergyEtaBandCells+=cells->GetAmplitude(iabsId); //should be Et
1162           }
1163         }
1164         else if (TMath::Abs(iCol-colCellLeadingClust)>nbConeSize && TMath::Abs(iCol+colCellLeadingClust)<nbConeSize){
1165           Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(inModule,iRowLoc,iColLoc);  // verifier les iRow et iCol
1166           sumEnergyConeCells+=cells->GetAmplitude(iabsId); //should be Et
1167         }
1168       }
1169     }
1170   }
1171   etIso = sumEnergyConeCells;
1172   etaBandcells = sumEnergyEtaBandCells;
1173 }
1174
1175
1176   //__________________________________________________________________________
1177 void AliAnalysisTaskEMCALPhotonIsolation::EtIsoClusPhiBand(TLorentzVector c, Float_t &ptIso, Float_t &phiBandclus, Int_t index){
1178     // Underlying events study with clusters in phi band
1179
1180   Float_t sumEnergyPhiBandClus=0., sumEnergyConeClus=0., sumpTConeCharged=0.;
1181
1182     //needs a check on the same cluster
1183   AliParticleContainer *clusters = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
1184   clusters->ResetCurrentID();
1185   Int_t localIndex=0;
1186   AliEmcalParticle *clust = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1187
1188   while(clust){ //check the position of other clusters in respect to the leading cluster
1189
1190     if(localIndex==index){
1191       localIndex++;
1192       clust = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1193       continue;
1194     }
1195     else{
1196       localIndex++;
1197
1198       AliVCluster *cluster= clust->GetCluster();
1199       if(!cluster){
1200         clust = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1201         continue;
1202       }
1203
1204       TLorentzVector nClust; //STILL NOT INITIALIZED
1205       cluster->GetMomentum(nClust,fVertex);
1206       Float_t phiClust =nClust.Phi();
1207       Float_t etaClust= nClust.Eta();
1208
1209       Double_t clustTOF = cluster->GetTOF()*1e9;
1210         if(clustTOF<-30 || clustTOF>30){
1211           clust=static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1212           continue;
1213         }
1214
1215         if(ClustTrackMatching(clust)){
1216           clust=static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1217           continue;
1218         }
1219
1220         //redefine phi/c.Eta() from the cluster we passed to the function
1221
1222       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
1223
1224       if(radius>fIsoConeRadius){ // the cluster is outside the isolation cone in this case study UE
1225
1226           // actually phi band here
1227         if(TMath::Abs(etaClust - c.Eta()) < fIsoConeRadius){
1228           sumEnergyPhiBandClus += nClust.Pt();
1229         }
1230       }
1231       else // if the cluster is in the isolation cone, add the cluster pT
1232         sumEnergyConeClus += nClust.Et();
1233
1234       clust = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1235     }
1236   }
1237
1238   fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1239    // name hard coded to use the defined tracks for analysis
1240
1241        if (!fTracksAna) {
1242       AliError(Form("Could not retrieve tracks !"));
1243       return;
1244     }
1245     const Int_t nbTracks = fTracksAna->GetEntries();
1246     Int_t iTracks = 0;
1247
1248 //
1249      while(iTracks<nbTracks){
1250          AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1251          if(!track){
1252             AliError(Form("No tracks in collection"));
1253               iTracks++;
1254               continue;
1255          }
1256       //CHECK IF TRACK IS IN BOUNDARIES
1257     Float_t phiTrack = track->Phi();
1258     Float_t etaTrack = track->Eta();
1259
1260       Float_t  radius = TMath::Sqrt(TMath::Power(phiTrack - c.Phi(),2)+TMath::Power(etaTrack - c.Eta(),2));
1261
1262       if(radius<fIsoConeRadius){ // if tracks are outside the isolation cone study
1263         sumpTConeCharged+=track->Pt(); // should not double count if the track matching is already done
1264       }
1265     iTracks++;
1266  } // end of tracks loop
1267
1268   fTestEnergyCone->Fill(sumEnergyConeClus,sumpTConeCharged);
1269
1270
1271   ptIso = sumEnergyConeClus + sumpTConeCharged;
1272   phiBandclus = sumEnergyPhiBandClus;
1273 }
1274
1275
1276   //__________________________________________________________________________
1277 void AliAnalysisTaskEMCALPhotonIsolation::EtIsoClusEtaBand(TLorentzVector c, Float_t &ptIso, Float_t &etaBandclus, Int_t index){
1278     // Underlying events study with clusters in eta band
1279
1280
1281
1282   Float_t sumEnergyEtaBandClus =0., sumEnergyConeClus=0., sumpTConeCharged=0;
1283   Double_t clustTOF=0;
1284
1285   AliParticleContainer *clusters = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
1286
1287 //  clusters->ResetCurrentID();
1288   Int_t localIndex=0;
1289   AliEmcalParticle *clust = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1290
1291   while(clust){ //check the position of other clusters in respect to the leading cluster
1292
1293     if(localIndex==index){
1294       localIndex++;
1295       clust = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1296       continue;
1297     }
1298
1299     else{
1300       localIndex++;
1301       AliVCluster *cluster= clust->GetCluster();
1302       if(!cluster){
1303         clust = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1304         continue;
1305       }
1306
1307       TLorentzVector nClust; //STILL NOT INITIALIZED
1308       cluster->GetMomentum(nClust,fVertex);
1309
1310       Float_t phiClust =nClust.Phi();
1311       Float_t etaClust= nClust.Eta();
1312       Float_t eTcluster=0;
1313
1314
1315       clustTOF = cluster->GetTOF()*1e9;
1316         if(clustTOF<-30 || clustTOF>30){
1317           clust=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1318           continue;
1319         }
1320
1321        if(ClustTrackMatching(clust)){
1322           clust=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1323           continue;
1324        }
1325         //redefine phi/c.Eta() from the cluster we passed to the function
1326
1327         // define the radius between the leading cluster and the considered cluster
1328       Float_t  radius = TMath::Sqrt(TMath::Power(phiClust-c.Phi(),2)+TMath::Power(etaClust-c.Eta(),2));
1329
1330
1331
1332       if(radius>fIsoConeRadius){ // the cluster is outside the isolation cone in this case study UE
1333
1334           // actually eta band here
1335         if(TMath::Abs(etaClust - c.Eta()) < fIsoConeRadius){
1336           sumEnergyEtaBandClus += nClust.Et();
1337         }
1338       }
1339       else if(radius<fIsoConeRadius && radius!=0.){  // if the cluster is in the isolation cone, add the cluster pT
1340       eTcluster=nClust.Pt();
1341
1342
1343         sumEnergyConeClus += nClust.Pt();
1344                 fTestEtaPhiCone->Fill(c.Eta(),c.Phi());
1345                  fTestIndex->Fill(index,localIndex);
1346
1347         fTestLocalIndexE->Fill(eTcluster,localIndex);
1348     }
1349       clust = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1350     }
1351   } // end of clusters loop
1352
1353
1354
1355    fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1356    // name hard coded to use the defined tracks for analysis
1357
1358        if (!fTracksAna) {
1359       AliError(Form("Could not retrieve tracks !"));
1360       return;
1361     }
1362     const Int_t nbTracks = fTracksAna->GetEntries();
1363     Int_t iTracks = 0;
1364
1365 //
1366      while(iTracks<nbTracks){
1367          AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1368          if(!track){
1369             AliError(Form("No tracks in collection"));
1370               iTracks++;
1371               continue;
1372          }
1373       //CHECK IF TRACK IS IN BOUNDARIES
1374     Float_t phiTrack = track->Phi();
1375     Float_t etaTrack = track->Eta();
1376
1377       Float_t  radius = TMath::Sqrt(TMath::Power(phiTrack - c.Phi(),2)+TMath::Power(etaTrack - c.Eta(),2));
1378
1379       if(radius<fIsoConeRadius){ // if tracks are outside the isolation cone study
1380         sumpTConeCharged+=track->Pt(); // should not double count if the track matching is already done
1381       }
1382     iTracks++;
1383  } // end of tracks loop
1384
1385   fTestEnergyCone->Fill(sumEnergyConeClus,sumpTConeCharged);
1386
1387   ptIso = sumEnergyConeClus + sumpTConeCharged;
1388   etaBandclus = sumEnergyEtaBandClus;
1389 }
1390
1391
1392
1393   //__________________________________________________________________________
1394 void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackPhiBand(TLorentzVector c, Float_t &ptIso, Float_t &phiBandtrack){
1395     // Underlying events study with tracks in phi band in EMCAL acceptance
1396
1397     //INSERT BOUNDARIES ACCORDING TO THE FLAG WE WANT TO USE
1398   Float_t sumpTConeCharged=0., sumpTPhiBandTrack=0.;
1399   Float_t minPhi= 0., maxPhi= 2*TMath::Pi(), minEta = -0.9, maxEta= 0.9;
1400
1401   if(!fTPC4Iso)
1402     {
1403     minEta = -0.7;
1404     maxEta = 0.7;
1405     minPhi = 1.4;
1406     maxPhi = TMath::Pi();
1407     }
1408
1409
1410     fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1411
1412     if (!fTracksAna)
1413     {
1414       AliError(Form("Could not retrieve tracks !"));
1415       return;
1416     }
1417     const Int_t nbTracks = fTracksAna->GetEntries();
1418     Int_t iTracks = 0;
1419
1420
1421   while(iTracks<nbTracks)
1422   {
1423       AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1424   if(!track)
1425     {
1426         AliError(Form("No tracks in collection"));
1427         iTracks++;
1428         continue;
1429     }
1430       //CHECK IF TRACK IS IN BOUNDARIES
1431       Float_t phiTrack = track->Phi();
1432       Float_t etaTrack = track->Eta();
1433       // define the radius between the leading cluster and the considered cluster
1434       //redefine phi/c.Eta() from the cluster we passed to the function
1435       if(phiTrack < maxPhi && phiTrack > minPhi && etaTrack < maxEta && etaTrack > minEta)
1436       {
1437             Float_t  radius = TMath::Sqrt(TMath::Power(phiTrack - c.Phi(),2)+TMath::Power(etaTrack - c.Eta(),2));
1438
1439             if(radius>fIsoConeRadius)
1440             { // if tracks are outside the isolation cone study
1441
1442           // actually phi band here --- ADD Boundaries conditions
1443           if(TMath::Abs(etaTrack - c.Eta()) < fIsoConeRadius)
1444             {
1445                 sumpTPhiBandTrack += track->Pt();
1446             }
1447             }
1448           else
1449             sumpTConeCharged+=track->Pt(); // should not double count if the track matching is already done
1450         }
1451           iTracks++;
1452     }
1453   ptIso = sumpTConeCharged;
1454   phiBandtrack = sumpTPhiBandTrack;
1455 }
1456
1457
1458   //__________________________________________________________________________
1459 void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackEtaBand(TLorentzVector c, Float_t &ptIso, Float_t &etaBandtrack){
1460     // Underlying events study with tracks in eta band in EMCAL acceptance
1461
1462     //INSERT BOUNDARIES ACCORDING TO THE FLAG WE WANT TO USE
1463   Float_t sumpTConeCharged=0., sumpTEtaBandTrack=0.;
1464   Float_t minPhi= 0., maxPhi= 2*TMath::Pi(), minEta = -0.9, maxEta= 0.9;
1465
1466   if(!fTPC4Iso){
1467     minEta = -0.7;
1468     maxEta = 0.7;
1469     minPhi = 1.4;
1470     maxPhi = TMath::Pi();
1471   }
1472
1473   fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1474
1475   if (!fTracksAna)
1476     {
1477       AliError(Form("Could not retrieve tracks !"));
1478       return;
1479     }
1480     const Int_t nbTracks = fTracksAna->GetEntries();
1481     Int_t iTracks = 0;
1482
1483
1484   while(iTracks<nbTracks)
1485     {
1486      AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1487
1488   if(!track)
1489     {
1490         AliError(Form("No tracks in collection"));
1491         iTracks++;
1492         continue;
1493     }
1494
1495     Float_t phiTrack = track->Phi();
1496     Float_t etaTrack = track->Eta();
1497       //redefine phi/c.Eta() from the cluster we passed to the function
1498     if(phiTrack < maxPhi && phiTrack > minPhi && etaTrack < maxEta && etaTrack > minEta)
1499         {
1500             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
1501
1502             if(radius>fIsoConeRadius)
1503                 { // if tracks are outside the isolation cone study UE
1504
1505           // actually eta band here --- ADD Boundaries conditions
1506                 if(TMath::Abs(phiTrack - c.Phi()) < fIsoConeRadius)
1507                 {
1508                     sumpTEtaBandTrack += track->Pt();
1509                 }
1510                 }
1511                 else sumpTConeCharged += track->Pt(); // should not double count if the track matching is already done
1512         }
1513      iTracks++;
1514     }
1515
1516   ptIso = sumpTConeCharged;
1517   etaBandtrack = sumpTEtaBandTrack;
1518 }
1519
1520
1521   //__________________________________________________________________________
1522 void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackOrthCones(TLorentzVector c, Float_t &ptIso, Float_t &cones){
1523     // Underlying events study with tracks in orthogonal cones in TPC
1524
1525   Float_t sumpTConeCharged=0., sumpTPerpConeTrack=0.;
1526
1527   Float_t etaClus = c.Eta();
1528   Float_t phiClus = c.Phi();
1529   Float_t phiCone1 = phiClus - TMath::PiOver2();
1530   Float_t phiCone2 = phiClus + TMath::PiOver2();
1531
1532   if (phiCone1 < 0.) phiCone1 += 2*TMath::Pi();
1533
1534   fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1535
1536     if (!fTracksAna)
1537     {
1538       AliError(Form("Could not retrieve tracks !"));
1539       return;
1540     }
1541
1542     const Int_t nbTracks = fTracksAna->GetEntries();
1543     Int_t iTracks = 0;
1544
1545   while(iTracks<nbTracks)
1546   {
1547
1548   AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1549
1550   if(!track)
1551     {
1552         AliError(Form("No tracks in collection"));
1553         iTracks++;
1554         continue;
1555     }
1556
1557     Float_t phiTrack = track->Phi();
1558     Float_t etaTrack = track->Eta();
1559     Float_t dist2Clust = TMath::Sqrt(TMath::Power(etaTrack-etaClus, 2)+TMath::Power(phiTrack-phiClus, 2));
1560
1561     if (dist2Clust<fIsoConeRadius) sumpTConeCharged += track->Pt(); // tracks are inside the isolation cone
1562
1563     else
1564     {//tracks outside the IsoCone
1565         //Distances from the centres of the two Orthogonal Cones
1566       Float_t dist2Cone1 = TMath::Sqrt(TMath::Power(etaTrack-etaClus, 2)+TMath::Power(phiTrack-phiCone1, 2));
1567       Float_t dist2Cone2 = TMath::Sqrt(TMath::Power(etaTrack-etaClus, 2)+TMath::Power(phiTrack-phiCone2, 2));
1568
1569         //Is the Track Inside one of the two Cones ->Add to UE
1570       if((dist2Cone1 < fIsoConeRadius) || (dist2Cone2 < fIsoConeRadius)) sumpTPerpConeTrack += track->Pt();
1571
1572     }
1573
1574     iTracks++;
1575
1576   }
1577
1578   ptIso = sumpTConeCharged;
1579   cones = sumpTPerpConeTrack;
1580 }
1581
1582   //__________________________________________________________________________
1583 void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackFullTPC(TLorentzVector c, Float_t &ptIso, Float_t &full){
1584     // Underlying events study with tracks in full TPC except back to back bands
1585
1586   Float_t sumpTConeCharged=0., sumpTTPCexceptB2B=0.;
1587
1588   fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1589
1590     if (!fTracksAna)
1591     {
1592       AliError(Form("Could not retrieve tracks !"));
1593       return;
1594     }
1595
1596     const Int_t nbTracks = fTracksAna->GetEntries();
1597     Int_t iTracks = 0;
1598
1599   while(iTracks<nbTracks)
1600   {
1601
1602   AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1603
1604   if(!track)
1605     {
1606         AliError(Form("No tracks in collection"));
1607         iTracks++;
1608         continue;
1609     }
1610
1611     Float_t phiTrack = track->Phi();
1612     Float_t etaTrack = track->Eta();
1613       //redefine phi/c.Eta() from the cluster we passed to the function
1614     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
1615
1616     if(radius>fIsoConeRadius)
1617     { // if tracks are outside the isolation cone study UE
1618       Double_t dphiUp = c.Phi() + TMath::Pi() - fIsoConeRadius;
1619       Double_t dphiDown = c.Phi() + TMath::Pi() + fIsoConeRadius;
1620         // TPC except B2B
1621       if(phiTrack < dphiDown && phiTrack> dphiUp) sumpTTPCexceptB2B += track->Pt();
1622
1623     }
1624
1625     else sumpTConeCharged += track->Pt(); // should not double count if the track matching is already done
1626
1627     iTracks++;
1628
1629   }
1630
1631   ptIso = sumpTConeCharged;
1632   full = sumpTTPCexceptB2B;
1633 }
1634
1635   //__________________________________________________________________________
1636 Bool_t AliAnalysisTaskEMCALPhotonIsolation::CheckBoundaries(TLorentzVector vecCOI){
1637     // Check if the cone around the considered cluster is in EMCAL acceptance
1638     //AliInfo("Inside CheckBoundaries\n");
1639
1640   Float_t minPhiBound= 1.5 , minEtaBound= -0.5, maxPhiBound= TMath::Pi()-0.1, maxEtaBound= 0.5;
1641   Bool_t isINBoundaries;
1642
1643 if(!fTPC4Iso)
1644     {
1645     minPhiBound = 1.4+fIsoConeRadius;
1646     maxPhiBound = 3.15-fIsoConeRadius; // normally 110° but shorter cut to avoid EMCAL border
1647     minEtaBound = -0.67+fIsoConeRadius; // ""
1648     maxEtaBound = 0.67-fIsoConeRadius; // ""
1649
1650  //   minPhiBound = 1.8; //to be changed with fIsoConeR
1651  //   maxPhiBound = 2.75;
1652  //   minEtaBound = -0.27;
1653  //   maxEtaBound = 0.27;
1654 }
1655
1656
1657   if(vecCOI.Eta() > maxEtaBound || vecCOI.Eta() < minEtaBound || vecCOI.Phi() > maxPhiBound || vecCOI.Phi() <minPhiBound)
1658     isINBoundaries=kFALSE;
1659   else
1660     isINBoundaries=kTRUE;
1661
1662
1663   return isINBoundaries;
1664 }
1665
1666   //_________________________________________________________________________
1667 void AliAnalysisTaskEMCALPhotonIsolation::LookforParticle(Int_t clusterlabel, Double_t energyCLS, Double_t phiCLS, Double_t etaCLS, Double_t time, Double_t ss, Int_t ncells){
1668
1669     //time and ncells not used for the moment
1670
1671     //  AliInfo("Inside AnalyzeMC");
1672   if (!fIsMC)
1673   {
1674     cout<<"not a montecarlo run!!!!!!"<<endl;
1675     return;
1676   } //AliInfo(Form("It's a MC analysis %e",fAODMCParticles));
1677   if(!fStack && !fAODMCParticles){
1678     cout<<"No Particle Stack !!!!!"<<endl;
1679     return;
1680   }
1681     //AliInfo("there's a List of particles");
1682     //DO THIS ALSO FOR ESDs
1683
1684   if(fAODMCParticles->GetEntries() < 1){
1685     cout<<"number of tracks insufficient"<<endl;
1686     return;
1687   }
1688
1689
1690   Int_t ndimsMCmix = fMCQAdim;
1691
1692
1693   Double_t outputvalueMCmix[ndimsMCmix];
1694     //cout<<"dimensions of the array: "<<ndimsMCmix<<endl;
1695
1696
1697   Int_t npart=fAODMCParticles->GetEntries();
1698     //cout<<"Number of particles in the event: "<<npart<<endl;
1699
1700   AliAODMCParticle *particle2Check, *MomP2Check;
1701   // *partMC;
1702
1703   Int_t clustPDG, p2clabel;
1704   Double_t enTrue,phiTrue, etaTrue;
1705   Double_t dPhi,dEta ;
1706   bool found=kFALSE;
1707   for(int b=0; b<npart && found!= kTRUE ;b++){
1708     particle2Check = static_cast<AliAODMCParticle*>(fAODMCParticles->At(b));
1709     p2clabel = particle2Check->Label();
1710
1711     if(clusterlabel==p2clabel){
1712       found=kTRUE;
1713       clustPDG = particle2Check->GetPdgCode();
1714       int mom2checkidx = particle2Check->GetMother();
1715       MomP2Check = static_cast<AliAODMCParticle*>(fAODMCParticles->At(mom2checkidx));
1716         //if(energyCLS>=40.)
1717           //cout<<"PDG associated: "<<clustPDG<<" Mother PDG: "<<MomP2Check->GetPdgCode()<<endl;
1718         if(clustPDG==22 || (TMath::Abs(clustPDG)==11 && MomP2Check->GetPdgCode()==22)) //continue;
1719         {
1720           phiTrue = particle2Check->Phi();
1721           etaTrue = particle2Check->Eta();
1722           enTrue = particle2Check->E()*TMath::Sin(particle2Check->Theta());
1723             //if(energyCLS>=40.)
1724               //cout<<"Energy of the single particle associated with the cluster: "<<enTrue<<endl;
1725             if(clustPDG==22){
1726               if(MomP2Check->GetPdgCode()==111){
1727
1728                 Int_t idxdaug1 = MomP2Check->GetFirstDaughter();
1729                 if (idxdaug1<npart){
1730                   if(idxdaug1==clusterlabel){
1731                     Int_t idxdaug2 = MomP2Check->GetLastDaughter();
1732                     if(idxdaug2<npart){
1733                       AliAODMCParticle *daug2 = static_cast<AliAODMCParticle*>(fAODMCParticles->At(idxdaug2));
1734                       if(daug2->GetPdgCode()==22 && (daug2->Phi()-phiTrue)<0.2 && (daug2->Eta()-etaTrue)<0.2){
1735                           //if(energyCLS >= 40.){
1736                             //cout<<"CASE1\nPDG of the other particle VERY close: "<<daug2->GetPdgCode()<<" with Label: "<<daug2->Label();
1737                             //cout<<" Energy of the other particle VERY close: "<<daug2->E()*TMath::Sin(daug2->Theta())<<endl;
1738                           //}
1739                         enTrue += daug2->E()*TMath::Sin(daug2->Theta());
1740                       }
1741                     }
1742                   }
1743                   else{
1744                     AliAODMCParticle *daug1 = static_cast<AliAODMCParticle*>(fAODMCParticles->At(idxdaug1));
1745
1746                     if(daug1->GetPdgCode()==22 && (daug1->Phi()-phiTrue)<0.2 && (daug1->Eta()-etaTrue)<0.2){
1747                         //if(energyCLS >= 40.){
1748                           //cout<<"CASE2\nPDG of the other particle VERY close: "<<daug1->GetPdgCode()<<" with Label: "<<daug1->Label();
1749                           //cout<<" Energy of the other particle VERY close: "<<daug1->E()*TMath::Sin(daug1->Theta())<<endl;
1750                         //}
1751                       enTrue += daug1->E()*TMath::Sin(daug1->Theta());
1752                     }
1753                   }
1754                 }
1755               }
1756             }
1757             else{
1758               Int_t firstidx=MomP2Check->GetFirstDaughter();
1759               if(firstidx< npart){
1760                 if(firstidx==clusterlabel){
1761                   Int_t lastidx=MomP2Check->GetLastDaughter();
1762                   if(lastidx<npart){
1763                     AliAODMCParticle *last=static_cast<AliAODMCParticle*>(fAODMCParticles->At(lastidx));
1764                     if((last->Phi()-phiTrue)<0.03 && (last->Eta()-etaTrue)<0.02){
1765                         //if(energyCLS >= 40.){
1766                           //cout<<"CASE3\nPDG of the other particle VERY close: "<<last->GetPdgCode()<<" with Label: "<<last->Label();
1767                           //cout<<" Energy of the other particle VERY close: "<<last->E()*TMath::Sin(last->Theta())<<endl;
1768                         //}
1769                       enTrue += last->E()*TMath::Sin(last->Theta());
1770                     }
1771                   }
1772                 }
1773                 else{
1774                   AliAODMCParticle *first=static_cast<AliAODMCParticle*>(fAODMCParticles->At(firstidx));
1775                   if((first->Phi()-phiTrue)<0.03 && (first->Eta()-etaTrue)<0.02){
1776                       //if(energyCLS >= 40.){
1777                         //cout<<"CASE4\nPDG of the other particle VERY close: "<<first->GetPdgCode()<<" with Label: "<<first->Label();
1778                         //cout<<" Energy of the other particle VERY close: "<<first->E()*TMath::Sin(first->Theta())<<endl;
1779                       //}
1780                     enTrue += first->E()*TMath::Sin(first->Theta());
1781                   }
1782                 }
1783               }
1784               Int_t idxgrandma = MomP2Check->GetMother();
1785               AliAODMCParticle *grandma=static_cast<AliAODMCParticle*>(fAODMCParticles->At(idxgrandma));
1786               if(grandma->GetPdgCode()==111){
1787                   //if(energyCLS >= 40.){
1788                     //cout<<"Energy of the pi0 grandmother: "<<grandma->E()*TMath::Sin(grandma->Theta())<<endl;
1789                   //}
1790                 Int_t idxaunt = grandma->GetFirstDaughter();
1791                 if(idxaunt<npart){
1792                   if(idxaunt == mom2checkidx){
1793                     Int_t auntid = grandma->GetLastDaughter();
1794                     if(auntid<npart){
1795                       AliAODMCParticle *lastaunt=static_cast<AliAODMCParticle*>(fAODMCParticles->At(auntid));
1796                       if((lastaunt->Phi()-phiTrue)<0.03 && (lastaunt->Eta()-etaTrue)<0.02){
1797                           //if(energyCLS >= 40.){
1798                             //cout<<"CASE5\nPDG of the other particle VERY close: "<<lastaunt->GetPdgCode()<<" with Label: "<<lastaunt->Label();
1799                             //cout<<" Energy of the other particle VERY close: "<<lastaunt->E()*TMath::Sin(lastaunt->Theta())<<endl;
1800                           //}
1801                         enTrue += lastaunt->E()*TMath::Sin(lastaunt->Theta());
1802                       }
1803                     }
1804                   }
1805                   else{
1806                     AliAODMCParticle *aunt =static_cast<AliAODMCParticle*>(fAODMCParticles->At(idxaunt));
1807                     if(aunt->GetPdgCode()==22 && (aunt->Phi()-phiTrue)<0.03 && (aunt->Eta()-etaTrue)<0.02){
1808                         //if(energyCLS >= 40.){
1809                           //cout<<"CASE6\nPDG of the other particle VERY close: "<<aunt->GetPdgCode()<<" with Label: "<<aunt->Label();
1810                           //cout<<" Energy of the other particle VERY close: "<<aunt->E()*TMath::Sin(aunt->Theta())<<endl;
1811                         //}
1812                       enTrue += aunt->E()*TMath::Sin(aunt->Theta());
1813                     }
1814                   }
1815                 }
1816               }
1817             }
1818
1819           dPhi = phiCLS-phiTrue;
1820           dEta = etaCLS-etaTrue;
1821
1822             //      if(fcount==388)
1823             //        AliInfo(Form("Found Particle with same label as cluster !!!! at position %d",b));
1824             //      if(fcount==388){
1825             //        AliInfo(Form(""));
1826             //        particle2Check->Print();
1827             //        cout<<"Energy of the Particle: "<<enTrue<<"  Mother PDG: "<<MomP2Check->GetPdgCode()<<"  Eta: "<<etaTrue<<"  Phi: "<<phiTrue<<endl;
1828             //if(energyCLS >= 40.){
1829               //cout<<"Transverse Energy of all the Particle VERY CLOSE TO THe ClusterLabel Particle: "<<enTrue<<endl;
1830               //cout<<endl;
1831             //}
1832           outputvalueMCmix[0] = energyCLS;
1833           outputvalueMCmix[1] = ss;
1834           outputvalueMCmix[2] = clustPDG;
1835           outputvalueMCmix[3] = MomP2Check->GetPdgCode();
1836           outputvalueMCmix[4] = enTrue;
1837             //outputvalueMCmix[5] = dPhi;
1838             //outputvalueMCmix[6] = dEta;
1839             //outputvalueMCmix[7] = ncells;
1840           fOutClustMC->Fill(outputvalueMCmix);
1841         }
1842         //      }
1843         //fPDGM02->Fill(clustPDG);
1844         //fEtrueEclustM02->Fill(energyCLS,enTrue);
1845         //fDphiDetaM02->Fill(dEta,dPhi);
1846         //fMomPDGM02->Fill(MomP2Check->GetPdgCode());
1847
1848         //if(TMath::Abs(enTrue-energyCLS)>0.2){
1849         //cout<<"Time of the cluster with energy mismatch: "<<time<<" energy of the cluster: "<<energyCLS<<" true energy: "<<enTrue<<" PDG "<<clustPDG<<" mother of the particle: "<<MomP2Check->GetPdgCode()<<endl;
1850         //fTvsE_MismatchEM02->Fill(enTrue,time);
1851         //break;
1852         //}
1853     }
1854   }
1855   if(found==kFALSE)
1856     printf("not a particle!!! Look at the STACK DOWN HERE!!!!\n\n");
1857   /*clustPDG=0;
1858    dPhi=5.;
1859    dEta=5.;
1860    enTrue = -1.;
1861    for(int b=0; b<npart; b++){
1862
1863    partMC=static_cast<AliAODMCParticle*>(fAODMCParticles->At(b));
1864    cout<<"particle "<<b<<endl;
1865    partMC->Print();
1866
1867    }
1868    fPDGM02->Fill(clustPDG);
1869    fEtrueEclustM02->Fill(energyCLS,enTrue);
1870    fDphiDetaM02->Fill(dEta,dPhi);
1871    return;
1872    }*/
1873
1874
1875
1876     //cout<<"EnergyT: "<<particle2Check->E()*TMath::Sin(particle2Check->Theta())<<"\tPDGCode: "<<particle2Check->GetPdgCode()<<"\tMotherPDG: "<<MomP2Check->GetPdgCode()<<"\tEta: "<<particle2Check->Eta()<<"\tPhi: "<<particle2Check->Phi()<<endl;
1877     //cout<<"\n\n";
1878
1879
1880   return;
1881 }
1882
1883   //__________________________________________________________________________
1884 Bool_t AliAnalysisTaskEMCALPhotonIsolation::FillGeneralHistograms(AliVCluster *coi, TLorentzVector vecCOI, Int_t index){
1885     //AliInfo("Inside FillGeneralHistograms\n");
1886
1887     // Fill the histograms for underlying events and isolation studies
1888
1889 // I would like to remove this part and fill the tracks multiplicity histogram in FillQAHistograms, is that ok for thnSparses? (especially cause here the histogram is filled several times per event)
1890   AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
1891   AliEmcalParticle *emcTrack = static_cast<AliEmcalParticle*>(tracks->GetNextAcceptParticle(0));
1892
1893   int nTracks=0;
1894   tracks->ResetCurrentID();
1895   while (emcTrack) {
1896     AliVTrack *track = emcTrack->GetTrack();
1897     if(!track) continue;
1898       // if(!(track->TestFilterBit("kHybrid"))) continue;
1899     nTracks++;
1900     emcTrack = static_cast<AliEmcalParticle*>(tracks->GetNextAcceptParticle());
1901   }
1902   fTrackMult->Fill(nTracks);
1903
1904
1905   Double_t eTCOI = 0., m02COI = 0.;
1906     //Int_t Ntracks;
1907     //Definition of the Array for Davide's Output
1908   const Int_t ndims =   fNDimensions;
1909   Double_t outputValues[ndims];
1910
1911   eTCOI = vecCOI.Et();
1912   m02COI = coi->GetM02();
1913
1914     //AliInfo(Form("M02 value: %lf\n",m02COI));
1915
1916     // ******** Isolation and UE calculation with different methods *********
1917
1918   Double_t eTThreshold = 5;
1919
1920   switch(fEtIsoMethod)
1921   {
1922     case 0:  // SumEt<EtThr
1923       eTThreshold = fEtIsoThreshold;
1924       break;
1925
1926     case 1:  // SumEt<%Ephoton
1927       eTThreshold = fEtIsoThreshold * eTCOI;
1928       break;
1929
1930     case 2: // Etmax<eTThreshold
1931       eTThreshold = fEtIsoThreshold;
1932       if( eTCOI<eTThreshold) // photon candidate, cuts have to be decided after studies
1933       {
1934         fEtIsolatedClust->Fill(eTCOI);
1935       }
1936       break;
1937   }
1938
1939     //DO NOT CHANGE EVER AGAIN THE FOLLOWING DEFINITIONS
1940   Float_t isoConeArea = TMath::Pi()*fIsoConeRadius*fIsoConeRadius;
1941   Float_t etaBandArea = (1.4 - 2.*fIsoConeRadius)*2.*fIsoConeRadius;
1942   Float_t phiBandArea = ((5./9.)*TMath::Pi()- 2.*fIsoConeRadius)*2.*fIsoConeRadius;
1943   Float_t etaBandAreaTr = (1.8 - 2.*fIsoConeRadius)*2.*fIsoConeRadius;
1944   Float_t phiBandAreaTr = (2.*TMath::Pi() - 4.*fIsoConeRadius)*2.*fIsoConeRadius;//there is a 2 more because of the JET CONE B2B
1945   Float_t perpConesArea = 2.*isoConeArea;
1946   Float_t fullTPCArea = 1.8*2.*TMath::Pi()-fIsoConeRadius*(TMath::Pi()*fIsoConeRadius + 3.6);
1947
1948   Float_t isolation, ue;
1949
1950   if(!fTPC4Iso){ //EMCAL Only for Acceptance of Clusters
1951     switch(fIsoMethod)
1952     {
1953       case 0: //EMCAL CELLS
1954
1955         switch (fUEMethod)
1956       {
1957         case 0: //phi band
1958           EtIsoCellPhiBand(vecCOI, isolation, ue);
1959             //Normalisation ue wrt Area - TO DO-
1960           ue = ue * (isoConeArea / phiBandArea);
1961           fPhiBandUECells->Fill(vecCOI.Pt() , ue);
1962           fEtIsoCells->Fill(isolation);
1963           if(isolation<eTThreshold)
1964           {
1965             fEtIsolatedCells->Fill(eTCOI);
1966             fEtisolatedT=eTCOI;
1967             fPtisolatedT=vecCOI.Pt();
1968           }
1969           break;
1970         case 1: //eta band
1971           EtIsoCellEtaBand(vecCOI, isolation, ue);
1972             //Normalisation ue wrt Area - TO DO-
1973           ue = ue * (isoConeArea / etaBandArea);
1974           fEtaBandUECells->Fill(vecCOI.Pt() , ue);
1975           fEtIsoCells->Fill(isolation);
1976           if(isolation<eTThreshold)
1977           {
1978             fEtIsolatedCells->Fill(eTCOI);
1979             fEtisolatedT=eTCOI;
1980             fPtisolatedT=vecCOI.Pt();
1981           }
1982           break;
1983       }
1984         break;
1985
1986       case 1: //EMCAL CLUSTERS
1987
1988         switch (fUEMethod)
1989       {
1990         case 0: //phi band
1991           EtIsoClusPhiBand(vecCOI, isolation, ue,index);
1992             //Normalisation ue wrt Area - TO DO-
1993           ue = ue * (isoConeArea / phiBandArea);
1994           fPhiBandUEClust->Fill(vecCOI.Pt() , ue);
1995           break;
1996         case 1: //eta band
1997           EtIsoClusEtaBand(vecCOI, isolation, ue,index);
1998             //Normalisation ue wrt Area - TO DO-
1999           ue = ue * (isoConeArea / etaBandArea);
2000           fEtaBandUEClust->Fill(vecCOI.Pt() , ue);
2001           break;
2002
2003           fEtIsoClust->Fill(vecCOI.Pt(),isolation);
2004           if(isolation<eTThreshold)
2005           {
2006            fPtvsM02iso->Fill(vecCOI.Pt(),coi->GetM02());
2007            fPtIsolatedNClust->Fill(vecCOI.Pt());
2008            fPtisoT=vecCOI.Pt();
2009            fM02isoT=m02COI;
2010
2011             if(fM02mincut < m02COI && m02COI < fM02maxcut)
2012             {
2013             fEtIsolatedClust->Fill(eTCOI);
2014             fEtisolatedT=eTCOI;
2015             fPtisolatedT=vecCOI.Pt();
2016             }
2017           }
2018
2019           else
2020           {
2021               fPtvsM02noiso->Fill(vecCOI.Pt(),coi->GetM02());
2022               fPtnoisoT=vecCOI.Pt();
2023               fM02noisoT=m02COI;
2024           }
2025       }
2026       break;
2027
2028       case 2: //TRACKS (TBD which tracks) //EMCAL tracks
2029         switch (fUEMethod)
2030       {
2031         case 0: //phi band
2032           PtIsoTrackPhiBand(vecCOI, isolation, ue);
2033             //Normalisation ue wrt Area - TO DO-
2034           ue = ue * (isoConeArea / phiBandAreaTr);
2035           fPhiBandUETracks->Fill(vecCOI.Pt() , ue);
2036         case 1: //eta band
2037           PtIsoTrackEtaBand(vecCOI, isolation, ue);
2038             //Normalisation ue wrt Area - TO DO-
2039           ue = ue * (isoConeArea / etaBandAreaTr);
2040           fEtaBandUETracks->Fill(vecCOI.Pt() , ue);
2041           break;
2042             // case 2: //Cones
2043             // PtIsoTrackOrthCones(vecCOI, absId, isolation, ue);
2044             // break;
2045             // case 3: //Full TPC
2046             // PtIsoTrackFullTPC(vecCOI, absId, isolation, ue);
2047             // break;
2048       }
2049       break;
2050
2051     // Fill histograms for isolation
2052           fPtIsoTrack->Fill(vecCOI.Pt() , isolation);
2053           if(isolation<eTThreshold)
2054           {
2055               fPtvsM02iso->Fill(vecCOI.Pt(),coi->GetM02());
2056               fPtIsolatedNTracks->Fill(vecCOI.Pt());
2057               fPtisoT=vecCOI.Pt();
2058               fM02isoT=m02COI;
2059
2060             if(fM02mincut < m02COI && m02COI < fM02maxcut)
2061             {
2062             fEtIsolatedTracks->Fill(eTCOI);
2063             fEtisolatedT=eTCOI;
2064             fPtisolatedT=vecCOI.Pt();
2065             }
2066           }
2067           else
2068           {
2069               fPtvsM02noiso->Fill(vecCOI.Pt(),coi->GetM02());
2070               fPtnoisoT=vecCOI.Pt();
2071               fM02noisoT=m02COI;
2072           }
2073     }
2074   }
2075   else{  //EMCAL + TPC (Only tracks for the Isolation since IsoCone Goes Out of EMCAL)
2076     switch (fUEMethod)
2077     {
2078       case 0: //phi band
2079         PtIsoTrackPhiBand(vecCOI, isolation, ue);
2080           //Normalisation ue wrt Area - TO DO-
2081         ue = ue * (isoConeArea / phiBandAreaTr);
2082         fPhiBandUETracks->Fill(vecCOI.Pt() , ue);
2083         break;
2084       case 1: //eta band
2085         PtIsoTrackEtaBand(vecCOI, isolation, ue);
2086           //Normalisation ue wrt Area - TO DO-
2087         ue = ue * (isoConeArea / etaBandAreaTr);
2088         fEtaBandUETracks->Fill(vecCOI.Pt() , ue);
2089         break;
2090       case 2: //Cones
2091         PtIsoTrackOrthCones(vecCOI, isolation, ue);
2092           //Normalisation ue wrt Area - TO DO-
2093         ue = ue * (isoConeArea / perpConesArea);
2094         fPerpConesUETracks ->Fill(vecCOI.Pt() , ue);
2095         break;
2096       case 3: //Full TPC
2097         PtIsoTrackFullTPC(vecCOI, isolation, ue);
2098           //Normalisation ue wrt Area - TO DO-
2099         ue = ue * (isoConeArea / fullTPCArea);
2100         fTPCWithoutIsoConeB2BbandUE->Fill(vecCOI.Pt() , ue);
2101         break;
2102
2103 // fill histograms for isolation
2104         fPtIsoTrack->Fill(vecCOI.Pt() , isolation);
2105         if(isolation<eTThreshold)
2106         {
2107             fPtvsM02iso->Fill(vecCOI.Pt(),coi->GetM02());
2108             fPtIsolatedNTracks->Fill(vecCOI.Pt());
2109             fPtisoT=vecCOI.Pt();
2110             fM02isoT=m02COI;
2111
2112             if(fM02mincut < m02COI && m02COI < fM02maxcut)
2113             {
2114                 fEtIsolatedTracks->Fill(eTCOI);
2115                 fEtisolatedT=eTCOI;
2116                 fPtisolatedT=vecCOI.Pt();
2117             }
2118         }
2119         else
2120         {
2121             fPtvsM02noiso->Fill(vecCOI.Pt(),coi->GetM02());
2122             fPtnoisoT=vecCOI.Pt();
2123             fM02noisoT=m02COI;
2124         }
2125     }
2126   }
2127
2128
2129   /*  Here we should call something to know the number of tracks...
2130    Soon I'll put in this version the "old way"; please let me know if
2131    any of you could do the same with the JET framework*/
2132
2133   switch(fWho) {
2134     case 0:
2135       flambda0T=m02COI; // for all neutral clusters
2136       fEtT=vecCOI.Et(); // for all neutral clusters
2137       fPtT=vecCOI.Pt(); // for all neutral clusters
2138       fetaT=vecCOI.Eta(); // for all neutral clusters
2139       fphiT=vecCOI.Phi(); //for all neutral clusters
2140       fsumEtisoconeT=isolation;
2141         //         AliError(Form("lambda 0 %f",flambda0T));
2142       fsumEtUE=ue;
2143
2144       fOutputTree->Fill();
2145       break;
2146
2147     case 1:
2148       outputValues[0] = nTracks;
2149       outputValues[1] = eTCOI;
2150       outputValues[2] = m02COI;
2151       outputValues[3] = isolation;
2152       outputValues[4] = ue;
2153       outputValues[5] = isolation-ue;
2154       outputValues[6] = vecCOI.Eta();
2155       outputValues[7] = vecCOI.Phi();
2156       /*if (fIsMC) {
2157        outputValues[8] = ptmc;
2158        outputValues[9] = mcptsum;
2159        }*/
2160       fOutputTHnS -> Fill(outputValues);
2161       break;
2162         //                 //            fOutPTMAX -> Fill(outputValues[1],outputValues[2],);
2163   }
2164   return kTRUE;
2165 }
2166
2167
2168   //_________________________________________________________________________
2169
2170 void AliAnalysisTaskEMCALPhotonIsolation::AnalyzeMC(){
2171
2172   if (!fIsMC)
2173     return;
2174     //AliInfo(Form("It's a MC analysis %e",fAODMCParticles));
2175   if(!fStack && !fAODMCParticles)
2176   {cout<<"no stack saved\n"; return;}
2177     //AliInfo("there's a List of particles");
2178     //DO THIS ALSO FOR ESDs
2179
2180   Double_t eT, sumEiso, sumUE,phi, eta, distance, phip, etap, mcfirstEnergy;
2181
2182   if(fAODMCParticles->GetEntries() < 1){
2183     AliError("number of tracks insufficient");
2184     return;
2185   }
2186   int nDimMC = fMCDimensions;
2187   Double_t outputValuesMC[nDimMC];
2188
2189   Int_t nTracks = fAODMCParticles->GetEntriesFast();
2190   Int_t nFSParticles = 0;
2191   AliAODMCParticle *multTracks;
2192
2193   for(int a=0; a<nTracks; a++){
2194
2195     multTracks = static_cast<AliAODMCParticle*>(fAODMCParticles->At(a));
2196
2197     if(multTracks->IsPrimary() && multTracks->IsPhysicalPrimary() && multTracks->GetStatus()<10){
2198       if(TMath::Abs(multTracks->Eta())<=0.9 && multTracks->Charge()!= 0)
2199         nFSParticles++;
2200       else
2201         continue;
2202     }//implement final state particle condition
2203     else
2204       continue;
2205   }
2206     //AliInfo(Form("number of particles in the array %d",nTracks));
2207   AliAODMCParticle *mcpart, *mom, *mcpp,*mcsearch, *mcfirst, *mcfirstmom,*matchingtrack;
2208
2209   Bool_t prompt=kFALSE;
2210   Double_t mcEnergy, maxE, energy;
2211   Int_t pdg, mompdg, photonlabel;
2212   Double_t mcFirstEta=0., mcFirstPhi=0.;
2213
2214   Float_t isoConeArea = TMath::Pi()*fIsoConeRadius*fIsoConeRadius;
2215   Float_t etaBandArea = (1.4 - 2.*fIsoConeRadius)*2.*fIsoConeRadius;
2216   Float_t phiBandArea = ((5./9.)*TMath::Pi()- 2.*fIsoConeRadius)*2.*fIsoConeRadius;
2217   Float_t etaBandAreaTr = (1.8 - 2.*fIsoConeRadius)*2.*fIsoConeRadius;
2218   Float_t phiBandAreaTr = (2.*TMath::Pi() - 4.*fIsoConeRadius)*2.*fIsoConeRadius;//there is a 2 more because of the JET CONE B2B
2219   Float_t perpConesArea = 2.*isoConeArea;
2220   Float_t fullTPCArea = 1.8*2.*TMath::Pi()-fIsoConeRadius*(TMath::Pi()*fIsoConeRadius + 3.6);
2221
2222     // AliAODMCParticle *mcfirst = static_cast<AliAODMCParticle*>(fAODMCParticles->At(0));
2223     //AliAODMCParticle *mcp, *mcpmaxE, *mcpp, *mom;
2224   if(!fisLCAnalysis){
2225       //Loop on the event
2226     for(int iTr=0;iTr<nTracks;iTr++){
2227
2228       mcEnergy=0.;energy =0;
2229       eT=0.; phi=0.; eta=0.;
2230
2231       mcpart = static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTr));
2232       /*if(mcpart->GetStatus()<10 /*&& mcpart->IsPrimary()==0){
2233        //if(mcpart->GetMCProcessCode()!=0){
2234        if(mcpart->IsPrimary() && mcpart->IsPhysicalPrimary()){
2235
2236        if(fcount==388){
2237        cout<<"Particle in Stack: "<<iTr<<"\tLabel : "<<mcpart->Label();
2238        mcpart->Print();
2239        cout<<"IsPrimary: "<<mcpart->IsPrimary()<<"\t IsSecondaryfromDecay: "<<mcpart->IsSecondaryFromWeakDecay()<<"\t Is SecondaryfromMaterial: "<<mcpart->IsSecondaryFromMaterial()<<endl;
2240        int momidx = mcpart->GetMother();
2241        mom = static_cast<AliAODMCParticle*>(fAODMCParticles->At(momidx));
2242        cout<<"NOW THE MOTHER\n"<<momidx<<"\t";
2243        mom->Print();
2244        cout<<"IsPrimary: "<<mom->IsPrimary()<<"\t IsSecondaryfromDecay: "<<mom->IsSecondaryFromWeakDecay()<<"\t Is SecondaryfromMaterial: "<<mom->IsSecondaryFromMaterial()<<endl;
2245
2246        cout<<"\n\n"<<endl;
2247        }
2248        }
2249        }*/
2250       if(mcpart->GetStatus()>10) continue;
2251       if(!mcpart->IsPrimary()) continue;
2252       if(!mcpart->IsPhysicalPrimary()) continue;
2253
2254       pdg = mcpart->GetPdgCode();
2255
2256       if(pdg!=22)
2257         continue;
2258
2259       eta = mcpart->Eta();
2260       phi = mcpart->Phi();
2261
2262         //check photons in EMCAL //to be redefined with fIsoConeR
2263       if((TMath::Abs(eta)>0.3) || (phi<1.8 || phi>(TMath::Pi()-0.4)))
2264         continue;
2265
2266       photonlabel = iTr;
2267       int momidx = mcpart->GetMother();
2268
2269       mom = static_cast<AliAODMCParticle*>(fAODMCParticles->At(momidx));
2270       mompdg= TMath::Abs(mom->GetPdgCode());
2271
2272       eT= mcpart->E()*TMath::Sin(mcpart->Theta()); //transform to transverse Energy
2273
2274       fphietaPhotons->Fill(eta,phi,eT);
2275
2276
2277       bool foundmatch=kFALSE;
2278       for(int m=0;m<nTracks && foundmatch==kFALSE;m++){
2279         if(m==iTr) continue;
2280
2281         matchingtrack = static_cast<AliAODMCParticle*>(fAODMCParticles->At(m));
2282
2283         if(! matchingtrack->IsPrimary()) continue;
2284         if(! matchingtrack->IsPhysicalPrimary()) continue;
2285         if(matchingtrack->GetStatus()> 10 ) continue;
2286
2287         Double_t etamatching = matchingtrack->Eta();
2288         Double_t phimatching = matchingtrack->Phi();
2289
2290         if(TMath::Abs(eta-etamatching)<=0.02 && TMath::Abs(phi-phimatching)<=0.03){
2291           foundmatch=kTRUE;
2292           fphietaOthers->Fill(matchingtrack->Eta(),matchingtrack->Phi(),eT);
2293           fphietaOthersBis->Fill(matchingtrack->Eta(),matchingtrack->Phi(),matchingtrack->Pt());
2294         }
2295       }
2296
2297         //int grandmaidx = mom->GetMother();
2298
2299       /*if((mcpart->IsPrimary()) || (mompdg==22 && grandmaidx== -1)){
2300        prompt = kTRUE;
2301        }
2302        else{
2303        prompt = kFALSE;
2304        }
2305        cout<<iTr<<"\t";
2306        mcpart->Print();
2307        cout<<"IsPrimary: "<<mcpart->IsPrimary()<<"\t IsSecondaryfromDecay: "<<mcpart->IsSecondaryFromWeakDecay()<<"\t Is SecondaryfromMaterial: "<<mcpart->IsSecondaryFromMaterial()<<endl;
2308        cout<<"NOW THE MOTHER\n"<<momidx<<"\t";
2309        mom->Print();
2310        cout<<"IsPrimary: "<<mom->IsPrimary()<<"\t IsSecondaryfromDecay: "<<mom->IsSecondaryFromWeakDecay()<<"\t Is SecondaryfromMaterial: "<<mom->IsSecondaryFromMaterial()<<endl;
2311        cout<<"Coordinates of the photon (eta,phi)= ("<<eta<<","<<phi<<")"<<endl;
2312        cout<<"\n\n"<<endl;
2313        */
2314
2315       distance=0.;
2316       phip=0., etap=0.;
2317       sumEiso=0,sumUE=0;
2318
2319       for(int iTrack=0;iTrack<nTracks;iTrack++){
2320
2321         if(iTrack==photonlabel)continue;
2322
2323         mcpp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
2324
2325         if(!mcpp) continue;
2326
2327         if((mcpp->GetPdgCode())==22) continue;
2328
2329         if(mcpp->GetStatus()>10) continue;
2330
2331         phip = mcpp->Phi();
2332         etap = mcpp->Eta();
2333
2334           //Depending on which Isolation method and UE method is considered.
2335
2336         distance= TMath::Sqrt((phi-phip)*(phi-phip) + (eta-etap)*(eta-etap));
2337           //cout<<iTrack<<endl;
2338           //cout<<"Coordinates of this particle (eta,phi)= ("<<etap<<","<<phip<<")"<<endl;
2339           //cout<<"distance of this particle from the photon: "<<distance<<endl;
2340
2341         if(distance <= 0.4){ //to be changed with fIsoConeR
2342           sumEiso += mcpp->E()*TMath::Sin(mcpp->Theta());
2343             //cout<<"\n\n Transverse Energy of this particle : "<<mcpp->E()*TMath::Sin(mcpp->Theta())<<endl;
2344         }
2345         else{
2346           if(!fTPC4Iso){
2347             if(TMath::Abs(etap)>=0.7 || (phip<=1.4 || phip>= TMath::Pi()))
2348               continue;
2349             else{
2350               switch(fUEMethod){
2351                 case 0: //Phi band
2352                   if(TMath::Abs(eta-etap)<0.4) //to be changed with fIsoConeRadius
2353                     sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2354                   else continue;
2355
2356                   break;
2357                 case 1: //Eta band
2358                   if(TMath::Abs(phi-phip)<0.4)
2359                     sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2360                   else continue;
2361
2362                   break;
2363               }
2364             }
2365           }
2366           else{
2367             if(TMath::Abs(etap)>=1.0)
2368               continue;
2369             else{
2370               switch(fUEMethod){
2371                 case 0: //Phi band
2372                 {if(TMath::Abs(eta-etap)<0.4) //to be changed with fIsoConeRadius
2373                   sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2374                 else continue;
2375                   break;
2376                 }
2377                 case 1: //Eta band
2378                 {  if(TMath::Abs(phi-phip)<0.4)
2379                   sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2380                 else continue;
2381
2382                   break;
2383                 }
2384                 case 2: //Orthogonal Cones
2385                 { double etacone1= eta;
2386                   double etacone2= eta;
2387                   double phicone1= phi - TMath::PiOver2();
2388                   double phicone2= phi + TMath::PiOver2();
2389
2390                   if (phicone1 < 0.) phicone1 += 2*TMath::Pi();
2391
2392                   if(TMath::Sqrt(TMath::Power(etap-etacone1,2)+TMath::Power(phip-phicone1,2))< 0.4 ||
2393                      TMath::Sqrt(TMath::Power(etap-etacone2,2)+TMath::Power(phip-phicone2,2))< 0.4) //to be changed with fIsoConeRadius
2394                   {sumUE += mcpp->Pt();}
2395                   else continue;
2396
2397                   break;
2398                 }
2399                 case 3: //Full TPC
2400                 {    //                  Double_t phiup= phi +TMath::Pi()+fIsoConeRadius;
2401                     //                  Double_t phidown= phi +TMath::Pi()-fIsoConeRadius;
2402                     //
2403                     //                  if(phip < phidown || phip > phiup ) //TO BE CHECKED
2404                     //                    continue;
2405                   break;
2406                 }
2407               }
2408             }
2409           }
2410         }
2411       }
2412         //cout<<"\n\nTotal Energy inside the Isolation Cone : "<<sumEiso<<endl;
2413         //cout<<"Total UE Energy : "<<sumUE<<endl;
2414       if(!fTPC4Iso){
2415         switch (fUEMethod){
2416           case 0:
2417             sumUE = sumUE * (isoConeArea / phiBandArea);
2418             break;
2419           case 1:
2420             sumUE = sumUE * (isoConeArea / etaBandArea);
2421             break;
2422         }
2423       }
2424       else{
2425         switch (fUEMethod){
2426           case 0:
2427             sumUE = sumUE * (isoConeArea / phiBandAreaTr);
2428             break;
2429           case 1:
2430             sumUE = sumUE * (isoConeArea / etaBandAreaTr);
2431             break;
2432           case 2:
2433             sumUE = sumUE * (isoConeArea / perpConesArea);
2434             break;
2435           case 3:
2436             sumUE = sumUE * (isoConeArea / fullTPCArea);
2437             break;
2438         }
2439       }
2440         // cout<<"Total SCALED UE Energy : "<<sumUE<<" calculated with method "<<fUEMethod<<"which brings a normalisation factor: "<<phiBandArea<<endl;
2441
2442       outputValuesMC[0] = nFSParticles;
2443       outputValuesMC[1] = eT;
2444         //outputValuesMC[2] = sumEiso;
2445         //outputValuesMC[3] = sumUE;
2446       outputValuesMC[2] = mompdg;
2447         //outputValuesMC[5] = eta;
2448         //outputValuesMC[6] = phi;
2449       outputValuesMC[3] = mcpart->GetLabel();
2450         // EtaPhiMCPhoton
2451         // EtMC
2452         // EtIsoCone
2453         // EtMother
2454         // UE Et
2455         // Mother PDG
2456         //fill some histograms or a THnSparse or a TTree.
2457       fOutMCTruth -> Fill(outputValuesMC);
2458
2459
2460     }
2461   }
2462   else{
2463     maxE=0.;
2464     int indexmaxE=0;
2465       //getting the index of the particle with the maximum energy.
2466     for(int iTr=0;iTr<nTracks;iTr++){
2467       mcsearch= static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTr));
2468
2469       if(!mcsearch) continue;
2470
2471       if(mcsearch->GetStatus()>10) continue;
2472       if(mcsearch->GetPdgCode()!=22) continue;
2473       if(TMath::Abs(mcsearch->Eta())>0.3) continue;
2474       if(mcsearch->Phi()<= 1.8 ||mcsearch->Phi()>= TMath::Pi()) continue;
2475
2476       mcfirstEnergy= mcsearch->E();
2477       if(mcfirstEnergy>maxE){
2478         maxE=mcfirstEnergy;
2479         indexmaxE=iTr;
2480       }
2481       else continue;
2482     }
2483     mcfirst= static_cast<AliAODMCParticle*>(fAODMCParticles->At(indexmaxE));
2484     mcfirstEnergy=mcfirst->E()*TMath::Sin(mcfirst->Theta());
2485
2486     int momidx= mcfirst->GetMother();
2487     mcfirstmom =  static_cast<AliAODMCParticle*>(fAODMCParticles->At(momidx));
2488     mompdg= TMath::Abs(mcfirstmom->GetPdgCode());
2489     mcFirstEta = mcfirst->Eta();
2490     mcFirstPhi = mcfirst->Phi();
2491
2492     phip=0., etap=0.;
2493     sumEiso=0,sumUE=0;
2494
2495     for(Int_t iTrack=1;iTrack<nTracks ;iTrack++){
2496       if(iTrack==indexmaxE) continue;
2497
2498       mcpp= static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
2499       phip = mcpp->Phi();
2500       etap = mcpp->Eta();
2501       if(!mcpp)
2502         continue;
2503
2504       if(mcpp->GetStatus()>10) continue;
2505       if(mcpp->GetPdgCode()==22)continue;
2506       if(TMath::Abs(etap>0.7)) continue;
2507       if(phip<=1.4 || phip>= TMath::Pi()) continue;
2508       distance=0.;
2509       distance= TMath::Sqrt((mcFirstPhi- phip)*(mcFirstPhi- phip) + (mcFirstEta- etap)*(mcFirstEta- etap));
2510
2511       if(distance<=0.4){
2512         sumEiso += mcpp->E()*TMath::Sin(mcpp->Theta());
2513       }
2514       else{
2515         if(!fTPC4Iso){
2516           if(TMath::Abs(etap)>=0.7 || (phip<=1.4 || phip>= TMath::Pi()))
2517             continue;
2518           else{
2519             switch(fUEMethod){
2520               case 0: //Phi band
2521                 if(TMath::Abs(mcFirstEta-etap)<0.4) //to be changed with fIsoConeRadius
2522                   sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2523                 else continue;
2524
2525                 break;
2526               case 1: //Eta band
2527                 if(TMath::Abs(mcFirstPhi-phip)<0.4)
2528                   sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2529                 else continue;
2530
2531                 break;
2532             }
2533           }
2534         }
2535         else{
2536           if(TMath::Abs(etap)>=1.0)
2537             continue;
2538           else{
2539             switch(fUEMethod){
2540               case 0: //Phi band
2541               {  if(TMath::Abs(mcFirstEta-etap)<0.4) //to be changed with fIsoConeRadius
2542                 sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2543               else continue;
2544                 break;
2545               }
2546               case 1: //Eta band
2547               {if(TMath::Abs(mcFirstPhi-phip)<0.4)
2548                 sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2549               else continue;
2550
2551                 break;
2552               }
2553               case 2: //Orthogonal Cones
2554               {    double phicone1= mcFirstPhi - TMath::PiOver2();
2555                 double phicone2= mcFirstPhi + TMath::PiOver2();
2556
2557                 if (phicone1 < 0.) phicone1 += 2*TMath::Pi();
2558
2559                 if(TMath::Sqrt(TMath::Power(etap-mcFirstEta,2)+TMath::Power(phip-phicone1,2))< 0.4 ||
2560                    TMath::Sqrt(TMath::Power(etap-mcFirstEta,2)+TMath::Power(phip-phicone2,2))< 0.4) //to be changed with fIsoConeRadius
2561                 {sumUE += mcpp->Pt();}
2562                 else continue;
2563                 break;
2564               }
2565               case 3: //Full TPC
2566               {    //                  Double_t phiup= phi +TMath::Pi()+fIsoConeRadius;
2567                   //                  Double_t phidown= phi +TMath::Pi()-fIsoConeRadius;
2568                   //
2569                   //                  if(phip < phidown || phip > phiup ) //TO BE CHECKED
2570                   //                    continue;
2571                 break;
2572               }
2573             }
2574           }
2575         }
2576       }
2577     }
2578       //  cout<<"\n\nTotal Energy inside the Isolation Cone : "<<sumEiso<<endl;
2579     if(!fTPC4Iso){
2580       switch (fUEMethod){
2581         case 0:
2582           sumUE = sumUE * (isoConeArea / phiBandArea);
2583           break;
2584         case 1:
2585           sumUE = sumUE * (isoConeArea / etaBandArea);
2586           break;
2587       }
2588     }
2589     else{
2590       switch (fUEMethod){
2591         case 0:
2592           sumUE = sumUE * (isoConeArea / phiBandAreaTr);
2593           break;
2594         case 1:
2595           sumUE = sumUE * (isoConeArea / etaBandAreaTr);
2596           break;
2597         case 2:
2598           sumUE = sumUE * (isoConeArea / perpConesArea);
2599           break;
2600         case 3:
2601           sumUE = sumUE * (isoConeArea / fullTPCArea);
2602           break;
2603       }
2604     }
2605       //cout<<"Total UE Energy : "<<sumUE<<" calculated with method "<<fUEMethod<<endl;
2606
2607       //Fill the Output TTree for MC Truth
2608   }
2609
2610   return;
2611 }
2612
2613
2614 /*
2615  else{
2616
2617  eT = mcpmaxE->E(); //transform to transverse Energy
2618  phi = mcpmaxE->Phi();
2619  eta = mcpmaxE->Eta();
2620  distance=0.;
2621  phip=0., etap=0.;
2622  for(iTrack=0;iTrack<nTracks;iTrack++){
2623
2624  if(iTrack==maxindex)
2625  continue;
2626
2627  mcpp=static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
2628  if(!mcpp)
2629  continue;
2630
2631  phip = mcpp->Phi();
2632  etap = mcpp->Eta();
2633  distance= TMath::Sqrt((phi-phip)*(phi-phip)+(eta-etap)*(eta-etap));
2634
2635  if(distance <= 0.4) //to be changed with fIsoConeR
2636  sum += mcpp->Pt();
2637
2638  else
2639  continue;
2640  }
2641  //fill some histograms (PDG, ET, Eta/phi distributions, sum in pT)
2642  }
2643  */
2644