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