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 <THnSparse.h>
  17 #include "AliAnalysisManager.h"
  18 #include "AliCentrality.h"
  19 #include "AliEMCALGeometry.h"
  20 #include "AliESDEvent.h"
  21 #include "AliAODEvent.h"
  22 #include "AliLog.h"
  23 #include "AliVCluster.h"
  24 #include "AliVEventHandler.h"
  25 #include "AliVParticle.h"
  26 #include "AliClusterContainer.h"
  27 #include "AliVTrack.h"
  28 #include "AliEmcalParticle.h"
  29 #include "AliParticleContainer.h"
  30 #include "AliAODCaloCluster.h"
  31 #include "AliESDCaloCluster.h"
  32 #include "AliVCaloCells.h"
  33 #include "AliPicoTrack.h"
  34
  35 #include "AliAnalysisTaskEMCALPhotonIsolation.h"
  36
  37 using std::cout;
  38 using std::endl;
  39
  40 ClassImp(AliAnalysisTaskEMCALPhotonIsolation)
  41
  42 //________________________________________________________________________
  43 AliAnalysisTaskEMCALPhotonIsolation::AliAnalysisTaskEMCALPhotonIsolation() :
  44 AliAnalysisTaskEmcal("AliAnalysisTaskEMCALPhotonIsolation",kTRUE),
  45 //fParticleCollArray(),
  46 fNCluster(0),
  47 fWho(-1),
  48 //fOutputList(0),
  49 fTrackMult(0),
  50 fTrackMultEMCAL(0),
  51 fClustMult(0),
  52 fPVZBefore(0),
  53 fEtaPhiCell(0),
  54 fEtaPhiClus(0),
  55 fClusEvsClusT(0),
  56 fGoodEventsOnPVZ(0),
  57 fPT(0),
  58 fM02(0),
  59 fNLM(0),
  60 fDeltaETAClusTrackVSpT(0),
  61 fDeltaPHIClusTrackVSpT(0),
  62 fEtIsoCells(0),
  63 fEtIsoClust(0),
  64 fPtIsoTrack(0),
  65 fPtEtIsoTC(0),
  66 fPhiBandUEClust(0),
  67 fEtaBandUEClust(0),
  68 fPhiBandUECells(0),
  69 fEtaBandUECells(0),
  70 fPhiBandUETracks(0),
  71 fEtaBandUETracks(0),
  72 fPerpConesUETracks(0),
  73 fTPCWithoutIsoConeB2BbandUE(0),
  74 fNTotClus10GeV(0),
  75 fRecoPV(0),
  76 fEtIsolatedCells(0),
  77 fEtIsolatedClust(0),
  78 fEtIsolatedTracks(0),
  79 fTest(0),
  80 fOutputTHnS(0),
  81 fOutPTMAX(0),
  82 fOutputTree(0),
  83 fIsoConeRadius(0.4),
  84 fEtIsoMethod(0),
  85 fEtIsoThreshold(4),
  86 fdetacut(0.025),
  87 fdphicut(0.03),
  88 fM02mincut(0.1),
  89 fM02maxcut(0.3),
  90 fQA(0),
  91 fIsMC(0),
  92 fTPC4Iso(0),
  93 fIsoMethod(0),
  94 fUEMethod(0),
  95 fVertex(0),
  96 fNDimensions(0),
  97 fisLCAnalysis(0),
  98 fevents(0),
  99 flambda0T(0),
 100 fEtT(0),
 101 fPtT(0),
 102 fEtisoT(0),
 103 fPtisoT(0),
 104 fetaT(0),
 105 fphiT(0),
 106 fsumEtisoconeT(0),
 107 fsumEtUE(0)
 108 //tracks(0),
 109 //clusters(0)
 110
 111 {
 112   // Default constructor.
 113
 114   //fParticleCollArray.SetOwner(kTRUE);
 115   // for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
 116
 117   SetMakeGeneralHistograms(kTRUE);
 118 }
 119
 120 //________________________________________________________________________
 121 AliAnalysisTaskEMCALPhotonIsolation::AliAnalysisTaskEMCALPhotonIsolation(const char *name, Bool_t histo) :
 122 AliAnalysisTaskEmcal(name, histo),
 123 //fParticleCollArray(),
 124 fNCluster(0),
 125 fWho(-1),
 126 //fOutputList(0),
 127 fTrackMult(0),
 128 fTrackMultEMCAL(0),
 129 fClustMult(0),
 130 fPVZBefore(0),
 131 fEtaPhiCell(0),
 132 fEtaPhiClus(0),
 133 fClusEvsClusT(0),
 134 fGoodEventsOnPVZ(0),
 135 fPT(0),
 136 fM02(0),
 137 fNLM(0),
 138 fDeltaETAClusTrackVSpT(0),
 139 fDeltaPHIClusTrackVSpT(0),
 140 fEtIsoCells(0),
 141 fEtIsoClust(0),
 142 fPtIsoTrack(0),
 143 fPtEtIsoTC(0),
 144 fPhiBandUEClust(0),
 145 fEtaBandUEClust(0),
 146 fPhiBandUECells(0),
 147 fEtaBandUECells(0),
 148 fPhiBandUETracks(0),
 149 fEtaBandUETracks(0),
 150 fPerpConesUETracks(0),
 151 fTPCWithoutIsoConeB2BbandUE(0),
 152 fNTotClus10GeV(0),
 153 fRecoPV(0),
 154 fEtIsolatedCells(0),
 155 fEtIsolatedClust(0),
 156 fEtIsolatedTracks(0),
 157 fTest(0),
 158 fOutputTHnS(0),
 159 fOutPTMAX(0),
 160 fOutputTree(0),
 161 fIsoConeRadius(0.4),
 162 fEtIsoMethod(0),
 163 fEtIsoThreshold(4),
 164 fdetacut(0.025),
 165 fdphicut(0.03),
 166 fM02mincut(0.1),
 167 fM02maxcut(0.3),
 168 fQA(0),
 169 fIsMC(0),
 170 fTPC4Iso(0),
 171 fIsoMethod(0),
 172 fUEMethod(0),
 173 fVertex(0),
 174 fNDimensions(0),
 175 fisLCAnalysis(0),
 176 fevents(0),
 177 flambda0T(0),
 178 fEtT(0),
 179 fPtT(0),
 180 fEtisoT(0),
 181 fPtisoT(0),
 182 fetaT(0),
 183 fphiT(0),
 184 fsumEtisoconeT(0),
 185 fsumEtUE(0)
 186 //tracks(0),
 187 //clusters(0)
 188
 189 {
 190   // Standard constructor.
 191
 192   //fParticleCollArray.SetOwner(kTRUE);
 193   //    for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
 194
 195   SetMakeGeneralHistograms(kTRUE);
 196 }
 197
 198 //________________________________________________________________________
 199 AliAnalysisTaskEMCALPhotonIsolation::~AliAnalysisTaskEMCALPhotonIsolation(){
 200   // Destructor
 201 }
 202
 203
 204 //________________________________________________________________________
 205 void AliAnalysisTaskEMCALPhotonIsolation::UserCreateOutputObjects(){
 206   // Create ouput histograms and THnSparse and TTree
 207
 208   AliAnalysisTaskEmcal::UserCreateOutputObjects();
 209
 210
 211   if ((fIsoMethod == 0 || fIsoMethod == 1) && fTPC4Iso){
 212     cout<<"Error: Iso_Methods with CELLS and CLUSTERS work only within EMCAL "<<endl;
 213     cout<<"Please Set Iso_Method and TPC4Iso Accordingly!!"<<endl;
 214     return;
 215   }
 216   if ((fIsoMethod == 0 || fIsoMethod == 1) && fUEMethod> 1) {
 217     cout<<"Error: UE_Methods with CELLS and CLUSTERS work only within EMCAL"<<endl;
 218     cout<<"Please Set Iso_Method and UE_Method Accordingly!!"<<endl;
 219     return;
 220   }
 221
 222
 223   TString sIsoMethod="\0",sUEMethod="\0",sBoundaries="\0";
 224
 225   if(fIsoMethod==0)
 226     sIsoMethod = "Cells";
 227   else if(fIsoMethod==1)
 228     sIsoMethod = "Clust";
 229   else if(fIsoMethod==2)
 230     sIsoMethod = "Tracks";
 231
 232   if(fUEMethod==0)
 233     sUEMethod = "PhiBand";
 234   else if(fUEMethod==1)
 235     sUEMethod = "EtaBand";
 236   else if(fUEMethod==2)
 237     sUEMethod = "PerpCones";
 238   else if(fUEMethod==3)
 239     sUEMethod = "FullTPC";
 240
 241   if(fTPC4Iso)
 242     sBoundaries = "TPC Acceptance";
 243   else
 244     sBoundaries = "EMCAL Acceptance";
 245
 246   if(fWho>1 || fWho==-1){
 247     cout<<"Error!!! OutputMode Can Only Be 0: TTree; 1: THnSparse"<<endl;
 248     return;
 249   }
 250   else{
 251     fOutput = new TList();
 252     fOutput->SetOwner();
 253     //Initialize the common Output histograms
 254     switch (fWho)
 255     {
 256       case 0:
 257         //Initialization by Lucile and Marco
 258         fOutputTree = new TTree("OutTree","OutTree");
 259         fOutputTree->Branch("fevents",&fevents);
 260         fOutputTree->Branch("flambda0T",&flambda0T);
 261         fOutputTree->Branch("fEtT",&fEtT);
 262         fOutputTree->Branch("fPtT",&fPtT);
 263         fOutputTree->Branch("fEtisoT",&fEtisoT);
 264         fOutputTree->Branch("fPtTiso",&fPtisoT);
 265         fOutputTree->Branch("fetaT",&fetaT);
 266         fOutputTree->Branch("fphiT",&fphiT);
 267         fOutputTree->Branch("fsumEtisoconeT",&fsumEtisoconeT);
 268         fOutputTree->Branch("fsumEtUE",&fsumEtUE);
 269
 270         fOutput->Add(fOutputTree);
 271
 272         break;
 273       case 1:
 274         //Initialization by Davide;
 275
 276         TString sTitle;
 277         Int_t binTrackMult=100, binPT=70, binM02=100, binETiso=100, binETUE=110, binETisoUE=110, binetacl=140,binphicl=100, binPTMC=70;
 278         Int_t bins[] = {binTrackMult, binPT, binM02, binETiso, binETUE, binETisoUE, binetacl, binphicl, binPTMC, binPTMC};
 279
 280         fNDimensions = sizeof(bins)/sizeof(Int_t);
 281         const Int_t ndims =   fNDimensions;
 282
 283         Double_t xmin[]= {  0.,  0., 0., -10., -10., -10.,-1.0, 1. ,-10.,-10.};
 284
 285         Double_t xmax[]= {1000., 70., 2., 100., 100., 100., 1.0, 3.5, 60., 60.};
 286
 287         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());
 288
 289         fOutputTHnS =  new THnSparseF("fHnOutput",sTitle.Data(), ndims, bins, xmin, xmax);
 290
 291         fOutput->Add(fOutputTHnS);
 292
 293         Int_t binsbis[] = {binPT, binM02, binETisoUE};
 294         Double_t xminbis[] = {  0., 0., -10.};
 295         Double_t xmaxbis[] = {100., 2., 100.};
 296
 297         fOutPTMAX = new THnSparseF ("fOutPTMAX","3D matrix E_{#gamma} VS M02 VS pT_{max}^{cone}; E_{T}^{#gamma} (GeV/c); M02; p_{T}^{Iso}(GeV/c)",3,binsbis,xminbis,xmaxbis);
 298
 299         fOutput->Add(fOutPTMAX);
 300         break;
 301     }
 302   }
 303
 304   if(fQA){
 305     //Include QA plots to the OutputList //DEFINE BETTER THE BINNING AND THE AXES LIMITS
 306     fTrackMult = new TH1D ("hTrackMult","Tracks multiplicity Distribution",250,0.,1000.);
 307     fTrackMult->Sumw2();
 308     fOutput->Add(fTrackMult);
 309
 310     fTrackMultEMCAL = new TH1D ("hTrackMultEMCAL","Tracks multiplicity Distribution inside EMCAL acceptance",250,0.,1000.);
 311     fTrackMultEMCAL->Sumw2();
 312     fOutput->Add(fTrackMultEMCAL);
 313
 314     fClustMult = new TH1D ("hClustMult","Clusters multiplicity Distribution",250,0.,1000.);
 315     fClustMult->Sumw2();
 316     fOutput->Add(fClustMult);
 317
 318     fRecoPV = new TH1D("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5);
 319     fRecoPV->Sumw2();
 320     fOutput->Add(fRecoPV);
 321
 322     fPVZBefore = new TH1D ("hPVZDistr","Z Distribution for the Reconstructed Vertex",200,0.,40.);
 323     fPVZBefore->Sumw2();
 324     fOutput->Add(fPVZBefore);
 325
 326     fEtaPhiCell = new TH2D ("hEtaPhiCellActivity","",250,0.,1000., 250, 0., 1000.);
 327     fEtaPhiCell->Sumw2();
 328     fOutput->Add(fEtaPhiCell);
 329
 330     fEtaPhiClus = new TH2D ("hEtaPhiClusActivity","",250,0.,1000., 250, 0., 1000.);
 331     fEtaPhiClus->Sumw2();
 332     fOutput->Add(fEtaPhiClus);
 333
 334     fClusEvsClusT = new TH2D ("hEnVSTime","",250,0.,1000., 250,0.,1000.);
 335     fClusEvsClusT->Sumw2();
 336     fOutput->Add(fClusEvsClusT);
 337
 338     fDeltaETAClusTrackVSpT = new TH2D("hTC_Dz","Track-Cluster Dz vs pT of Cluster",100,0.,100.,100,-0.05,0.05);
 339     fDeltaETAClusTrackVSpT->Sumw2();
 340     fOutput->Add(fDeltaETAClusTrackVSpT);
 341
 342     fDeltaPHIClusTrackVSpT = new TH2D("hTC_Dx","Track-Cluster Dx vs pT of Cluster",100,0.,100.,100,-0.05,0.05);
 343     fDeltaPHIClusTrackVSpT->Sumw2();
 344     fOutput->Add(fDeltaPHIClusTrackVSpT);
 345
 346   }
 347   //Initialize only the Common THistos for the Three different output
 348
 349   fGoodEventsOnPVZ = new TH1D ("hGOODwrtPVZ","Number of Selected Events wrt Cut on Primary Vertex Z (0=disregarded,1=selected)",2,0.,2.);
 350   fGoodEventsOnPVZ->Sumw2();
 351   fOutput->Add(fGoodEventsOnPVZ);
 352
 353   fPT = new TH1D("hPt_NC","P_{T} distribution for Neutral Clusters",100,0.,100.);
 354   fPT->Sumw2();
 355   fOutput->Add(fPT);
 356
 357   fM02 = new TH2D("hM02_NC","M02 distribution for Neutral Clusters vs E",100,0.,100.,500,0.,5.);
 358   fM02->Sumw2();
 359   fOutput->Add(fM02);
 360
 361   fNLM = new TH1D("hNLM_NC","NLM distribution for Neutral Clusters",200,0.,4.);
 362   fNLM->Sumw2();
 363   fOutput->Add(fNLM);
 364
 365   fEtIsoCells = new TH1D("hEtIsoCell_NC","E_{T}^{iso cone} in iso cone distribution for Neutral Clusters with EMCAL Cells",100,0.,100.);
 366   fEtIsoCells->SetXTitle("#Sigma E_{T}^{iso cone} (GeV/c)");
 367   fEtIsoCells->Sumw2();
 368   fOutput->Add(fEtIsoCells);
 369
 370   fEtIsoClust = new TH1D("hEtIsoClus_NC","#Sigma E_{T}^{iso cone} in iso cone distribution for Neutral Clusters with EMCAL Clusters",100,0.,100.);
 371   fEtIsoClust->SetXTitle("#Sigma E_{T}^{iso cone} (GeV/c)");
 372   fEtIsoClust->Sumw2();
 373   fOutput->Add(fEtIsoClust);
 374
 375   fPtIsoTrack = new TH1D("hPtIsoTrack_NC"," #Sigma P_{T}^{iso cone} in iso cone distribution for Neutral Clusters with Tracks",100,0.,100.);
 376   fPtIsoTrack->SetXTitle("#Sigma P_{T}^{iso cone} (GeV/c)");
 377   fPtIsoTrack->Sumw2();
 378   fOutput->Add(fPtIsoTrack);
 379
 380   fPtEtIsoTC = new TH1D("hPtEtIsoTrackClust_NC","#Sigma P_{T}^{iso cone} + #Sigma E_{T}^{iso cone} in iso cone distribution for Neutral Clusters with Tracks and Clusters",100,0.,100.);
 381   fPtEtIsoTC->SetXTitle("#Sigma P_{T}^{iso cone} + #Sigma E_{T}^{iso cone} (GeV/c)");
 382   fPtEtIsoTC->Sumw2();
 383   fOutput->Add(fPtEtIsoTC);
 384
 385   fPhiBandUEClust = new TH2D(Form("hPhiBandUE_Cluster"),Form("UE Estimation with Phi Band Clusters"),100,0.,100.,100,0.,100.);
 386   fPhiBandUEClust->SetXTitle("E_{T}");
 387   fPhiBandUEClust->SetYTitle("#Sigma E_{T}^{UE}");
 388   fPhiBandUEClust->Sumw2();
 389   fOutput->Add(fPhiBandUEClust);
 390
 391   fEtaBandUEClust = new TH2D(Form("hEtaBandUE_Cluster"),Form("UE Estimation with Phi Band Clusters"),100,0.,100.,100,0.,100.);
 392   fEtaBandUEClust->SetXTitle("E_{T}");
 393   fEtaBandUEClust->SetYTitle("#Sigma E_{T}^{UE}");
 394   fEtaBandUEClust->Sumw2();
 395   fOutput->Add(fEtaBandUEClust);
 396
 397   fPhiBandUECells = new TH2D(Form("hPhiBandUE_CELLS"),Form("UE Estimation with Phi Band CELLS"),100,0.,100.,100,0.,100.);
 398   fPhiBandUECells->SetXTitle("E_{T}");
 399   fPhiBandUECells->SetYTitle("#Sigma E_{T}^{UE}");
 400   fPhiBandUECells->Sumw2();
 401   fOutput->Add(fPhiBandUECells);
 402
 403   fEtaBandUECells = new TH2D(Form("hEtaBandUE_CELLS"),Form("UE Estimation with Phi Band and CELLS"),100,0.,100.,100,0.,100.);
 404   fEtaBandUECells->SetXTitle("E_{T}");
 405   fEtaBandUECells->SetYTitle("#Sigma E_{T}^{UE}");
 406   fEtaBandUECells->Sumw2();
 407   fOutput->Add(fEtaBandUECells);
 408
 409   fPhiBandUETracks = new TH2D(Form("hPhiBandUE_TPC"),Form("UE Estimation with Phi Band TPC "),100,0.,100.,100,0.,100.);
 410   fPhiBandUETracks->SetXTitle("E_{T}");
 411   fPhiBandUETracks->SetYTitle("#Sigma P_{T}^{UE}");
 412   fPhiBandUETracks->Sumw2();
 413   fOutput->Add(fPhiBandUETracks);
 414
 415   fEtaBandUETracks = new TH2D(Form("hEtaBandUE_TPC"),Form("UE Estimation with Phi Band and TPC"),100,0.,100.,100,0.,100.);
 416   fEtaBandUETracks->SetXTitle("E_{T}");
 417   fEtaBandUETracks->SetYTitle("#Sigma P_{T}^{UE}");
 418   fEtaBandUETracks->Sumw2();
 419   fOutput->Add(fEtaBandUETracks);
 420
 421   fPerpConesUETracks = new TH2D("hConesUE","UE Estimation with Perpendicular Cones in TPC",100,0.,100.,100,0.,100.);
 422   fPerpConesUETracks->SetXTitle("E_{T}");
 423   fPerpConesUETracks->SetYTitle("#Sigma P_{T}^{UE}");
 424   fPerpConesUETracks->Sumw2();
 425   fOutput->Add(fPerpConesUETracks);
 426
 427   fTPCWithoutIsoConeB2BbandUE = new TH2D("hFullTPCUE","UE Estimation with almost Full TPC",100,0.,100.,100,0.,100.);
 428   fPhiBandUEClust->SetXTitle("E_{T}");
 429   fPhiBandUEClust->SetYTitle("#Sigma E_{T}^{UE}");
 430   fTPCWithoutIsoConeB2BbandUE->Sumw2();
 431   fOutput->Add(fTPCWithoutIsoConeB2BbandUE);
 432
 433   fEtIsolatedClust = new TH1D("hEtIsolatedClust","E_{T} distribution for Isolated Photons with clusters; #Sigma E_{T}^{iso cone}<Ethres",100,0.,100.);
 434   fEtIsolatedClust->SetXTitle("E_{T}^{iso}");
 435   fEtIsolatedClust->Sumw2();
 436   fOutput->Add(fEtIsolatedClust);
 437
 438   fEtIsolatedCells = new TH1D("hEtIsolatedCells","E_{T} distribution for Isolated Photons with cells; #Sigma E_{T}^{iso cone}<Ethres",100,0.,100.);
 439   fEtIsolatedCells->SetXTitle("E_{T}^{iso}");
 440   fEtIsolatedCells->Sumw2();
 441   fOutput->Add(fEtIsolatedCells);
 442
 443   fEtIsolatedTracks = new TH1D("hEtIsolatedTracks","E_{T} distribution for Isolated Photons with tracks; #Sigma P_{T}^{iso cone}<Pthres",100,0.,100.);
 444   fEtIsolatedTracks->SetXTitle("E_{T}^{iso}");
 445   fEtIsolatedTracks->Sumw2();
 446   fOutput->Add(fEtIsolatedTracks);
 447
 448   fTest = new TH1D ("hTest","test cluster collection",100,-2.,6.);
 449   fTest->Sumw2();
 450   fOutput->Add(fTest);
 451
 452
 453   PostData(1, fOutput);
 454   //   return;
 455 }
 456
 457 //________________________________________________________________________
 458 Float_t* AliAnalysisTaskEMCALPhotonIsolation::GenerateFixedBinArray(Int_t n, Float_t min, Float_t max) const
 459 {
 460   // Generate the bin array for the ThnSparse
 461
 462   Float_t *bins = new Float_t[n+1];
 463
 464   Float_t binWidth = (max-min)/n;
 465   bins[0] = min;
 466   for (Int_t i = 1; i <= n; i++) {
 467     bins[i] = bins[i-1]+binWidth;
 468   }
 469
 470   return bins;
 471 }
 472
 473 //________________________________________________________________________
 474 void AliAnalysisTaskEMCALPhotonIsolation::ExecOnce()
 475 {
 476   //   Init the analysis.
 477
 478
 479
 480   if (fParticleCollArray.GetEntriesFast()<2) {
 481     AliError(Form("Wrong number of particle collections (%d), required 2",fParticleCollArray.GetEntriesFast()));
 482     return;
 483   }
 484
 485
 486   //  for (Int_t i = 0; i < 2; i++) {
 487   //    AliParticleContainer *contain = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
 488   //    // contain->GetClassName("AliEmcalParticle");
 489   //  }
 490
 491
 492
 493   AliAnalysisTaskEmcal::ExecOnce();
 494   if (!fInitialized) {
 495
 496     AliError(Form("AliAnalysisTask not initialized"));
 497     return;
 498   }
 499
 500   fevents+=1;
 501
 502 }
 503
 504 //______________________________________________________________________________________
 505 Bool_t AliAnalysisTaskEMCALPhotonIsolation::Run()
 506 {
 507   // Run the analysis
 508
 509
 510   AliParticleContainer *clusters = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
 511   AliEmcalParticle *emccluster = 0;
 512
 513   // delete output USEFUL LATER FOR CONTAINER CREATION !!
 514   //fOutClusters->Delete();
 515
 516   //Int_t clusCount = 0;    AliError(Form("Should be here each time"));
 517   // loop over all clusters
 518   clusters->ResetCurrentID();
 519   Int_t index=-1;
 520
 521   //Double_t ETleadingclust = 0., M02leadingcluster = 0., lambda0cluster = 0., phileadingclust = 0., etaleadingclust = 0., ptmc = 0.,mcptsum = 0.;
 522   //Int_t Ntracks;
 523   //Definition of the Array for Davide's Output
 524   //const Int_t ndims =   fNDimensions;
 525   //Double_t outputValues[ndims];
 526
 527
 528   if (fisLCAnalysis) {
 529     //   AliError(Form("Check the loop"));
 530     // get the leading particle
 531
 532
 533     emccluster = static_cast<AliEmcalParticle*>(clusters->GetLeadingParticle());
 534
 535     if(!emccluster){
 536
 537       AliError(Form("no leading one"));
 538       return kTRUE;
 539     }
 540
 541
 542
 543     //     emccluster = clusters->GetLeadingParticle();
 544     index = emccluster->IdInCollection();
 545
 546     //add a command to get the index of the leading cluster!
 547     if (!emccluster->IsEMCAL()) return kTRUE;
 548
 549     AliVCluster *coi = emccluster->GetCluster();
 550     if (!coi) return kTRUE;
 551
 552     TLorentzVector vecCOI;
 553     coi->GetMomentum(vecCOI,fVertex);
 554
 555
 556     if(!CheckBoundaries(vecCOI))
 557       return kTRUE;
 558
 559     if(ClustTrackMatching(coi))
 560       return kTRUE;
 561
 562     else
 563       FillGeneralHistograms(coi,vecCOI, index);
 564   }
 565   else{
 566     //get the entries of the Cluster Container
 567     //whatever is a RETURN in LCAnalysis here is a CONTINUE,
 568     //since there are more than 1 Cluster per Event
 569
 570     while((emccluster = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle()))){
 571       index++;
 572       if (!emccluster->IsEMCAL()) return kTRUE;
 573
 574       AliVCluster *coi = emccluster->GetCluster();
 575       if(!coi) return kTRUE;
 576
 577       TLorentzVector vecCOI;
 578       coi->GetMomentum(vecCOI,fVertex);
 579
 580       if(!CheckBoundaries(vecCOI))
 581         return kTRUE;
 582
 583       if(ClustTrackMatching(coi))
 584         continue;
 585       else
 586         FillGeneralHistograms(coi,vecCOI, index);
 587
 588     }
 589   }
 590   //  PostData(1, fOutput);
 591   return kTRUE;
 592 }
 593
 594
 595 //__________________________________________________________________________
 596 Bool_t AliAnalysisTaskEMCALPhotonIsolation::ClustTrackMatching(AliVCluster *cluster) {
 597   // Check if the cluster match to a track
 598
 599   Double_t deta=999;
 600   Double_t dphi=999;
 601
 602   AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
 603   // for an incoming cluster reject it if there is a track corresponding with a deta and dphi lower than the cuts
 604   TLorentzVector nPart;
 605   cluster->GetMomentum(nPart, fVertex);
 606
 607   AliVTrack *mt = NULL;
 608   AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
 609   if(acl) {
 610     if(acl->GetNTracksMatched()>1)
 611       mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
 612   }
 613   else {
 614     AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
 615     if(!ecl){
 616       AliError("ClusTrack matching did not work");
 617       return kTRUE;
 618     }
 619     Int_t im = ecl->GetTrackMatchedIndex();
 620     if(tracks && im>=0) {
 621       mt = static_cast<AliVTrack*>(tracks->GetParticle(im));
 622     }
 623   }
 624   //  if(mt && mt->TestFilterBit(768)) { //Hybrid tracks with AOD
 625   AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
 626   fDeltaETAClusTrackVSpT->Fill(nPart.Pt(), deta);
 627   fDeltaPHIClusTrackVSpT->Fill(nPart.Pt(), dphi);
 628
 629
 630   if(mt) return kTRUE;
 631
 632   else return kFALSE;
 633
 634 }
 635
 636
 637
 638 //__________________________________________________________________________
 639 void AliAnalysisTaskEMCALPhotonIsolation::EtIsoCellPhiBand(TLorentzVector c, Float_t &etIso, Float_t &phiBandcells){
 640   // Underlying events study with EMCAL cells in phi band
 641
 642   AliEMCALGeometry* emcalGeom = AliEMCALGeometry::GetInstance();
 643
 644   Float_t sumEnergyPhiBandCells=0., sumEnergyConeCells=0.;
 645
 646
 647   // check the cell corresponding to the leading cluster
 648   Int_t absId = 999;
 649   //maybe best to call it LeadingCellIdinClus or better maxId since it won't be used anywhere else ???
 650   Bool_t cellLeadingClustID = emcalGeom->GetAbsCellIdFromEtaPhi(c.Eta(),c.Phi(),absId);
 651   if(!cellLeadingClustID) return;
 652
 653   else{
 654     Int_t iTower = -1;
 655     Int_t iModule = -1;
 656     Int_t imEta=-1, imPhi=-1;
 657     Int_t iEta =-1, iPhi =-1;
 658
 659     emcalGeom->GetCellIndex(absId,iModule,iTower,imPhi,imEta); // to get the module, the tower, eta and phi for the cell corresponding to the leading cluster
 660     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
 661
 662     // Get the row and the column of the cell corresponding to the leading cluster in EMCAL
 663     Int_t colCellLeadingClust = iEta;
 664     if(iModule % 2) colCellLeadingClust = AliEMCALGeoParams::fgkEMCALCols + iEta ; // if the SM number is even you need to translate to have the one corresponding in EMCAL
 665     Int_t rowCellLeadingClust = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(iModule/2); // to have the corresponding row in EMCAL
 666
 667     // total number or rows and columns in EMCAL
 668     Int_t nTotalRows = AliEMCALGeoParams::fgkEMCALRows*16/3 ;  // 5 + 2/3 supermodules in a row
 669     Int_t nTotalCols = 2*AliEMCALGeoParams::fgkEMCALCols; // 2 supermodules in a column
 670
 671     Int_t nbConeSize = int(fIsoConeRadius/0.0143); //0.0143 tower acceptance
 672
 673     // Get the cells
 674     AliVCaloCells * cells =InputEvent()->GetEMCALCells();
 675
 676     // define the max and min row and column corresponding to the isolation cone around the seed cell from the leading cluster
 677     Int_t iRowMinCone = rowCellLeadingClust-nbConeSize;
 678     if(iRowMinCone<0) iRowMinCone=0;
 679
 680     Int_t iRowMaxCone = rowCellLeadingClust+nbConeSize;
 681     if(iRowMaxCone>AliEMCALGeoParams::fgkEMCALRows) iRowMaxCone=AliEMCALGeoParams::fgkEMCALRows;  // AliEMCALGeoParams::fgkEMCALRows = 24 in a supermodule
 682
 683     Int_t iColMinCone = colCellLeadingClust - nbConeSize;
 684     if(iColMinCone<0) iColMinCone = 0;
 685
 686     Int_t iColMaxCone = colCellLeadingClust+nbConeSize;
 687     if(iColMaxCone>AliEMCALGeoParams::fgkEMCALCols) iColMaxCone=AliEMCALGeoParams::fgkEMCALCols;  // AliEMCALGeoParams::fgkEMCALCols = 48 in a supermodule
 688
 689     // loop on all cells
 690     for(Int_t iCol=0; iCol<nTotalCols; iCol++){
 691       for(Int_t iRow=0; iRow<nTotalRows; iRow++){
 692         // now recover the cell indexes in a supermodule
 693         Int_t iSector = int(iRow/AliEMCALGeoParams::fgkEMCALRows); // check in which SM is the cell
 694         if(iSector==5) continue;  //
 695         Int_t inModule = -1;
 696         Int_t iColLoc  = -1;
 697         if(iCol < AliEMCALGeoParams::fgkEMCALCols){ // if the SM number is odd the column is the one corresponding in the supermodule
 698           inModule = 2*iSector + 1;
 699           iColLoc  = iCol;
 700         }
 701         else if(iCol > AliEMCALGeoParams::fgkEMCALCols - 1){ // if the SM number is even the column isn't the one corresponding in the supermodule
 702           inModule = 2*iSector;
 703           iColLoc  = iCol-AliEMCALGeoParams::fgkEMCALCols;
 704         }
 705
 706         Int_t iRowLoc  = iRow - AliEMCALGeoParams::fgkEMCALRows*iSector ;
 707
 708         if(TMath::Abs(iCol-colCellLeadingClust)<nbConeSize && TMath::Abs(iCol+colCellLeadingClust)>nbConeSize){
 709           if(iRow<iRowMaxCone && iRow>iRowMinCone){
 710             Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(inModule,iRowLoc,iColLoc);  // verifier les iRow et iCol
 711             sumEnergyPhiBandCells+=cells->GetAmplitude(iabsId); //should be Et
 712           }
 713         }
 714         else if (TMath::Abs(iCol-colCellLeadingClust)>nbConeSize && TMath::Abs(iCol+colCellLeadingClust)<nbConeSize){
 715           Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(inModule,iRowLoc,iColLoc);  // verifier les iRow et iCol
 716           sumEnergyConeCells+=cells->GetAmplitude(iabsId); //should be Et
 717         }
 718       }
 719     }
 720   }
 721   etIso = sumEnergyConeCells;
 722   phiBandcells = sumEnergyPhiBandCells;
 723 }
 724
 725
 726 //__________________________________________________________________________
 727 void AliAnalysisTaskEMCALPhotonIsolation::EtIsoCellEtaBand(TLorentzVector c, Float_t &etIso, Float_t &etaBandcells){
 728   // Underlying events study with EMCAL cell in eta band
 729
 730
 731   AliEMCALGeometry* emcalGeom = AliEMCALGeometry::GetInstance();
 732
 733   Float_t sumEnergyEtaBandCells=0., sumEnergyConeCells=0.;
 734
 735
 736
 737   // check the cell corresponding to the leading cluster
 738   Int_t absId = 999;
 739   //maybe best to call it LeadingCellIdinClus or better maxId since it won't be used anywhere else ???
 740   Bool_t cellLeadingClustID = emcalGeom->GetAbsCellIdFromEtaPhi(c.Eta(),c.Phi(),absId);
 741   if(!cellLeadingClustID) return;
 742
 743   else{
 744
 745     Int_t iTower = -1;
 746     Int_t iModule = -1;
 747     Int_t imEta=-1, imPhi=-1;
 748     Int_t iEta =-1, iPhi =-1;
 749
 750     emcalGeom->GetCellIndex(absId,iModule,iTower,imPhi,imEta); // to get the module, the tower, eta and phi for the cell corresponding to the leading cluster
 751     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
 752
 753     // Get the row and the column of the cell corresponding to the leading cluster in EMCAL
 754     Int_t colCellLeadingClust = iEta;
 755     if(iModule % 2) colCellLeadingClust = AliEMCALGeoParams::fgkEMCALCols + iEta ; // if the SM number is even you need to translate to have the one corresponding in EMCAL
 756     Int_t rowCellLeadingClust = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(iModule/2); // to have the corresponding row in EMCAL
 757
 758     // total number or rows and columns in EMCAL
 759     Int_t nTotalRows = AliEMCALGeoParams::fgkEMCALRows*16/3 ;  // 5 + 2/3 supermodules in a row
 760     Int_t nTotalCols = 2*AliEMCALGeoParams::fgkEMCALCols; // 2 supermodules in a column
 761
 762     Int_t nbConeSize = int(fIsoConeRadius/0.0143); //0.0143 tower acceptance
 763
 764     // Get the cells
 765     AliVCaloCells * cells =InputEvent()->GetEMCALCells();
 766
 767     // define the max and min row and column corresponding to the isolation cone around the seed cell from the leading cluster
 768     Int_t iRowMinCone = rowCellLeadingClust-nbConeSize;
 769     if(iRowMinCone<0) iRowMinCone=0;
 770
 771     Int_t iRowMaxCone = rowCellLeadingClust+nbConeSize;
 772     if(iRowMaxCone>AliEMCALGeoParams::fgkEMCALRows) iRowMaxCone=AliEMCALGeoParams::fgkEMCALRows;  // AliEMCALGeoParams::fgkEMCALRows = 24 in a supermodule
 773
 774     Int_t iColMinCone = colCellLeadingClust-nbConeSize;
 775     if(iColMinCone<0) iColMinCone = 0;
 776
 777     Int_t iColMaxCone = colCellLeadingClust+nbConeSize;
 778     if(iColMaxCone>AliEMCALGeoParams::fgkEMCALCols) iColMaxCone=AliEMCALGeoParams::fgkEMCALCols;  // AliEMCALGeoParams::fgkEMCALCols = 48 in a supermodule
 779
 780     // loop on all cells
 781     for(Int_t iCol=0; iCol<nTotalCols; iCol++)
 782     {
 783       for(Int_t iRow=0; iRow<nTotalRows; iRow++)
 784       {
 785         // now recover the cell indexes in a supermodule
 786         Int_t iSector = int(iRow/AliEMCALGeoParams::fgkEMCALRows); // check in which SM is the cell
 787         if(iSector==5) continue;  //
 788         Int_t inModule = -1;
 789         Int_t iColLoc  = -1;
 790         if(iCol < AliEMCALGeoParams::fgkEMCALCols){ // if the SM number is odd the column is the one corresponding in the supermodule
 791           inModule = 2*iSector + 1;
 792           iColLoc  = iCol;
 793         }
 794         else if(iCol > AliEMCALGeoParams::fgkEMCALCols - 1){ // if the SM number is even the column isn't the one corresponding in the supermodule
 795           inModule = 2*iSector;
 796           iColLoc  = iCol-AliEMCALGeoParams::fgkEMCALCols;
 797         }
 798
 799         Int_t iRowLoc  = iRow - AliEMCALGeoParams::fgkEMCALRows*iSector ;
 800
 801         if(TMath::Abs(iCol-colCellLeadingClust)<nbConeSize && TMath::Abs(iCol+colCellLeadingClust)>nbConeSize){
 802           if(iCol<iColMaxCone && iCol>iColMinCone){
 803             Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(inModule,iRowLoc,iColLoc);  // verifier les iRow et iCol
 804             sumEnergyEtaBandCells+=cells->GetAmplitude(iabsId); //should be Et
 805           }
 806         }
 807         else if (TMath::Abs(iCol-colCellLeadingClust)>nbConeSize && TMath::Abs(iCol+colCellLeadingClust)<nbConeSize){
 808           Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(inModule,iRowLoc,iColLoc);  // verifier les iRow et iCol
 809           sumEnergyConeCells+=cells->GetAmplitude(iabsId); //should be Et
 810         }
 811       }
 812     }
 813   }
 814   etIso = sumEnergyConeCells;
 815   etaBandcells = sumEnergyEtaBandCells;
 816 }
 817
 818
 819 //__________________________________________________________________________
 820 void AliAnalysisTaskEMCALPhotonIsolation::EtIsoClusPhiBand(TLorentzVector c, Float_t &etIso, Float_t &phiBandclus, Int_t index){
 821   // Underlying events study with clusters in phi band
 822
 823   Float_t sumEnergyPhiBandClus=0., sumEnergyConeClus=0.;
 824
 825   //needs a check on the same cluster
 826   AliParticleContainer *clusters = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
 827   AliEmcalParticle *clust;
 828
 829   Int_t localIndex=-1;
 830   while((clust = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle()))){ //check the position of other clusters in respect to the leading cluster
 831     localIndex++;
 832     if(localIndex==index) continue;
 833
 834     AliVCluster *cluster= clust->GetCluster();
 835
 836     TLorentzVector nClust; //STILL NOT INITIALIZED
 837     cluster->GetMomentum(nClust,fVertex);
 838     Float_t phiClust =nClust.Phi();
 839     Float_t etaClust= nClust.Eta();
 840     //redefine phi/c.Eta() from the cluster we passed to the function
 841
 842     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
 843
 844     if(radius>fIsoConeRadius){ // the cluster is outside the isolation cone in this case study UE
 845
 846       // actually phi band here
 847       if(TMath::Abs(etaClust - c.Eta()) < fIsoConeRadius){
 848         sumEnergyPhiBandClus += nClust.Pt();
 849       }
 850     }
 851     else // if the cluster is in the isolation cone, add the cluster pT
 852       sumEnergyConeClus += nClust.Pt();
 853
 854   }
 855   etIso = sumEnergyConeClus;
 856   phiBandclus = sumEnergyPhiBandClus;
 857 }
 858
 859
 860 //__________________________________________________________________________
 861 void AliAnalysisTaskEMCALPhotonIsolation::EtIsoClusEtaBand(TLorentzVector c, Float_t &etIso, Float_t &etaBandclus, Int_t index){
 862   // Underlying events study with clusters in eta band
 863
 864   Float_t sumEnergyEtaBandClus =0., sumEnergyConeClus=0.;
 865   AliParticleContainer *clusters = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
 866
 867   AliEmcalParticle *clust;
 868
 869   Int_t localIndex=-1;
 870   while((clust = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle()))){ //check the position of other clusters in respect to the leading cluster
 871     localIndex++;
 872     if(localIndex==index) continue;
 873
 874     AliVCluster *cluster= clust->GetCluster();
 875
 876     TLorentzVector nClust; //STILL NOT INITIALIZED
 877     cluster->GetMomentum(nClust,fVertex);
 878
 879     Float_t phiClust =nClust.Phi();
 880     Float_t etaClust= nClust.Eta();
 881     //redefine phi/c.Eta() from the cluster we passed to the function
 882
 883     // define the radius between the leading cluster and the considered cluster
 884     Float_t  radius = TMath::Sqrt(TMath::Power(phiClust-c.Phi(),2)+TMath::Power(etaClust-c.Eta(),2));
 885
 886     if(radius>fIsoConeRadius){ // the cluster is outside the isolation cone in this case study UE
 887
 888       // actually eta band here
 889       if(TMath::Abs(etaClust - c.Eta()) < fIsoConeRadius){
 890         sumEnergyEtaBandClus += nClust.Pt();
 891       }
 892     }
 893     else  // if the cluster is in the isolation cone, add the cluster pT
 894       sumEnergyConeClus += nClust.Pt();
 895
 896   }
 897   etIso = sumEnergyConeClus;
 898   etaBandclus = sumEnergyEtaBandClus;
 899 }
 900
 901
 902 //__________________________________________________________________________
 903 void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackPhiBand(TLorentzVector c, Float_t &ptIso, Float_t &phiBandtrack){
 904   // Underlying events study with tracks in phi band in EMCAL acceptance
 905
 906   //INSERT BOUNDARIES ACCORDING TO THE FLAG WE WANT TO USE
 907   Float_t sumpTConeCharged=0., sumpTPhiBandTrack=0.;
 908   Float_t minPhi= 0., maxPhi= 2*TMath::Pi(), minEta = -0.9, maxEta= 0.9;
 909
 910   AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
 911
 912   if(!fTPC4Iso){
 913     minEta = -0.7;
 914     maxEta = 0.7;
 915     minPhi = 1.4;
 916     maxPhi = TMath::Pi();
 917   }
 918
 919   AliVParticle *track = tracks->GetNextAcceptParticle(0);
 920
 921   while(track){
 922
 923     //CHECK IF TRACK IS IN BOUNDARIES
 924     Float_t phiTrack = track->Phi();
 925     Float_t etaTrack = track->Eta();
 926     // define the radius between the leading cluster and the considered cluster
 927     //redefine phi/c.Eta() from the cluster we passed to the function
 928     if(phiTrack < maxPhi && phiTrack > minPhi && etaTrack < maxEta && etaTrack > minEta){
 929       Float_t  radius = TMath::Sqrt(TMath::Power(phiTrack - c.Phi(),2)+TMath::Power(etaTrack - c.Eta(),2));
 930
 931       if(radius>fIsoConeRadius){ // if tracks are outside the isolation cone study
 932
 933         // actually phi band here --- ADD Boundaries conditions
 934         if(TMath::Abs(etaTrack - c.Eta()) < fIsoConeRadius){
 935           sumpTPhiBandTrack += track->Pt();
 936         }
 937       }
 938       else
 939         sumpTConeCharged+=track->Pt(); // should not double count if the track matching is already done
 940     }
 941     track=tracks->GetNextAcceptParticle();
 942   }
 943   ptIso = sumpTConeCharged;
 944   phiBandtrack = sumpTPhiBandTrack;
 945 }
 946
 947
 948 //__________________________________________________________________________
 949 void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackEtaBand(TLorentzVector c, Float_t &ptIso, Float_t &etaBandtrack){
 950   // Underlying events study with tracks in eta band in EMCAL acceptance
 951
 952   //INSERT BOUNDARIES ACCORDING TO THE FLAG WE WANT TO USE
 953   Float_t sumpTConeCharged=0., sumpTEtaBandTrack=0.;
 954   Float_t minPhi= 0., maxPhi= 2*TMath::Pi(), minEta = -0.9, maxEta= 0.9;
 955
 956   AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
 957
 958
 959   if(!fTPC4Iso){
 960     minEta = -0.7;
 961     maxEta = 0.7;
 962     minPhi = 1.4;
 963     maxPhi = TMath::Pi();
 964   }
 965
 966
 967   AliVParticle *track = tracks->GetNextAcceptParticle(0);
 968
 969   while(track){
 970
 971     Float_t phiTrack = track->Phi();
 972     Float_t etaTrack = track->Eta();
 973     //redefine phi/c.Eta() from the cluster we passed to the function
 974     if(phiTrack < maxPhi && phiTrack > minPhi && etaTrack < maxEta && etaTrack > minEta){
 975       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
 976       if(radius>fIsoConeRadius){ // if tracks are outside the isolation cone study UE
 977
 978         // actually eta band here --- ADD Boundaries conditions
 979         if(TMath::Abs(phiTrack - c.Phi()) < fIsoConeRadius){
 980           sumpTEtaBandTrack += track->Pt();
 981         }
 982       }
 983       else
 984         sumpTConeCharged += track->Pt(); // should not double count if the track matching is already done
 985     }
 986     track=tracks->GetNextAcceptParticle();
 987   }
 988   ptIso = sumpTConeCharged;
 989   etaBandtrack = sumpTEtaBandTrack;
 990 }
 991
 992
 993 //__________________________________________________________________________
 994 void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackOrthCones(TLorentzVector c, Float_t &ptIso, Float_t &cones){
 995   // Underlying events study with tracks in orthogonal cones in TPC
 996
 997   Float_t sumpTConeCharged=0., sumpTPerpConeTrack=0.;
 998
 999   AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
1000
1001   Float_t etaClus = c.Eta();
1002   Float_t phiClus = c.Phi();
1003   Float_t phiCone1 = phiClus - TMath::PiOver2();
1004   Float_t phiCone2 = phiClus + TMath::PiOver2();
1005
1006   if (phiCone1 < 0.) phiCone1 += 2*TMath::Pi();
1007
1008   AliVParticle *track = tracks->GetNextAcceptParticle(0);
1009
1010   while(track){
1011
1012     Float_t phiTrack = track->Phi();
1013     Float_t etaTrack = track->Eta();
1014
1015     Float_t dist2Clust = TMath::Sqrt(TMath::Power(etaTrack-etaClus, 2)+TMath::Power(phiTrack-phiClus, 2));
1016     if (dist2Clust<fIsoConeRadius)// if tracks are inside the IsoCone
1017       sumpTConeCharged += track->Pt();
1018
1019     else{//tracks outside the IsoCone
1020       //Distances from the centres of the two Orthogonal Cones
1021       Float_t dist2Cone1 = TMath::Sqrt(TMath::Power(etaTrack-etaClus, 2)+TMath::Power(phiTrack-phiCone1, 2));
1022       Float_t dist2Cone2 = TMath::Sqrt(TMath::Power(etaTrack-etaClus, 2)+TMath::Power(phiTrack-phiCone2, 2));
1023
1024       //Is the Track Inside one of the two Cones ->Add to UE
1025       if((dist2Cone1 < fIsoConeRadius) || (dist2Cone2 < fIsoConeRadius))
1026         sumpTPerpConeTrack += track->Pt();
1027     }
1028     track=tracks->GetNextAcceptParticle();
1029   }
1030   ptIso = sumpTConeCharged;
1031   cones = sumpTPerpConeTrack;
1032 }
1033
1034 //__________________________________________________________________________
1035 void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackFullTPC(TLorentzVector c, Float_t &ptIso, Float_t &full){
1036   // Underlying events study with tracks in full TPC except back to back bands
1037
1038   Float_t sumpTConeCharged=0., sumpTTPCexceptB2B=0.;
1039
1040   AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
1041   AliVParticle *track = tracks->GetNextAcceptParticle(0);
1042
1043   while(track){
1044
1045     Float_t phiTrack = track->Phi();
1046     Float_t etaTrack = track->Eta();
1047     //redefine phi/c.Eta() from the cluster we passed to the function
1048
1049     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
1050     if(radius>fIsoConeRadius){ // if tracks are outside the isolation cone study UE
1051       Double_t dphiUp = c.Phi() + TMath::Pi() - fIsoConeRadius;
1052       Double_t dphiDown = c.Phi() + TMath::Pi() + fIsoConeRadius;
1053       // TPC except B2B
1054       if(TMath::Power(phiTrack-c.Phi(),2) +TMath::Power(etaTrack-c.Eta(),2) > TMath::Power((fIsoConeRadius+0.1),2) && phiTrack < dphiDown && phiTrack> dphiUp)
1055         sumpTTPCexceptB2B += track->Pt();
1056     }
1057     else
1058       sumpTConeCharged += track->Pt(); // should not double count if the track matching is already done
1059
1060     track=tracks->GetNextAcceptParticle();
1061   }
1062   ptIso = sumpTConeCharged;
1063   full = sumpTTPCexceptB2B;
1064 }
1065
1066 //__________________________________________________________________________
1067 Bool_t AliAnalysisTaskEMCALPhotonIsolation::CheckBoundaries(TLorentzVector vecCOI){
1068   // Check if the cone around the considered cluster is in EMCAL acceptance
1069
1070   Float_t minPhiBound= 1.5 , minEtaBound= -0.5, maxPhiBound= TMath::Pi()-0.1, maxEtaBound= 0.5;
1071   Bool_t isINBoundaries;
1072
1073   if(!fTPC4Iso){
1074     minPhiBound = 1.4+0.4; //to be changed with fIsoConeR
1075     maxPhiBound = TMath::Pi()-0.4;
1076     minEtaBound = -0.7+0.4;
1077     maxEtaBound = 0.7-0.4;
1078   }
1079
1080   if(vecCOI.Eta() > maxEtaBound || vecCOI.Eta() < minEtaBound || vecCOI.Phi() > maxPhiBound || vecCOI.Phi() >minPhiBound)
1081     isINBoundaries=kFALSE;
1082   else
1083     isINBoundaries=kTRUE;
1084
1085   return isINBoundaries;
1086 }
1087
1088 //__________________________________________________________________________
1089 Bool_t AliAnalysisTaskEMCALPhotonIsolation::FillGeneralHistograms(AliVCluster *coi, TLorentzVector vecCOI, Int_t index){
1090   // Fill the histograms for underlying events and isolation studies
1091
1092   AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
1093
1094   AliEmcalParticle *emcTrack = 0;
1095
1096   int nTracks=0;
1097   tracks->ResetCurrentID();
1098   while ((emcTrack = static_cast<AliEmcalParticle*>(tracks->GetNextAcceptParticle()))) {
1099     AliVTrack *track = emcTrack->GetTrack();
1100     if(!track) continue;
1101     // if(!(track->TestFilterBit("kHybrid"))) continue;
1102
1103     nTracks++;
1104   }
1105
1106   Double_t eTCOI = 0., m02COI = 0., lambda0cluster = 0., ptmc = 0., mcptsum = 0.;
1107   //Int_t Ntracks;
1108   //Definition of the Array for Davide's Output
1109   const Int_t ndims =   fNDimensions;
1110   Double_t outputValues[ndims];
1111
1112   eTCOI = vecCOI.Et();
1113   m02COI = coi->GetM02();
1114
1115
1116   // ******** Isolation and UE calculation with different methods *********
1117
1118   Double_t eTThreshold = 5;
1119
1120   switch(fEtIsoMethod)
1121   {
1122     case 0:  // SumEt<EtThr
1123       if(fM02mincut < m02COI && m02COI < fM02maxcut)  // photon candidate, cuts have to be decided after studies
1124       {
1125         eTThreshold = fEtIsoThreshold;
1126       }
1127       break;
1128
1129     case 1:  // SumEt<%Ephoton
1130       if(fM02mincut < m02COI && m02COI < fM02maxcut) // photon candidate, cuts have to be decided after studies
1131       {
1132         eTThreshold = fEtIsoThreshold * eTCOI;
1133       }
1134       break;
1135
1136     case 2: // Etmax<eTThreshold
1137       eTThreshold = fEtIsoThreshold;
1138       if(fM02mincut < m02COI && m02COI < fM02maxcut && eTCOI<eTThreshold) // photon candidate, cuts have to be decided after studies
1139       {
1140         fEtIsolatedClust->Fill(eTCOI);
1141       }
1142       break;
1143   }
1144
1145   //DO NOT CHANGE EVER AGAIN THE FOLLOWING DEFINITIONS
1146   Float_t isoConeArea = TMath::Pi()*fIsoConeRadius*fIsoConeRadius;
1147   Float_t etaBandArea = (1.4 - 2.*fIsoConeRadius)*2.*fIsoConeRadius;
1148   Float_t phiBandArea = ((5./9.)*TMath::Pi()- 2.*fIsoConeRadius)*2.*fIsoConeRadius;
1149   Float_t etaBandAreaTr = (1.8 - 2.*fIsoConeRadius)*2.*fIsoConeRadius;
1150   Float_t phiBandAreaTr = (2.*TMath::Pi() - 4.*fIsoConeRadius)*2.*fIsoConeRadius;//there is a 2 more because of the JET CONE B2B
1151   Float_t perpConesArea = 2.*isoConeArea;
1152   Float_t fullTPCArea = 1.8*2.*TMath::Pi()-fIsoConeRadius*(TMath::Pi()*fIsoConeRadius + 3.6);
1153
1154   Float_t isolation, ue;
1155
1156   if(!fTPC4Iso){ //EMCAL Only for Acceptance of Clusters
1157     switch(fIsoMethod)
1158     {
1159       case 0: //EMCAL CELLS
1160
1161         switch (fUEMethod)
1162       {
1163         case 0: //phi band
1164           EtIsoCellPhiBand(vecCOI, isolation, ue);
1165           //Normalisation ue wrt Area - TO DO-
1166           ue = ue * (isoConeArea / phiBandArea);
1167           fPhiBandUECells->Fill(vecCOI.Pt() , ue);
1168           fEtIsoCells->Fill(isolation);
1169           if(isolation<eTThreshold)
1170           {
1171             fEtIsolatedCells->Fill(eTCOI);
1172             fEtisoT=eTCOI;
1173             fPtisoT=vecCOI.Pt();
1174           }
1175           break;
1176         case 1: //eta band
1177           EtIsoCellEtaBand(vecCOI, isolation, ue);
1178           //Normalisation ue wrt Area - TO DO-
1179           ue = ue * (isoConeArea / etaBandArea);
1180           fEtaBandUECells->Fill(vecCOI.Pt() , ue);
1181           fEtIsoCells->Fill(isolation);
1182           if(isolation<eTThreshold)
1183           {
1184             fEtIsolatedCells->Fill(eTCOI);
1185             fEtisoT=eTCOI;
1186             fPtisoT=vecCOI.Pt();
1187           }
1188           break;
1189       }
1190         break;
1191
1192       case 1: //EMCAL CLUSTERS
1193
1194         switch (fUEMethod)
1195       {
1196         case 0: //phi band
1197           EtIsoClusPhiBand(vecCOI, isolation, ue,index);
1198           //Normalisation ue wrt Area - TO DO-
1199           ue = ue * (isoConeArea / phiBandArea);
1200           fPhiBandUEClust->Fill(vecCOI.Pt() , ue);
1201           fEtIsoClust->Fill(isolation);
1202           break;
1203         case 1: //eta band
1204           EtIsoClusEtaBand(vecCOI, isolation, ue,index);
1205           //Normalisation ue wrt Area - TO DO-
1206           ue = ue * (isoConeArea / etaBandArea);
1207           fEtaBandUEClust->Fill(vecCOI.Pt() , ue);
1208           fEtIsoClust->Fill(isolation);
1209           if(isolation<eTThreshold)
1210           {
1211             fEtIsolatedClust->Fill(eTCOI);
1212             fEtisoT=eTCOI;
1213             fPtisoT=vecCOI.Pt();
1214           }
1215           break;
1216       }
1217       case 2: //TRACKS (TBD which tracks) //EMCAL tracks
1218         switch (fUEMethod)
1219       {
1220         case 0: //phi band
1221           PtIsoTrackPhiBand(vecCOI, isolation, ue);
1222           //Normalisation ue wrt Area - TO DO-
1223           ue = ue * (isoConeArea / phiBandAreaTr);
1224           fPhiBandUETracks->Fill(vecCOI.Pt() , ue);
1225           fPtIsoTrack->Fill(isolation);
1226           if(isolation<eTThreshold)
1227           {
1228             fEtIsolatedTracks->Fill(eTCOI);
1229             fEtisoT=eTCOI;
1230             fPtisoT=vecCOI.Pt();
1231           }
1232           break;
1233         case 1: //eta band
1234           PtIsoTrackEtaBand(vecCOI, isolation, ue);
1235           //Normalisation ue wrt Area - TO DO-
1236           ue = ue * (isoConeArea / etaBandAreaTr);
1237           fEtaBandUETracks->Fill(vecCOI.Pt() , ue);
1238           fPtIsoTrack->Fill(isolation);
1239           if(isolation<eTThreshold)
1240           {
1241             fEtIsolatedTracks->Fill(eTCOI);
1242             fEtisoT=eTCOI;
1243             fPtisoT=vecCOI.Pt();
1244           }
1245           break;
1246           // case 2: //Cones
1247           // PtIsoTrackOrthCones(vecCOI, absId, isolation, ue);
1248           // break;
1249           // case 3: //Full TPC
1250           // PtIsoTrackFullTPC(vecCOI, absId, isolation, ue);
1251           // break;
1252       }
1253     }
1254   }
1255   else{  //EMCAL + TPC (Only tracks for the Isolation since IsoCone Goes Out of EMCAL)
1256     switch (fUEMethod)
1257     {
1258       case 0: //phi band
1259         PtIsoTrackPhiBand(vecCOI, isolation, ue);
1260         //Normalisation ue wrt Area - TO DO-
1261         ue = ue * (isoConeArea / phiBandAreaTr);
1262         fPhiBandUETracks->Fill(vecCOI.Pt() , ue);
1263         fPtIsoTrack->Fill(isolation);
1264         if(isolation<eTThreshold)
1265         {
1266           fEtIsolatedTracks->Fill(eTCOI);
1267           fEtisoT=eTCOI;
1268           fPtisoT=vecCOI.Pt();
1269         }
1270         break;
1271       case 1: //eta band
1272         PtIsoTrackEtaBand(vecCOI, isolation, ue);
1273         //Normalisation ue wrt Area - TO DO-
1274         ue = ue * (isoConeArea / etaBandAreaTr);
1275         fEtaBandUETracks->Fill(vecCOI.Pt() , ue);
1276         fPtIsoTrack->Fill(isolation);
1277         if(isolation<eTThreshold)
1278         {
1279           fEtIsolatedTracks->Fill(eTCOI);
1280           fEtisoT=eTCOI;
1281           fPtisoT=vecCOI.Pt();
1282         }
1283         break;
1284       case 2: //Cones
1285         PtIsoTrackOrthCones(vecCOI, isolation, ue);
1286         //Normalisation ue wrt Area - TO DO-
1287         ue = ue * (isoConeArea / perpConesArea);
1288         fPerpConesUETracks ->Fill(vecCOI.Pt() , ue);
1289         fPtIsoTrack->Fill(isolation);
1290         if(isolation<eTThreshold)
1291         {
1292           fEtIsolatedTracks->Fill(eTCOI);
1293           fEtisoT=eTCOI;
1294           fPtisoT=vecCOI.Pt();
1295         }
1296         break;
1297       case 3: //Full TPC
1298         PtIsoTrackFullTPC(vecCOI, isolation, ue);
1299         //Normalisation ue wrt Area - TO DO-
1300         ue = ue * (isoConeArea / fullTPCArea);
1301         fTPCWithoutIsoConeB2BbandUE->Fill(vecCOI.Pt() , ue);
1302         fPtIsoTrack->Fill(isolation);
1303         if(isolation<eTThreshold)
1304         {
1305           fEtIsolatedTracks->Fill(eTCOI);
1306           fEtisoT=eTCOI;
1307           fPtisoT=vecCOI.Pt();
1308         }
1309         break;
1310     }
1311   }
1312
1313
1314   //Here we should call something to know the number of tracks...
1315   //Soon I'll put in this version the "old way"; please let me know if
1316   //any of you could do the same with the JET framework
1317
1318   switch(fWho) {
1319     case 0:
1320       flambda0T=m02COI;
1321       fEtT=vecCOI.Et();
1322       fPtT=vecCOI.Pt();
1323       fetaT=vecCOI.Eta();
1324       fphiT=vecCOI.Phi();
1325       fsumEtisoconeT=isolation;
1326       //           AliError(Form("lambda 0 %f",flambda0T));
1327       fsumEtUE=ue;
1328
1329       fOutputTree->Fill();
1330       break;
1331
1332     case 1:
1333       outputValues[0] = nTracks;
1334       outputValues[1] = eTCOI;
1335       outputValues[2] = lambda0cluster;
1336       outputValues[3] = isolation;
1337       outputValues[4] = ue;
1338       outputValues[5] = isolation-ue;
1339       outputValues[6] = vecCOI.Eta();
1340       outputValues[7] = vecCOI.Phi();
1341       if (fIsMC) {
1342         outputValues[8] = ptmc;
1343         outputValues[9] = mcptsum;
1344       }
1345       fOutputTHnS -> Fill(outputValues);
1346       break;
1347       //                 //            fOutPTMAX -> Fill(outputValues[1],outputValues[2],);
1348   }
1349   return kTRUE;
1350 }
1351
1352
1353