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