]> 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
Fixed new version
[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
754 fAODMCParticles = dynamic_cast <TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
755 mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
756
757 AnalyzeMC();
758 }
759
1c662fe8 760 if (fisLCAnalysis) {
7c97fe2d 761
762 // get the leading particle
763 AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(clusters->GetLeadingParticle());
764
1c662fe8 765 if(!emccluster){
7c97fe2d 766
767 AliError(Form("No leading cluster"));
768 return kFALSE;
1c662fe8 769 }
7c97fe2d 770
771
1c662fe8 772 index = emccluster->IdInCollection();
7c97fe2d 773 //add a command to get the index of the leading cluster!
774 if (!emccluster->IsEMCAL()) return kFALSE;
775
bab35745 776 AliVCluster *coi = emccluster->GetCluster();
7c97fe2d 777 if (!coi) return kFALSE;
778
bab35745 779 TLorentzVector vecCOI;
780 coi->GetMomentum(vecCOI,fVertex);
7c97fe2d 781
782
783 Double_t coiTOF = coi->GetTOF()*1e9;
784 if(coiTOF<-30 || coiTOF>30)
785 return kFALSE;
786
787
788 if(ClustTrackMatching(emccluster))
789 return kFALSE;
790
bab35745 791 if(!CheckBoundaries(vecCOI))
7c97fe2d 792 return kFALSE;
793
1c662fe8 794 else
bab35745 795 FillGeneralHistograms(coi,vecCOI, index);
1c662fe8 796 }
797 else{
7c97fe2d 798 //get the entries of the Cluster Container
799 //whatever is a RETURN in LCAnalysis here is a CONTINUE,
800 //since there are more than 1 Cluster per Event
801 // clusters->ResetCurrentID();
802
803 // AliError("fonctionne bien");
804 AliEmcalParticle *emccluster=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(0));
805 index=0;
806 // AliError("fonctionne bien");
807
808 while(emccluster){
809
bab35745 810 AliVCluster *coi = emccluster->GetCluster();
7c97fe2d 811 if(!coi) {
812 index++;
813 emccluster = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
814 continue;
815 }
816 //
bab35745 817 TLorentzVector vecCOI;
818 coi->GetMomentum(vecCOI,fVertex);
7c97fe2d 819 Double_t coiTOF = coi->GetTOF()*1e9;
820 // Double_t coiM02 = coi->GetM02();
821
822 FillQAHistograms(coi,vecCOI);
823 //AliInfo(Form("Cluster number: %d; \t Cluster ToF: %lf ;\tCluster M02:%lf\n",index,coiTOF,coiM02));
824
825 // if(vecCOI.E()<0.3){ // normally already done
826 // index++;
827 // emccluster=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
828 // continue;
829 // }
830
831 if(!fIsMC){
832 if(coiTOF<-30 || coiTOF>30){
833 index++;
834 emccluster=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
835 continue;
836 }
837 }
838 else{//different timing cuts for NON CALIBRATED ToF Signal...
839 if(coiTOF<-570 || coiTOF>630){
840 index++;
841 emccluster=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
842 continue;
843 }
844 }
845 fPtaftTime->Fill(vecCOI.Pt());
846 if(ClustTrackMatching(emccluster)){
847 index++;
848 emccluster = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
1c662fe8 849 continue;
7c97fe2d 850 }
851 fPtaftTM->Fill(vecCOI.Pt());
852
853
854 if(!CheckBoundaries(vecCOI)){
855 index++;
856 emccluster = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
857 continue;
858 }
859
860 fPtaftFC->Fill(vecCOI.Pt());
861
862 if(vecCOI.Et()<5.){
863 index++;
864 emccluster = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
865 continue;
866
867 }
868
869fTestIndexE->Fill(vecCOI.Pt(),index);
870
871 //AliInfo("Passed the CheckBoundaries conditions");
872
873 FillGeneralHistograms(coi,vecCOI, index);
874 index++;
875 emccluster = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(index));
876
1c662fe8 877 }
7c97fe2d 878
1c662fe8 879 }
7c97fe2d 880
881 // PostData(1, fOutput);
1c662fe8 882 return kTRUE;
883}
884
885
7c97fe2d 886 //_________________________________________________________________________________
887void AliAnalysisTaskEMCALPhotonIsolation::FillQAHistograms(AliVCluster *coi, TLorentzVector vecCOI){
888
889 switch(fWho){
890
891 case 0:
892 fevents=0;
893 fEClustersT=vecCOI.E();
894 fPtClustersT=vecCOI.Pt();
895 fEtClustersT=vecCOI.Et();
896 fEtaClustersT=vecCOI.Eta();
897 fPhiClustersT=vecCOI.Phi();
898 fM02ClustersT=coi->GetM02();
899
900 fOutputQATree->Fill();
901
902 break;
903
904 case 1:
905
906 break;
907
908 }
909
910 fPT->Fill(vecCOI.Pt());
911 fE->Fill(vecCOI.E());
912 fM02->Fill(vecCOI.E(),coi->GetM02());
913 fEtaPhiClus->Fill(vecCOI.Eta(),vecCOI.Phi());
914
915 Double_t checktof = coi->GetTOF()*1e9;
916 if(checktof>-30 && checktof<30){
917 fClusTime->Fill(checktof);
918 // fPtaftTime->Fill(vecCOI.Pt());
919
920 // if(!ClustTrackMatching(coi)){
921 // fPtaftTM->Fill(vecCOI.Pt());
922
923 if(CheckBoundaries(vecCOI)){
924 // fPtaftFC->Fill(vecCOI.Pt());
925
926 Double_t checkM02=coi->GetM02();
927 if(fM02mincut < checkM02 && checkM02 < fM02maxcut){
928 fPtaftM02C->Fill(vecCOI.Pt());
bab35745 929 }
1c662fe8 930 }
7c97fe2d 931 // }
932}
933}
934
935
936 //__________________________________________________________________________
937Bool_t AliAnalysisTaskEMCALPhotonIsolation::ClustTrackMatching(AliEmcalParticle *partC) {
938 // Check if the cluster match to a track
939
940
941 AliVCluster *cluster = partC->GetCluster();
942 TLorentzVector nPart;
943 cluster->GetMomentum(nPart, fVertex);
944
945
946AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
947
948AliVCluster *clust = partC -> GetCluster();
949
950Int_t nbMObj = partC -> GetNumberOfMatchedObj();
951
952if (nbMObj == 0) return kFALSE;
953
954for(Int_t i=0;i<nbMObj;i++){
955Int_t imt = partC->GetMatchedObjId(i);
956
957if(imt<0) continue;
958
959AliEmcalParticle *partT = static_cast<AliEmcalParticle*>(tracks->GetAcceptParticle(imt));
960AliVTrack *mt = partT ->GetTrack();
961
962if(!mt) continue;
963
964Double_t deta = 999;
965Double_t dphi = 999;
966
967 Double_t veta = mt->GetTrackEtaOnEMCal();
968 Double_t vphi = mt->GetTrackPhiOnEMCal();
969
970 Float_t pos[3] = {0};
971 clust->GetPosition(pos);
972 TVector3 cpos(pos);
973 Double_t ceta = cpos.Eta();
974 Double_t cphi = cpos.Phi();
975 deta=veta-ceta;
976 dphi=TVector2::Phi_mpi_pi(vphi-cphi);
977
978 fDeltaETAClusTrack->Fill(deta);
979 fDeltaPHIClusTrack->Fill(dphi);
980//
981
982 if(TMath::Abs(dphi)<fdphicut && TMath::Abs(deta)<fdetacut){
983 fDeltaETAClusTrackMatch->Fill(deta);
984 fDeltaPHIClusTrackMatch->Fill(dphi);
985 return kTRUE;
1c662fe8 986 }
7c97fe2d 987
988}
989
990return kFALSE;
1c662fe8 991}
992
993
994
7c97fe2d 995 //__________________________________________________________________________
bab35745 996void AliAnalysisTaskEMCALPhotonIsolation::EtIsoCellPhiBand(TLorentzVector c, Float_t &etIso, Float_t &phiBandcells){
7c97fe2d 997 // Underlying events study with EMCAL cells in phi band // have to be tested
998
bab35745 999 AliEMCALGeometry* emcalGeom = AliEMCALGeometry::GetInstance();
7c97fe2d 1000
1c662fe8 1001 Float_t sumEnergyPhiBandCells=0., sumEnergyConeCells=0.;
7c97fe2d 1002
1003
1004 // check the cell corresponding to the leading cluster
1c662fe8 1005 Int_t absId = 999;
7c97fe2d 1006 //maybe best to call it LeadingCellIdinClus or better maxId since it won't be used anywhere else ???
bab35745 1007 Bool_t cellLeadingClustID = emcalGeom->GetAbsCellIdFromEtaPhi(c.Eta(),c.Phi(),absId);
1008 if(!cellLeadingClustID) return;
7c97fe2d 1009
1c662fe8 1010 else{
1011 Int_t iTower = -1;
1012 Int_t iModule = -1;
1013 Int_t imEta=-1, imPhi=-1;
bab35745 1014 Int_t iEta =-1, iPhi =-1;
7c97fe2d 1015
bab35745 1016 emcalGeom->GetCellIndex(absId,iModule,iTower,imPhi,imEta); // to get the module, the tower, eta and phi for the cell corresponding to the leading cluster
1017 emcalGeom->GetCellPhiEtaIndexInSModule(iModule,iTower,imPhi,imEta,iPhi,iEta); // to get the cell eta and phi in the super module for the cell coming from the leading cluster
7c97fe2d 1018
1019 // Get the row and the column of the cell corresponding to the leading cluster in EMCAL
bab35745 1020 Int_t colCellLeadingClust = iEta;
1021 if(iModule % 2) colCellLeadingClust = AliEMCALGeoParams::fgkEMCALCols + iEta ; // if the SM number is even you need to translate to have the one corresponding in EMCAL
1022 Int_t rowCellLeadingClust = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(iModule/2); // to have the corresponding row in EMCAL
7c97fe2d 1023
1024 // total number or rows and columns in EMCAL
bab35745 1025 Int_t nTotalRows = AliEMCALGeoParams::fgkEMCALRows*16/3 ; // 5 + 2/3 supermodules in a row
1026 Int_t nTotalCols = 2*AliEMCALGeoParams::fgkEMCALCols; // 2 supermodules in a column
7c97fe2d 1027
bab35745 1028 Int_t nbConeSize = int(fIsoConeRadius/0.0143); //0.0143 tower acceptance
7c97fe2d 1029
1030 // Get the cells
1c662fe8 1031 AliVCaloCells * cells =InputEvent()->GetEMCALCells();
7c97fe2d 1032
1033 // define the max and min row and column corresponding to the isolation cone around the seed cell from the leading cluster
bab35745 1034 Int_t iRowMinCone = rowCellLeadingClust-nbConeSize;
1035 if(iRowMinCone<0) iRowMinCone=0;
7c97fe2d 1036
bab35745 1037 Int_t iRowMaxCone = rowCellLeadingClust+nbConeSize;
1038 if(iRowMaxCone>AliEMCALGeoParams::fgkEMCALRows) iRowMaxCone=AliEMCALGeoParams::fgkEMCALRows; // AliEMCALGeoParams::fgkEMCALRows = 24 in a supermodule
7c97fe2d 1039
bab35745 1040 Int_t iColMinCone = colCellLeadingClust - nbConeSize;
1041 if(iColMinCone<0) iColMinCone = 0;
7c97fe2d 1042
bab35745 1043 Int_t iColMaxCone = colCellLeadingClust+nbConeSize;
1044 if(iColMaxCone>AliEMCALGeoParams::fgkEMCALCols) iColMaxCone=AliEMCALGeoParams::fgkEMCALCols; // AliEMCALGeoParams::fgkEMCALCols = 48 in a supermodule
7c97fe2d 1045
1046 // loop on all cells
bab35745 1047 for(Int_t iCol=0; iCol<nTotalCols; iCol++){
1048 for(Int_t iRow=0; iRow<nTotalRows; iRow++){
7c97fe2d 1049 // now recover the cell indexes in a supermodule
bab35745 1050 Int_t iSector = int(iRow/AliEMCALGeoParams::fgkEMCALRows); // check in which SM is the cell
1c662fe8 1051 if(iSector==5) continue; //
1052 Int_t inModule = -1;
bab35745 1053 Int_t iColLoc = -1;
1054 if(iCol < AliEMCALGeoParams::fgkEMCALCols){ // if the SM number is odd the column is the one corresponding in the supermodule
494479c2 1055 inModule = 2*iSector + 1;
bab35745 1056 iColLoc = iCol;
1c662fe8 1057 }
bab35745 1058 else if(iCol > AliEMCALGeoParams::fgkEMCALCols - 1){ // if the SM number is even the column isn't the one corresponding in the supermodule
1c662fe8 1059 inModule = 2*iSector;
bab35745 1060 iColLoc = iCol-AliEMCALGeoParams::fgkEMCALCols;
1c662fe8 1061 }
7c97fe2d 1062
bab35745 1063 Int_t iRowLoc = iRow - AliEMCALGeoParams::fgkEMCALRows*iSector ;
7c97fe2d 1064
bab35745 1065 if(TMath::Abs(iCol-colCellLeadingClust)<nbConeSize && TMath::Abs(iCol+colCellLeadingClust)>nbConeSize){
1066 if(iRow<iRowMaxCone && iRow>iRowMinCone){
7c97fe2d 1067 Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(iModule,iRow,iCol); // verifier les iRow et iCol
1c662fe8 1068 sumEnergyPhiBandCells+=cells->GetAmplitude(iabsId); //should be Et
1069 }
1070 }
bab35745 1071 else if (TMath::Abs(iCol-colCellLeadingClust)>nbConeSize && TMath::Abs(iCol+colCellLeadingClust)<nbConeSize){
494479c2 1072 Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(inModule,iRowLoc,iColLoc); // verifier les iRow et iCol
1c662fe8 1073 sumEnergyConeCells+=cells->GetAmplitude(iabsId); //should be Et
1074 }
1075 }
1076 }
1077 }
bab35745 1078 etIso = sumEnergyConeCells;
1079 phiBandcells = sumEnergyPhiBandCells;
1c662fe8 1080}
1081
1082
7c97fe2d 1083 //__________________________________________________________________________
bab35745 1084void AliAnalysisTaskEMCALPhotonIsolation::EtIsoCellEtaBand(TLorentzVector c, Float_t &etIso, Float_t &etaBandcells){
7c97fe2d 1085 // Underlying events study with EMCAL cell in eta band // have to be tested
1086
1087
bab35745 1088 AliEMCALGeometry* emcalGeom = AliEMCALGeometry::GetInstance();
7c97fe2d 1089
1c662fe8 1090 Float_t sumEnergyEtaBandCells=0., sumEnergyConeCells=0.;
7c97fe2d 1091
1092
1093
1094 // check the cell corresponding to the leading cluster
1c662fe8 1095 Int_t absId = 999;
7c97fe2d 1096 //maybe best to call it LeadingCellIdinClus or better maxId since it won't be used anywhere else ???
bab35745 1097 Bool_t cellLeadingClustID = emcalGeom->GetAbsCellIdFromEtaPhi(c.Eta(),c.Phi(),absId);
1098 if(!cellLeadingClustID) return;
7c97fe2d 1099
1c662fe8 1100 else{
7c97fe2d 1101
1c662fe8 1102 Int_t iTower = -1;
1103 Int_t iModule = -1;
1104 Int_t imEta=-1, imPhi=-1;
bab35745 1105 Int_t iEta =-1, iPhi =-1;
7c97fe2d 1106
bab35745 1107 emcalGeom->GetCellIndex(absId,iModule,iTower,imPhi,imEta); // to get the module, the tower, eta and phi for the cell corresponding to the leading cluster
1108 emcalGeom->GetCellPhiEtaIndexInSModule(iModule,iTower,imPhi,imEta,iPhi,iEta); // to get the cell eta and phi in the super module for the cell coming from the leading cluster
7c97fe2d 1109
1110 // Get the row and the column of the cell corresponding to the leading cluster in EMCAL
bab35745 1111 Int_t colCellLeadingClust = iEta;
1112 if(iModule % 2) colCellLeadingClust = AliEMCALGeoParams::fgkEMCALCols + iEta ; // if the SM number is even you need to translate to have the one corresponding in EMCAL
1113 Int_t rowCellLeadingClust = iPhi + AliEMCALGeoParams::fgkEMCALRows*int(iModule/2); // to have the corresponding row in EMCAL
7c97fe2d 1114
1115 // total number or rows and columns in EMCAL
bab35745 1116 Int_t nTotalRows = AliEMCALGeoParams::fgkEMCALRows*16/3 ; // 5 + 2/3 supermodules in a row
1117 Int_t nTotalCols = 2*AliEMCALGeoParams::fgkEMCALCols; // 2 supermodules in a column
7c97fe2d 1118
bab35745 1119 Int_t nbConeSize = int(fIsoConeRadius/0.0143); //0.0143 tower acceptance
7c97fe2d 1120
1121 // Get the cells
1c662fe8 1122 AliVCaloCells * cells =InputEvent()->GetEMCALCells();
7c97fe2d 1123
1124 // define the max and min row and column corresponding to the isolation cone around the seed cell from the leading cluster
bab35745 1125 Int_t iRowMinCone = rowCellLeadingClust-nbConeSize;
1126 if(iRowMinCone<0) iRowMinCone=0;
7c97fe2d 1127
bab35745 1128 Int_t iRowMaxCone = rowCellLeadingClust+nbConeSize;
1129 if(iRowMaxCone>AliEMCALGeoParams::fgkEMCALRows) iRowMaxCone=AliEMCALGeoParams::fgkEMCALRows; // AliEMCALGeoParams::fgkEMCALRows = 24 in a supermodule
7c97fe2d 1130
bab35745 1131 Int_t iColMinCone = colCellLeadingClust-nbConeSize;
1132 if(iColMinCone<0) iColMinCone = 0;
7c97fe2d 1133
bab35745 1134 Int_t iColMaxCone = colCellLeadingClust+nbConeSize;
1135 if(iColMaxCone>AliEMCALGeoParams::fgkEMCALCols) iColMaxCone=AliEMCALGeoParams::fgkEMCALCols; // AliEMCALGeoParams::fgkEMCALCols = 48 in a supermodule
7c97fe2d 1136
1137 // loop on all cells
bab35745 1138 for(Int_t iCol=0; iCol<nTotalCols; iCol++)
1c662fe8 1139 {
bab35745 1140 for(Int_t iRow=0; iRow<nTotalRows; iRow++)
1c662fe8 1141 {
7c97fe2d 1142 // now recover the cell indexes in a supermodule
bab35745 1143 Int_t iSector = int(iRow/AliEMCALGeoParams::fgkEMCALRows); // check in which SM is the cell
1c662fe8 1144 if(iSector==5) continue; //
1145 Int_t inModule = -1;
bab35745 1146 Int_t iColLoc = -1;
1147 if(iCol < AliEMCALGeoParams::fgkEMCALCols){ // if the SM number is odd the column is the one corresponding in the supermodule
494479c2 1148 inModule = 2*iSector + 1;
bab35745 1149 iColLoc = iCol;
1c662fe8 1150 }
bab35745 1151 else if(iCol > AliEMCALGeoParams::fgkEMCALCols - 1){ // if the SM number is even the column isn't the one corresponding in the supermodule
1c662fe8 1152 inModule = 2*iSector;
bab35745 1153 iColLoc = iCol-AliEMCALGeoParams::fgkEMCALCols;
1c662fe8 1154 }
7c97fe2d 1155
bab35745 1156 Int_t iRowLoc = iRow - AliEMCALGeoParams::fgkEMCALRows*iSector ;
7c97fe2d 1157
bab35745 1158 if(TMath::Abs(iCol-colCellLeadingClust)<nbConeSize && TMath::Abs(iCol+colCellLeadingClust)>nbConeSize){
1159 if(iCol<iColMaxCone && iCol>iColMinCone){
7c97fe2d 1160 Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(iModule,iRow,iCol); // verifier les iRow et iCol
1c662fe8 1161 sumEnergyEtaBandCells+=cells->GetAmplitude(iabsId); //should be Et
1162 }
1163 }
bab35745 1164 else if (TMath::Abs(iCol-colCellLeadingClust)>nbConeSize && TMath::Abs(iCol+colCellLeadingClust)<nbConeSize){
494479c2 1165 Int_t iabsId = emcalGeom->GetAbsCellIdFromCellIndexes(inModule,iRowLoc,iColLoc); // verifier les iRow et iCol
1c662fe8 1166 sumEnergyConeCells+=cells->GetAmplitude(iabsId); //should be Et
1167 }
1168 }
1169 }
1170 }
bab35745 1171 etIso = sumEnergyConeCells;
1172 etaBandcells = sumEnergyEtaBandCells;
1c662fe8 1173}
1174
1175
7c97fe2d 1176 //__________________________________________________________________________
1177void AliAnalysisTaskEMCALPhotonIsolation::EtIsoClusPhiBand(TLorentzVector c, Float_t &ptIso, Float_t &phiBandclus, Int_t index){
1178 // Underlying events study with clusters in phi band
1179
1180 Float_t sumEnergyPhiBandClus=0., sumEnergyConeClus=0., sumpTConeCharged=0.;
1181
1182 //needs a check on the same cluster
1c662fe8 1183 AliParticleContainer *clusters = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
7c97fe2d 1184 clusters->ResetCurrentID();
1185 Int_t localIndex=0;
1186 AliEmcalParticle *clust = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1187
1188 while(clust){ //check the position of other clusters in respect to the leading cluster
1189
1190 if(localIndex==index){
1191 localIndex++;
1192 clust = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1193 continue;
1194 }
1195 else{
1196 localIndex++;
1197
1198 AliVCluster *cluster= clust->GetCluster();
1199 if(!cluster){
1200 clust = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1201 continue;
1202 }
1203
1204 TLorentzVector nClust; //STILL NOT INITIALIZED
1205 cluster->GetMomentum(nClust,fVertex);
1206 Float_t phiClust =nClust.Phi();
1207 Float_t etaClust= nClust.Eta();
1208
1209 Double_t clustTOF = cluster->GetTOF()*1e9;
1210 if(clustTOF<-30 || clustTOF>30){
1211 clust=static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1212 continue;
1213 }
1214
1215 if(ClustTrackMatching(clust)){
1216 clust=static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1217 continue;
1218 }
1219
1220 //redefine phi/c.Eta() from the cluster we passed to the function
1221
1222 Float_t radius = TMath::Sqrt(TMath::Power(phiClust- c.Phi(),2)+TMath::Power(etaClust-c.Eta(),2)); // define the radius between the leading cluster and the considered cluster
1223
1224 if(radius>fIsoConeRadius){ // the cluster is outside the isolation cone in this case study UE
1225
1226 // actually phi band here
1227 if(TMath::Abs(etaClust - c.Eta()) < fIsoConeRadius){
1228 sumEnergyPhiBandClus += nClust.Pt();
1229 }
1230 }
1231 else // if the cluster is in the isolation cone, add the cluster pT
1232 sumEnergyConeClus += nClust.Et();
1233
1234 clust = static_cast<AliEmcalParticle*>(clusters->GetNextAcceptParticle(localIndex));
1235 }
1236 }
1237
1238 fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1239 // name hard coded to use the defined tracks for analysis
1240
1241 if (!fTracksAna) {
1242 AliError(Form("Could not retrieve tracks !"));
1243 return;
1244 }
1245 const Int_t nbTracks = fTracksAna->GetEntries();
1246 Int_t iTracks = 0;
1247
1248//
1249 while(iTracks<nbTracks){
1250 AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1251 if(!track){
1252 AliError(Form("No tracks in collection"));
1253 iTracks++;
1254 continue;
1255 }
1256 //CHECK IF TRACK IS IN BOUNDARIES
1257 Float_t phiTrack = track->Phi();
1258 Float_t etaTrack = track->Eta();
1259
1260 Float_t radius = TMath::Sqrt(TMath::Power(phiTrack - c.Phi(),2)+TMath::Power(etaTrack - c.Eta(),2));
1261
1262 if(radius<fIsoConeRadius){ // if tracks are outside the isolation cone study
1263 sumpTConeCharged+=track->Pt(); // should not double count if the track matching is already done
1264 }
1265 iTracks++;
1266 } // end of tracks loop
1267
1268 fTestEnergyCone->Fill(sumEnergyConeClus,sumpTConeCharged);
1269
1270
1271 ptIso = sumEnergyConeClus + sumpTConeCharged;
bab35745 1272 phiBandclus = sumEnergyPhiBandClus;
1c662fe8 1273}
1274
1275
7c97fe2d 1276 //__________________________________________________________________________
1277void AliAnalysisTaskEMCALPhotonIsolation::EtIsoClusEtaBand(TLorentzVector c, Float_t &ptIso, Float_t &etaBandclus, Int_t index){
1278 // Underlying events study with clusters in eta band
1279
1280
1281
1282 Float_t sumEnergyEtaBandClus =0., sumEnergyConeClus=0., sumpTConeCharged=0;
1283 Double_t clustTOF=0;
1284
1c662fe8 1285 AliParticleContainer *clusters = static_cast<AliParticleContainer*>(fParticleCollArray.At(1));
7c97fe2d 1286
1287// clusters->ResetCurrentID();
1288 Int_t localIndex=0;
1289 AliEmcalParticle *clust = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1290
1291 while(clust){ //check the position of other clusters in respect to the leading cluster
1292
1293 if(localIndex==index){
1294 localIndex++;
1295 clust = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1296 continue;
1297 }
1298
1299 else{
1300 localIndex++;
1301 AliVCluster *cluster= clust->GetCluster();
1302 if(!cluster){
1303 clust = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1304 continue;
1c662fe8 1305 }
7c97fe2d 1306
1307 TLorentzVector nClust; //STILL NOT INITIALIZED
1308 cluster->GetMomentum(nClust,fVertex);
1309
1310 Float_t phiClust =nClust.Phi();
1311 Float_t etaClust= nClust.Eta();
1312 Float_t eTcluster=0;
1313
1314
1315 clustTOF = cluster->GetTOF()*1e9;
1316 if(clustTOF<-30 || clustTOF>30){
1317 clust=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1318 continue;
1319 }
1320
1321 if(ClustTrackMatching(clust)){
1322 clust=static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1323 continue;
1324 }
1325 //redefine phi/c.Eta() from the cluster we passed to the function
1326
1327 // define the radius between the leading cluster and the considered cluster
1328 Float_t radius = TMath::Sqrt(TMath::Power(phiClust-c.Phi(),2)+TMath::Power(etaClust-c.Eta(),2));
1329
1330
1331
1332 if(radius>fIsoConeRadius){ // the cluster is outside the isolation cone in this case study UE
1333
1334 // actually eta band here
1335 if(TMath::Abs(etaClust - c.Eta()) < fIsoConeRadius){
1336 sumEnergyEtaBandClus += nClust.Et();
1337 }
1338 }
1339 else if(radius<fIsoConeRadius && radius!=0.){ // if the cluster is in the isolation cone, add the cluster pT
1340 eTcluster=nClust.Pt();
1341
1342
1343 sumEnergyConeClus += nClust.Pt();
1344 fTestEtaPhiCone->Fill(c.Eta(),c.Phi());
1345 fTestIndex->Fill(index,localIndex);
1346
1347 fTestLocalIndexE->Fill(eTcluster,localIndex);
1c662fe8 1348 }
7c97fe2d 1349 clust = static_cast<AliEmcalParticle*>(clusters->GetAcceptParticle(localIndex));
1350 }
1351 } // end of clusters loop
1352
1353
1354
1355 fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1356 // name hard coded to use the defined tracks for analysis
1357
1358 if (!fTracksAna) {
1359 AliError(Form("Could not retrieve tracks !"));
1360 return;
1361 }
1362 const Int_t nbTracks = fTracksAna->GetEntries();
1363 Int_t iTracks = 0;
1364
1365//
1366 while(iTracks<nbTracks){
1367 AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1368 if(!track){
1369 AliError(Form("No tracks in collection"));
1370 iTracks++;
1371 continue;
1372 }
1373 //CHECK IF TRACK IS IN BOUNDARIES
1374 Float_t phiTrack = track->Phi();
1375 Float_t etaTrack = track->Eta();
1376
1377 Float_t radius = TMath::Sqrt(TMath::Power(phiTrack - c.Phi(),2)+TMath::Power(etaTrack - c.Eta(),2));
1378
1379 if(radius<fIsoConeRadius){ // if tracks are outside the isolation cone study
1380 sumpTConeCharged+=track->Pt(); // should not double count if the track matching is already done
1381 }
1382 iTracks++;
1383 } // end of tracks loop
1384
1385 fTestEnergyCone->Fill(sumEnergyConeClus,sumpTConeCharged);
1386
1387 ptIso = sumEnergyConeClus + sumpTConeCharged;
bab35745 1388 etaBandclus = sumEnergyEtaBandClus;
1c662fe8 1389}
1390
1391
7c97fe2d 1392
1393 //__________________________________________________________________________
bab35745 1394void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackPhiBand(TLorentzVector c, Float_t &ptIso, Float_t &phiBandtrack){
7c97fe2d 1395 // Underlying events study with tracks in phi band in EMCAL acceptance
1396
1397 //INSERT BOUNDARIES ACCORDING TO THE FLAG WE WANT TO USE
1c662fe8 1398 Float_t sumpTConeCharged=0., sumpTPhiBandTrack=0.;
bab35745 1399 Float_t minPhi= 0., maxPhi= 2*TMath::Pi(), minEta = -0.9, maxEta= 0.9;
7c97fe2d 1400
1401 if(!fTPC4Iso)
1402 {
bab35745 1403 minEta = -0.7;
1404 maxEta = 0.7;
1405 minPhi = 1.4;
1406 maxPhi = TMath::Pi();
7c97fe2d 1407 }
1408
1409
1410 fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1411
1412 if (!fTracksAna)
1413 {
1414 AliError(Form("Could not retrieve tracks !"));
1415 return;
1416 }
1417 const Int_t nbTracks = fTracksAna->GetEntries();
1418 Int_t iTracks = 0;
1419
1420
1421 while(iTracks<nbTracks)
1422 {
1423 AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1424 if(!track)
1425 {
1426 AliError(Form("No tracks in collection"));
1427 iTracks++;
1428 continue;
1429 }
1430 //CHECK IF TRACK IS IN BOUNDARIES
1431 Float_t phiTrack = track->Phi();
1432 Float_t etaTrack = track->Eta();
1433 // define the radius between the leading cluster and the considered cluster
1434 //redefine phi/c.Eta() from the cluster we passed to the function
1435 if(phiTrack < maxPhi && phiTrack > minPhi && etaTrack < maxEta && etaTrack > minEta)
1436 {
1437 Float_t radius = TMath::Sqrt(TMath::Power(phiTrack - c.Phi(),2)+TMath::Power(etaTrack - c.Eta(),2));
1438
1439 if(radius>fIsoConeRadius)
1440 { // if tracks are outside the isolation cone study
1441
1442 // actually phi band here --- ADD Boundaries conditions
1443 if(TMath::Abs(etaTrack - c.Eta()) < fIsoConeRadius)
1444 {
1445 sumpTPhiBandTrack += track->Pt();
1446 }
1447 }
1448 else
1449 sumpTConeCharged+=track->Pt(); // should not double count if the track matching is already done
1c662fe8 1450 }
7c97fe2d 1451 iTracks++;
1c662fe8 1452 }
bab35745 1453 ptIso = sumpTConeCharged;
1454 phiBandtrack = sumpTPhiBandTrack;
1c662fe8 1455}
1456
1457
7c97fe2d 1458 //__________________________________________________________________________
bab35745 1459void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackEtaBand(TLorentzVector c, Float_t &ptIso, Float_t &etaBandtrack){
7c97fe2d 1460 // Underlying events study with tracks in eta band in EMCAL acceptance
1461
1462 //INSERT BOUNDARIES ACCORDING TO THE FLAG WE WANT TO USE
1c662fe8 1463 Float_t sumpTConeCharged=0., sumpTEtaBandTrack=0.;
bab35745 1464 Float_t minPhi= 0., maxPhi= 2*TMath::Pi(), minEta = -0.9, maxEta= 0.9;
7c97fe2d 1465
1c662fe8 1466 if(!fTPC4Iso){
bab35745 1467 minEta = -0.7;
1468 maxEta = 0.7;
1469 minPhi = 1.4;
1470 maxPhi = TMath::Pi();
1c662fe8 1471 }
7c97fe2d 1472
1473 fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1474
1475 if (!fTracksAna)
1476 {
1477 AliError(Form("Could not retrieve tracks !"));
1478 return;
1479 }
1480 const Int_t nbTracks = fTracksAna->GetEntries();
1481 Int_t iTracks = 0;
1482
1483
1484 while(iTracks<nbTracks)
1485 {
1486 AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1487
1488 if(!track)
1489 {
1490 AliError(Form("No tracks in collection"));
1491 iTracks++;
1492 continue;
1493 }
1494
bab35745 1495 Float_t phiTrack = track->Phi();
1496 Float_t etaTrack = track->Eta();
7c97fe2d 1497 //redefine phi/c.Eta() from the cluster we passed to the function
1498 if(phiTrack < maxPhi && phiTrack > minPhi && etaTrack < maxEta && etaTrack > minEta)
1499 {
1500 Float_t radius = TMath::Sqrt(TMath::Power(phiTrack-c.Phi(),2)+TMath::Power(etaTrack-c.Eta(),2)); // define the radius between the leading cluster and the considered cluster
1501
1502 if(radius>fIsoConeRadius)
1503 { // if tracks are outside the isolation cone study UE
1504
1505 // actually eta band here --- ADD Boundaries conditions
1506 if(TMath::Abs(phiTrack - c.Phi()) < fIsoConeRadius)
1507 {
1508 sumpTEtaBandTrack += track->Pt();
1509 }
1510 }
1511 else sumpTConeCharged += track->Pt(); // should not double count if the track matching is already done
1c662fe8 1512 }
7c97fe2d 1513 iTracks++;
1c662fe8 1514 }
7c97fe2d 1515
bab35745 1516 ptIso = sumpTConeCharged;
1517 etaBandtrack = sumpTEtaBandTrack;
1c662fe8 1518}
1519
1520
7c97fe2d 1521 //__________________________________________________________________________
bab35745 1522void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackOrthCones(TLorentzVector c, Float_t &ptIso, Float_t &cones){
7c97fe2d 1523 // Underlying events study with tracks in orthogonal cones in TPC
1524
1c662fe8 1525 Float_t sumpTConeCharged=0., sumpTPerpConeTrack=0.;
7c97fe2d 1526
bab35745 1527 Float_t etaClus = c.Eta();
1528 Float_t phiClus = c.Phi();
1529 Float_t phiCone1 = phiClus - TMath::PiOver2();
1530 Float_t phiCone2 = phiClus + TMath::PiOver2();
7c97fe2d 1531
bab35745 1532 if (phiCone1 < 0.) phiCone1 += 2*TMath::Pi();
7c97fe2d 1533
1534 fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1535
1536 if (!fTracksAna)
1537 {
1538 AliError(Form("Could not retrieve tracks !"));
1539 return;
1540 }
1541
1542 const Int_t nbTracks = fTracksAna->GetEntries();
1543 Int_t iTracks = 0;
1544
1545 while(iTracks<nbTracks)
1546 {
1547
1548 AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1549
1550 if(!track)
1551 {
1552 AliError(Form("No tracks in collection"));
1553 iTracks++;
1554 continue;
1555 }
1556
bab35745 1557 Float_t phiTrack = track->Phi();
1558 Float_t etaTrack = track->Eta();
bab35745 1559 Float_t dist2Clust = TMath::Sqrt(TMath::Power(etaTrack-etaClus, 2)+TMath::Power(phiTrack-phiClus, 2));
7c97fe2d 1560
1561 if (dist2Clust<fIsoConeRadius) sumpTConeCharged += track->Pt(); // tracks are inside the isolation cone
1562
1563 else
1564 {//tracks outside the IsoCone
1565 //Distances from the centres of the two Orthogonal Cones
bab35745 1566 Float_t dist2Cone1 = TMath::Sqrt(TMath::Power(etaTrack-etaClus, 2)+TMath::Power(phiTrack-phiCone1, 2));
1567 Float_t dist2Cone2 = TMath::Sqrt(TMath::Power(etaTrack-etaClus, 2)+TMath::Power(phiTrack-phiCone2, 2));
7c97fe2d 1568
1569 //Is the Track Inside one of the two Cones ->Add to UE
1570 if((dist2Cone1 < fIsoConeRadius) || (dist2Cone2 < fIsoConeRadius)) sumpTPerpConeTrack += track->Pt();
1571
1c662fe8 1572 }
7c97fe2d 1573
1574 iTracks++;
1575
1c662fe8 1576 }
7c97fe2d 1577
bab35745 1578 ptIso = sumpTConeCharged;
1579 cones = sumpTPerpConeTrack;
1c662fe8 1580}
1581
7c97fe2d 1582 //__________________________________________________________________________
bab35745 1583void AliAnalysisTaskEMCALPhotonIsolation::PtIsoTrackFullTPC(TLorentzVector c, Float_t &ptIso, Float_t &full){
7c97fe2d 1584 // Underlying events study with tracks in full TPC except back to back bands
1585
1c662fe8 1586 Float_t sumpTConeCharged=0., sumpTTPCexceptB2B=0.;
7c97fe2d 1587
1588 fTracksAna = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("ESDFilterTracksAna"));
1589
1590 if (!fTracksAna)
1591 {
1592 AliError(Form("Could not retrieve tracks !"));
1593 return;
1594 }
1595
1596 const Int_t nbTracks = fTracksAna->GetEntries();
1597 Int_t iTracks = 0;
1598
1599 while(iTracks<nbTracks)
1600 {
1601
1602 AliVTrack *track = static_cast<AliVTrack*>(fTracksAna->At(iTracks));
1603
1604 if(!track)
1605 {
1606 AliError(Form("No tracks in collection"));
1607 iTracks++;
1608 continue;
1609 }
1610
bab35745 1611 Float_t phiTrack = track->Phi();
1612 Float_t etaTrack = track->Eta();
7c97fe2d 1613 //redefine phi/c.Eta() from the cluster we passed to the function
bab35745 1614 Float_t radius = TMath::Sqrt(TMath::Power(phiTrack-c.Phi(),2)+TMath::Power(etaTrack-c.Eta(),2)); // define the radius between the leading cluster and the considered cluster
7c97fe2d 1615
1616 if(radius>fIsoConeRadius)
1617 { // if tracks are outside the isolation cone study UE
bab35745 1618 Double_t dphiUp = c.Phi() + TMath::Pi() - fIsoConeRadius;
1619 Double_t dphiDown = c.Phi() + TMath::Pi() + fIsoConeRadius;
7c97fe2d 1620 // TPC except B2B
1621 if(phiTrack < dphiDown && phiTrack> dphiUp) sumpTTPCexceptB2B += track->Pt();
1622
1c662fe8 1623 }
7c97fe2d 1624
1625 else sumpTConeCharged += track->Pt(); // should not double count if the track matching is already done
1626
1627 iTracks++;
1628
1c662fe8 1629 }
7c97fe2d 1630
bab35745 1631 ptIso = sumpTConeCharged;
1632 full = sumpTTPCexceptB2B;
1c662fe8 1633}
1634
7c97fe2d 1635 //__________________________________________________________________________
bab35745 1636Bool_t AliAnalysisTaskEMCALPhotonIsolation::CheckBoundaries(TLorentzVector vecCOI){
7c97fe2d 1637 // Check if the cone around the considered cluster is in EMCAL acceptance
1638 //AliInfo("Inside CheckBoundaries\n");
1639
bab35745 1640 Float_t minPhiBound= 1.5 , minEtaBound= -0.5, maxPhiBound= TMath::Pi()-0.1, maxEtaBound= 0.5;
1c662fe8 1641 Bool_t isINBoundaries;
7c97fe2d 1642
1643if(!fTPC4Iso)
1644 {
1645 minPhiBound = 1.4+fIsoConeRadius;
1646 maxPhiBound = 3.15-fIsoConeRadius; // normally 110° but shorter cut to avoid EMCAL border
1647 minEtaBound = -0.67+fIsoConeRadius; // ""
1648 maxEtaBound = 0.67-fIsoConeRadius; // ""
1649
1650 // minPhiBound = 1.8; //to be changed with fIsoConeR
1651 // maxPhiBound = 2.75;
1652 // minEtaBound = -0.27;
1653 // maxEtaBound = 0.27;
1654}
1655
1656
1657 if(vecCOI.Eta() > maxEtaBound || vecCOI.Eta() < minEtaBound || vecCOI.Phi() > maxPhiBound || vecCOI.Phi() <minPhiBound)
1c662fe8 1658 isINBoundaries=kFALSE;
1659 else
1660 isINBoundaries=kTRUE;
7c97fe2d 1661
1662
1c662fe8 1663 return isINBoundaries;
1664}
1665
7c97fe2d 1666 //_________________________________________________________________________
1667void AliAnalysisTaskEMCALPhotonIsolation::LookforParticle(Int_t clusterlabel, Double_t energyCLS, Double_t phiCLS, Double_t etaCLS, Double_t time, Double_t ss, Int_t ncells){
1668
1669 //time and ncells not used for the moment
1670
1671 // AliInfo("Inside AnalyzeMC");
1672 if (!fIsMC)
1673 {
1674 cout<<"not a montecarlo run!!!!!!"<<endl;
1675 return;
1676 } //AliInfo(Form("It's a MC analysis %e",fAODMCParticles));
1677 if(!fStack && !fAODMCParticles){
1678 cout<<"No Particle Stack !!!!!"<<endl;
1679 return;
1680 }
1681 //AliInfo("there's a List of particles");
1682 //DO THIS ALSO FOR ESDs
1683
1684 if(fAODMCParticles->GetEntries() < 1){
1685 cout<<"number of tracks insufficient"<<endl;
1686 return;
1687 }
1688
1689
1690 Int_t ndimsMCmix = fMCQAdim;
1691
1692
1693 Double_t outputvalueMCmix[ndimsMCmix];
1694 //cout<<"dimensions of the array: "<<ndimsMCmix<<endl;
1695
1696
1697 Int_t npart=fAODMCParticles->GetEntries();
1698 //cout<<"Number of particles in the event: "<<npart<<endl;
1699
1700 AliAODMCParticle *particle2Check, *MomP2Check;
1701 // *partMC;
1702
1703 Int_t clustPDG, p2clabel;
1704 Double_t enTrue,phiTrue, etaTrue;
1705 Double_t dPhi,dEta ;
1706 bool found=kFALSE;
1707 for(int b=0; b<npart && found!= kTRUE ;b++){
1708 particle2Check = static_cast<AliAODMCParticle*>(fAODMCParticles->At(b));
1709 p2clabel = particle2Check->Label();
1710
1711 if(clusterlabel==p2clabel){
1712 found=kTRUE;
1713 clustPDG = particle2Check->GetPdgCode();
1714 int mom2checkidx = particle2Check->GetMother();
1715 MomP2Check = static_cast<AliAODMCParticle*>(fAODMCParticles->At(mom2checkidx));
1716 //if(energyCLS>=40.)
1717 //cout<<"PDG associated: "<<clustPDG<<" Mother PDG: "<<MomP2Check->GetPdgCode()<<endl;
1718 if(clustPDG==22 || (TMath::Abs(clustPDG)==11 && MomP2Check->GetPdgCode()==22)) //continue;
1719 {
1720 phiTrue = particle2Check->Phi();
1721 etaTrue = particle2Check->Eta();
1722 enTrue = particle2Check->E()*TMath::Sin(particle2Check->Theta());
1723 //if(energyCLS>=40.)
1724 //cout<<"Energy of the single particle associated with the cluster: "<<enTrue<<endl;
1725 if(clustPDG==22){
1726 if(MomP2Check->GetPdgCode()==111){
1727
1728 Int_t idxdaug1 = MomP2Check->GetFirstDaughter();
1729 if (idxdaug1<npart){
1730 if(idxdaug1==clusterlabel){
1731 Int_t idxdaug2 = MomP2Check->GetLastDaughter();
1732 if(idxdaug2<npart){
1733 AliAODMCParticle *daug2 = static_cast<AliAODMCParticle*>(fAODMCParticles->At(idxdaug2));
1734 if(daug2->GetPdgCode()==22 && (daug2->Phi()-phiTrue)<0.2 && (daug2->Eta()-etaTrue)<0.2){
1735 //if(energyCLS >= 40.){
1736 //cout<<"CASE1\nPDG of the other particle VERY close: "<<daug2->GetPdgCode()<<" with Label: "<<daug2->Label();
1737 //cout<<" Energy of the other particle VERY close: "<<daug2->E()*TMath::Sin(daug2->Theta())<<endl;
1738 //}
1739 enTrue += daug2->E()*TMath::Sin(daug2->Theta());
1740 }
1741 }
1742 }
1743 else{
1744 AliAODMCParticle *daug1 = static_cast<AliAODMCParticle*>(fAODMCParticles->At(idxdaug1));
1745
1746 if(daug1->GetPdgCode()==22 && (daug1->Phi()-phiTrue)<0.2 && (daug1->Eta()-etaTrue)<0.2){
1747 //if(energyCLS >= 40.){
1748 //cout<<"CASE2\nPDG of the other particle VERY close: "<<daug1->GetPdgCode()<<" with Label: "<<daug1->Label();
1749 //cout<<" Energy of the other particle VERY close: "<<daug1->E()*TMath::Sin(daug1->Theta())<<endl;
1750 //}
1751 enTrue += daug1->E()*TMath::Sin(daug1->Theta());
1752 }
1753 }
1754 }
1755 }
1756 }
1757 else{
1758 Int_t firstidx=MomP2Check->GetFirstDaughter();
1759 if(firstidx< npart){
1760 if(firstidx==clusterlabel){
1761 Int_t lastidx=MomP2Check->GetLastDaughter();
1762 if(lastidx<npart){
1763 AliAODMCParticle *last=static_cast<AliAODMCParticle*>(fAODMCParticles->At(lastidx));
1764 if((last->Phi()-phiTrue)<0.03 && (last->Eta()-etaTrue)<0.02){
1765 //if(energyCLS >= 40.){
1766 //cout<<"CASE3\nPDG of the other particle VERY close: "<<last->GetPdgCode()<<" with Label: "<<last->Label();
1767 //cout<<" Energy of the other particle VERY close: "<<last->E()*TMath::Sin(last->Theta())<<endl;
1768 //}
1769 enTrue += last->E()*TMath::Sin(last->Theta());
1770 }
1771 }
1772 }
1773 else{
1774 AliAODMCParticle *first=static_cast<AliAODMCParticle*>(fAODMCParticles->At(firstidx));
1775 if((first->Phi()-phiTrue)<0.03 && (first->Eta()-etaTrue)<0.02){
1776 //if(energyCLS >= 40.){
1777 //cout<<"CASE4\nPDG of the other particle VERY close: "<<first->GetPdgCode()<<" with Label: "<<first->Label();
1778 //cout<<" Energy of the other particle VERY close: "<<first->E()*TMath::Sin(first->Theta())<<endl;
1779 //}
1780 enTrue += first->E()*TMath::Sin(first->Theta());
1781 }
1782 }
1783 }
1784 Int_t idxgrandma = MomP2Check->GetMother();
1785 AliAODMCParticle *grandma=static_cast<AliAODMCParticle*>(fAODMCParticles->At(idxgrandma));
1786 if(grandma->GetPdgCode()==111){
1787 //if(energyCLS >= 40.){
1788 //cout<<"Energy of the pi0 grandmother: "<<grandma->E()*TMath::Sin(grandma->Theta())<<endl;
1789 //}
1790 Int_t idxaunt = grandma->GetFirstDaughter();
1791 if(idxaunt<npart){
1792 if(idxaunt == mom2checkidx){
1793 Int_t auntid = grandma->GetLastDaughter();
1794 if(auntid<npart){
1795 AliAODMCParticle *lastaunt=static_cast<AliAODMCParticle*>(fAODMCParticles->At(auntid));
1796 if((lastaunt->Phi()-phiTrue)<0.03 && (lastaunt->Eta()-etaTrue)<0.02){
1797 //if(energyCLS >= 40.){
1798 //cout<<"CASE5\nPDG of the other particle VERY close: "<<lastaunt->GetPdgCode()<<" with Label: "<<lastaunt->Label();
1799 //cout<<" Energy of the other particle VERY close: "<<lastaunt->E()*TMath::Sin(lastaunt->Theta())<<endl;
1800 //}
1801 enTrue += lastaunt->E()*TMath::Sin(lastaunt->Theta());
1802 }
1803 }
1804 }
1805 else{
1806 AliAODMCParticle *aunt =static_cast<AliAODMCParticle*>(fAODMCParticles->At(idxaunt));
1807 if(aunt->GetPdgCode()==22 && (aunt->Phi()-phiTrue)<0.03 && (aunt->Eta()-etaTrue)<0.02){
1808 //if(energyCLS >= 40.){
1809 //cout<<"CASE6\nPDG of the other particle VERY close: "<<aunt->GetPdgCode()<<" with Label: "<<aunt->Label();
1810 //cout<<" Energy of the other particle VERY close: "<<aunt->E()*TMath::Sin(aunt->Theta())<<endl;
1811 //}
1812 enTrue += aunt->E()*TMath::Sin(aunt->Theta());
1813 }
1814 }
1815 }
1816 }
1817 }
1818
1819 dPhi = phiCLS-phiTrue;
1820 dEta = etaCLS-etaTrue;
1821
1822 // if(fcount==388)
1823 // AliInfo(Form("Found Particle with same label as cluster !!!! at position %d",b));
1824 // if(fcount==388){
1825 // AliInfo(Form(""));
1826 // particle2Check->Print();
1827 // cout<<"Energy of the Particle: "<<enTrue<<" Mother PDG: "<<MomP2Check->GetPdgCode()<<" Eta: "<<etaTrue<<" Phi: "<<phiTrue<<endl;
1828 //if(energyCLS >= 40.){
1829 //cout<<"Transverse Energy of all the Particle VERY CLOSE TO THe ClusterLabel Particle: "<<enTrue<<endl;
1830 //cout<<endl;
1831 //}
1832 outputvalueMCmix[0] = energyCLS;
1833 outputvalueMCmix[1] = ss;
1834 outputvalueMCmix[2] = clustPDG;
1835 outputvalueMCmix[3] = MomP2Check->GetPdgCode();
1836 outputvalueMCmix[4] = enTrue;
1837 //outputvalueMCmix[5] = dPhi;
1838 //outputvalueMCmix[6] = dEta;
1839 //outputvalueMCmix[7] = ncells;
1840 fOutClustMC->Fill(outputvalueMCmix);
1841 }
1842 // }
1843 //fPDGM02->Fill(clustPDG);
1844 //fEtrueEclustM02->Fill(energyCLS,enTrue);
1845 //fDphiDetaM02->Fill(dEta,dPhi);
1846 //fMomPDGM02->Fill(MomP2Check->GetPdgCode());
1847
1848 //if(TMath::Abs(enTrue-energyCLS)>0.2){
1849 //cout<<"Time of the cluster with energy mismatch: "<<time<<" energy of the cluster: "<<energyCLS<<" true energy: "<<enTrue<<" PDG "<<clustPDG<<" mother of the particle: "<<MomP2Check->GetPdgCode()<<endl;
1850 //fTvsE_MismatchEM02->Fill(enTrue,time);
1851 //break;
1852 //}
1853 }
1854 }
1855 if(found==kFALSE)
1856 printf("not a particle!!! Look at the STACK DOWN HERE!!!!\n\n");
1857 /*clustPDG=0;
1858 dPhi=5.;
1859 dEta=5.;
1860 enTrue = -1.;
1861 for(int b=0; b<npart; b++){
1862
1863 partMC=static_cast<AliAODMCParticle*>(fAODMCParticles->At(b));
1864 cout<<"particle "<<b<<endl;
1865 partMC->Print();
1866
1867 }
1868 fPDGM02->Fill(clustPDG);
1869 fEtrueEclustM02->Fill(energyCLS,enTrue);
1870 fDphiDetaM02->Fill(dEta,dPhi);
1871 return;
1872 }*/
1873
1874
1875
1876 //cout<<"EnergyT: "<<particle2Check->E()*TMath::Sin(particle2Check->Theta())<<"\tPDGCode: "<<particle2Check->GetPdgCode()<<"\tMotherPDG: "<<MomP2Check->GetPdgCode()<<"\tEta: "<<particle2Check->Eta()<<"\tPhi: "<<particle2Check->Phi()<<endl;
1877 //cout<<"\n\n";
1878
1879
1880 return;
1881}
1882
1883 //__________________________________________________________________________
bab35745 1884Bool_t AliAnalysisTaskEMCALPhotonIsolation::FillGeneralHistograms(AliVCluster *coi, TLorentzVector vecCOI, Int_t index){
7c97fe2d 1885 //AliInfo("Inside FillGeneralHistograms\n");
1886
1887 // Fill the histograms for underlying events and isolation studies
1888
1889// I would like to remove this part and fill the tracks multiplicity histogram in FillQAHistograms, is that ok for thnSparses? (especially cause here the histogram is filled several times per event)
bab35745 1890 AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
7c97fe2d 1891 AliEmcalParticle *emcTrack = static_cast<AliEmcalParticle*>(tracks->GetNextAcceptParticle(0));
1892
bab35745 1893 int nTracks=0;
1894 tracks->ResetCurrentID();
7c97fe2d 1895 while (emcTrack) {
bab35745 1896 AliVTrack *track = emcTrack->GetTrack();
1c662fe8 1897 if(!track) continue;
7c97fe2d 1898 // if(!(track->TestFilterBit("kHybrid"))) continue;
bab35745 1899 nTracks++;
7c97fe2d 1900 emcTrack = static_cast<AliEmcalParticle*>(tracks->GetNextAcceptParticle());
1c662fe8 1901 }
7c97fe2d 1902 fTrackMult->Fill(nTracks);
1903
1904
1905 Double_t eTCOI = 0., m02COI = 0.;
1906 //Int_t Ntracks;
1907 //Definition of the Array for Davide's Output
1c662fe8 1908 const Int_t ndims = fNDimensions;
1909 Double_t outputValues[ndims];
7c97fe2d 1910
bab35745 1911 eTCOI = vecCOI.Et();
1912 m02COI = coi->GetM02();
494479c2 1913
7c97fe2d 1914 //AliInfo(Form("M02 value: %lf\n",m02COI));
1915
1916 // ******** Isolation and UE calculation with different methods *********
1917
bab35745 1918 Double_t eTThreshold = 5;
7c97fe2d 1919
1c662fe8 1920 switch(fEtIsoMethod)
1921 {
1922 case 0: // SumEt<EtThr
7c97fe2d 1923 eTThreshold = fEtIsoThreshold;
1c662fe8 1924 break;
7c97fe2d 1925
1c662fe8 1926 case 1: // SumEt<%Ephoton
7c97fe2d 1927 eTThreshold = fEtIsoThreshold * eTCOI;
1c662fe8 1928 break;
7c97fe2d 1929
bab35745 1930 case 2: // Etmax<eTThreshold
1931 eTThreshold = fEtIsoThreshold;
7c97fe2d 1932 if( eTCOI<eTThreshold) // photon candidate, cuts have to be decided after studies
1c662fe8 1933 {
bab35745 1934 fEtIsolatedClust->Fill(eTCOI);
1c662fe8 1935 }
1936 break;
1937 }
7c97fe2d 1938
1939 //DO NOT CHANGE EVER AGAIN THE FOLLOWING DEFINITIONS
bab35745 1940 Float_t isoConeArea = TMath::Pi()*fIsoConeRadius*fIsoConeRadius;
1941 Float_t etaBandArea = (1.4 - 2.*fIsoConeRadius)*2.*fIsoConeRadius;
1942 Float_t phiBandArea = ((5./9.)*TMath::Pi()- 2.*fIsoConeRadius)*2.*fIsoConeRadius;
1943 Float_t etaBandAreaTr = (1.8 - 2.*fIsoConeRadius)*2.*fIsoConeRadius;
1944 Float_t phiBandAreaTr = (2.*TMath::Pi() - 4.*fIsoConeRadius)*2.*fIsoConeRadius;//there is a 2 more because of the JET CONE B2B
1945 Float_t perpConesArea = 2.*isoConeArea;
1946 Float_t fullTPCArea = 1.8*2.*TMath::Pi()-fIsoConeRadius*(TMath::Pi()*fIsoConeRadius + 3.6);
7c97fe2d 1947
bab35745 1948 Float_t isolation, ue;
7c97fe2d 1949
1c662fe8 1950 if(!fTPC4Iso){ //EMCAL Only for Acceptance of Clusters
1951 switch(fIsoMethod)
1952 {
1953 case 0: //EMCAL CELLS
7c97fe2d 1954
1c662fe8 1955 switch (fUEMethod)
1956 {
1957 case 0: //phi band
bab35745 1958 EtIsoCellPhiBand(vecCOI, isolation, ue);
7c97fe2d 1959 //Normalisation ue wrt Area - TO DO-
bab35745 1960 ue = ue * (isoConeArea / phiBandArea);
1961 fPhiBandUECells->Fill(vecCOI.Pt() , ue);
1962 fEtIsoCells->Fill(isolation);
1963 if(isolation<eTThreshold)
1c662fe8 1964 {
bab35745 1965 fEtIsolatedCells->Fill(eTCOI);
7c97fe2d 1966 fEtisolatedT=eTCOI;
1967 fPtisolatedT=vecCOI.Pt();
1c662fe8 1968 }
1969 break;
1970 case 1: //eta band
bab35745 1971 EtIsoCellEtaBand(vecCOI, isolation, ue);
7c97fe2d 1972 //Normalisation ue wrt Area - TO DO-
bab35745 1973 ue = ue * (isoConeArea / etaBandArea);
1974 fEtaBandUECells->Fill(vecCOI.Pt() , ue);
1975 fEtIsoCells->Fill(isolation);
1976 if(isolation<eTThreshold)
1c662fe8 1977 {
bab35745 1978 fEtIsolatedCells->Fill(eTCOI);
7c97fe2d 1979 fEtisolatedT=eTCOI;
1980 fPtisolatedT=vecCOI.Pt();
1c662fe8 1981 }
1982 break;
1983 }
1984 break;
7c97fe2d 1985
1c662fe8 1986 case 1: //EMCAL CLUSTERS
7c97fe2d 1987
1c662fe8 1988 switch (fUEMethod)
1989 {
1990 case 0: //phi band
bab35745 1991 EtIsoClusPhiBand(vecCOI, isolation, ue,index);
7c97fe2d 1992 //Normalisation ue wrt Area - TO DO-
bab35745 1993 ue = ue * (isoConeArea / phiBandArea);
1994 fPhiBandUEClust->Fill(vecCOI.Pt() , ue);
1c662fe8 1995 break;
1996 case 1: //eta band
bab35745 1997 EtIsoClusEtaBand(vecCOI, isolation, ue,index);
7c97fe2d 1998 //Normalisation ue wrt Area - TO DO-
bab35745 1999 ue = ue * (isoConeArea / etaBandArea);
2000 fEtaBandUEClust->Fill(vecCOI.Pt() , ue);
7c97fe2d 2001 break;
2002
2003 fEtIsoClust->Fill(vecCOI.Pt(),isolation);
bab35745 2004 if(isolation<eTThreshold)
1c662fe8 2005 {
7c97fe2d 2006 fPtvsM02iso->Fill(vecCOI.Pt(),coi->GetM02());
2007 fPtIsolatedNClust->Fill(vecCOI.Pt());
2008 fPtisoT=vecCOI.Pt();
2009 fM02isoT=m02COI;
2010
2011 if(fM02mincut < m02COI && m02COI < fM02maxcut)
2012 {
bab35745 2013 fEtIsolatedClust->Fill(eTCOI);
7c97fe2d 2014 fEtisolatedT=eTCOI;
2015 fPtisolatedT=vecCOI.Pt();
2016 }
2017 }
2018
2019 else
2020 {
2021 fPtvsM02noiso->Fill(vecCOI.Pt(),coi->GetM02());
2022 fPtnoisoT=vecCOI.Pt();
2023 fM02noisoT=m02COI;
1c662fe8 2024 }
1c662fe8 2025 }
7c97fe2d 2026 break;
2027
bab35745 2028 case 2: //TRACKS (TBD which tracks) //EMCAL tracks
1c662fe8 2029 switch (fUEMethod)
2030 {
2031 case 0: //phi band
bab35745 2032 PtIsoTrackPhiBand(vecCOI, isolation, ue);
7c97fe2d 2033 //Normalisation ue wrt Area - TO DO-
bab35745 2034 ue = ue * (isoConeArea / phiBandAreaTr);
2035 fPhiBandUETracks->Fill(vecCOI.Pt() , ue);
1c662fe8 2036 case 1: //eta band
bab35745 2037 PtIsoTrackEtaBand(vecCOI, isolation, ue);
7c97fe2d 2038 //Normalisation ue wrt Area - TO DO-
bab35745 2039 ue = ue * (isoConeArea / etaBandAreaTr);
2040 fEtaBandUETracks->Fill(vecCOI.Pt() , ue);
7c97fe2d 2041 break;
2042 // case 2: //Cones
2043 // PtIsoTrackOrthCones(vecCOI, absId, isolation, ue);
2044 // break;
2045 // case 3: //Full TPC
2046 // PtIsoTrackFullTPC(vecCOI, absId, isolation, ue);
2047 // break;
2048 }
2049 break;
2050
2051 // Fill histograms for isolation
2052 fPtIsoTrack->Fill(vecCOI.Pt() , isolation);
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 }
1c662fe8 2073 }
2074 }
bab35745 2075 else{ //EMCAL + TPC (Only tracks for the Isolation since IsoCone Goes Out of EMCAL)
1c662fe8 2076 switch (fUEMethod)
2077 {
2078 case 0: //phi band
bab35745 2079 PtIsoTrackPhiBand(vecCOI, isolation, ue);
7c97fe2d 2080 //Normalisation ue wrt Area - TO DO-
bab35745 2081 ue = ue * (isoConeArea / phiBandAreaTr);
2082 fPhiBandUETracks->Fill(vecCOI.Pt() , ue);
1c662fe8 2083 break;
2084 case 1: //eta band
bab35745 2085 PtIsoTrackEtaBand(vecCOI, isolation, ue);
7c97fe2d 2086 //Normalisation ue wrt Area - TO DO-
bab35745 2087 ue = ue * (isoConeArea / etaBandAreaTr);
2088 fEtaBandUETracks->Fill(vecCOI.Pt() , ue);
1c662fe8 2089 break;
2090 case 2: //Cones
bab35745 2091 PtIsoTrackOrthCones(vecCOI, isolation, ue);
7c97fe2d 2092 //Normalisation ue wrt Area - TO DO-
bab35745 2093 ue = ue * (isoConeArea / perpConesArea);
2094 fPerpConesUETracks ->Fill(vecCOI.Pt() , ue);
1c662fe8 2095 break;
2096 case 3: //Full TPC
bab35745 2097 PtIsoTrackFullTPC(vecCOI, isolation, ue);
7c97fe2d 2098 //Normalisation ue wrt Area - TO DO-
bab35745 2099 ue = ue * (isoConeArea / fullTPCArea);
2100 fTPCWithoutIsoConeB2BbandUE->Fill(vecCOI.Pt() , ue);
7c97fe2d 2101 break;
2102
2103// fill histograms for isolation
2104 fPtIsoTrack->Fill(vecCOI.Pt() , isolation);
bab35745 2105 if(isolation<eTThreshold)
1c662fe8 2106 {
7c97fe2d 2107 fPtvsM02iso->Fill(vecCOI.Pt(),coi->GetM02());
2108 fPtIsolatedNTracks->Fill(vecCOI.Pt());
2109 fPtisoT=vecCOI.Pt();
2110 fM02isoT=m02COI;
2111
2112 if(fM02mincut < m02COI && m02COI < fM02maxcut)
2113 {
2114 fEtIsolatedTracks->Fill(eTCOI);
2115 fEtisolatedT=eTCOI;
2116 fPtisolatedT=vecCOI.Pt();
2117 }
2118 }
2119 else
2120 {
2121 fPtvsM02noiso->Fill(vecCOI.Pt(),coi->GetM02());
2122 fPtnoisoT=vecCOI.Pt();
2123 fM02noisoT=m02COI;
1c662fe8 2124 }
1c662fe8 2125 }
2126 }
7c97fe2d 2127
2128
2129 /* Here we should call something to know the number of tracks...
2130 Soon I'll put in this version the "old way"; please let me know if
2131 any of you could do the same with the JET framework*/
2132
1c662fe8 2133 switch(fWho) {
2134 case 0:
7c97fe2d 2135 flambda0T=m02COI; // for all neutral clusters
2136 fEtT=vecCOI.Et(); // for all neutral clusters
2137 fPtT=vecCOI.Pt(); // for all neutral clusters
2138 fetaT=vecCOI.Eta(); // for all neutral clusters
2139 fphiT=vecCOI.Phi(); //for all neutral clusters
bab35745 2140 fsumEtisoconeT=isolation;
7c97fe2d 2141 // AliError(Form("lambda 0 %f",flambda0T));
bab35745 2142 fsumEtUE=ue;
7c97fe2d 2143
1c662fe8 2144 fOutputTree->Fill();
2145 break;
7c97fe2d 2146
1c662fe8 2147 case 1:
bab35745 2148 outputValues[0] = nTracks;
2149 outputValues[1] = eTCOI;
7c97fe2d 2150 outputValues[2] = m02COI;
bab35745 2151 outputValues[3] = isolation;
2152 outputValues[4] = ue;
2153 outputValues[5] = isolation-ue;
2154 outputValues[6] = vecCOI.Eta();
2155 outputValues[7] = vecCOI.Phi();
7c97fe2d 2156 /*if (fIsMC) {
2157 outputValues[8] = ptmc;
2158 outputValues[9] = mcptsum;
2159 }*/
1c662fe8 2160 fOutputTHnS -> Fill(outputValues);
2161 break;
7c97fe2d 2162 // // fOutPTMAX -> Fill(outputValues[1],outputValues[2],);
1c662fe8 2163 }
2164 return kTRUE;
2165}
2166
2167
7c97fe2d 2168 //_________________________________________________________________________
2169
2170void AliAnalysisTaskEMCALPhotonIsolation::AnalyzeMC(){
2171
2172 if (!fIsMC)
2173 return;
2174 //AliInfo(Form("It's a MC analysis %e",fAODMCParticles));
2175 if(!fStack && !fAODMCParticles)
2176 {cout<<"no stack saved\n"; return;}
2177 //AliInfo("there's a List of particles");
2178 //DO THIS ALSO FOR ESDs
2179
2180 Double_t eT, sumEiso, sumUE,phi, eta, distance, phip, etap, mcfirstEnergy;
2181
2182 if(fAODMCParticles->GetEntries() < 1){
2183 AliError("number of tracks insufficient");
2184 return;
2185 }
2186 int nDimMC = fMCDimensions;
2187 Double_t outputValuesMC[nDimMC];
2188
2189 Int_t nTracks = fAODMCParticles->GetEntriesFast();
2190 Int_t nFSParticles = 0;
2191 AliAODMCParticle *multTracks;
2192
2193 for(int a=0; a<nTracks; a++){
2194
2195 multTracks = static_cast<AliAODMCParticle*>(fAODMCParticles->At(a));
2196
2197 if(multTracks->IsPrimary() && multTracks->IsPhysicalPrimary() && multTracks->GetStatus()<10){
2198 if(TMath::Abs(multTracks->Eta())<=0.9 && multTracks->Charge()!= 0)
2199 nFSParticles++;
2200 else
2201 continue;
2202 }//implement final state particle condition
2203 else
2204 continue;
2205 }
2206 //AliInfo(Form("number of particles in the array %d",nTracks));
2207 AliAODMCParticle *mcpart, *mom, *mcpp,*mcsearch, *mcfirst, *mcfirstmom,*matchingtrack;
2208
2209 Bool_t prompt=kFALSE;
2210 Double_t mcEnergy, maxE, energy;
2211 Int_t pdg, mompdg, photonlabel;
2212 Double_t mcFirstEta=0., mcFirstPhi=0.;
2213
2214 Float_t isoConeArea = TMath::Pi()*fIsoConeRadius*fIsoConeRadius;
2215 Float_t etaBandArea = (1.4 - 2.*fIsoConeRadius)*2.*fIsoConeRadius;
2216 Float_t phiBandArea = ((5./9.)*TMath::Pi()- 2.*fIsoConeRadius)*2.*fIsoConeRadius;
2217 Float_t etaBandAreaTr = (1.8 - 2.*fIsoConeRadius)*2.*fIsoConeRadius;
2218 Float_t phiBandAreaTr = (2.*TMath::Pi() - 4.*fIsoConeRadius)*2.*fIsoConeRadius;//there is a 2 more because of the JET CONE B2B
2219 Float_t perpConesArea = 2.*isoConeArea;
2220 Float_t fullTPCArea = 1.8*2.*TMath::Pi()-fIsoConeRadius*(TMath::Pi()*fIsoConeRadius + 3.6);
2221
2222 // AliAODMCParticle *mcfirst = static_cast<AliAODMCParticle*>(fAODMCParticles->At(0));
2223 //AliAODMCParticle *mcp, *mcpmaxE, *mcpp, *mom;
2224 if(!fisLCAnalysis){
2225 //Loop on the event
2226 for(int iTr=0;iTr<nTracks;iTr++){
2227
2228 mcEnergy=0.;energy =0;
2229 eT=0.; phi=0.; eta=0.;
2230
2231 mcpart = static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTr));
3fe4753a 2232 /*if(mcpart->GetStatus()<10 && mcpart->IsPrimary()==0){
7c97fe2d 2233 //if(mcpart->GetMCProcessCode()!=0){
2234 if(mcpart->IsPrimary() && mcpart->IsPhysicalPrimary()){
2235
2236 if(fcount==388){
2237 cout<<"Particle in Stack: "<<iTr<<"\tLabel : "<<mcpart->Label();
2238 mcpart->Print();
2239 cout<<"IsPrimary: "<<mcpart->IsPrimary()<<"\t IsSecondaryfromDecay: "<<mcpart->IsSecondaryFromWeakDecay()<<"\t Is SecondaryfromMaterial: "<<mcpart->IsSecondaryFromMaterial()<<endl;
2240 int momidx = mcpart->GetMother();
2241 mom = static_cast<AliAODMCParticle*>(fAODMCParticles->At(momidx));
2242 cout<<"NOW THE MOTHER\n"<<momidx<<"\t";
2243 mom->Print();
2244 cout<<"IsPrimary: "<<mom->IsPrimary()<<"\t IsSecondaryfromDecay: "<<mom->IsSecondaryFromWeakDecay()<<"\t Is SecondaryfromMaterial: "<<mom->IsSecondaryFromMaterial()<<endl;
2245
2246 cout<<"\n\n"<<endl;
2247 }
2248 }
2249 }*/
2250 if(mcpart->GetStatus()>10) continue;
2251 if(!mcpart->IsPrimary()) continue;
2252 if(!mcpart->IsPhysicalPrimary()) continue;
2253
2254 pdg = mcpart->GetPdgCode();
2255
2256 if(pdg!=22)
2257 continue;
2258
2259 eta = mcpart->Eta();
2260 phi = mcpart->Phi();
2261
2262 //check photons in EMCAL //to be redefined with fIsoConeR
2263 if((TMath::Abs(eta)>0.3) || (phi<1.8 || phi>(TMath::Pi()-0.4)))
2264 continue;
2265
2266 photonlabel = iTr;
2267 int momidx = mcpart->GetMother();
2268
2269 mom = static_cast<AliAODMCParticle*>(fAODMCParticles->At(momidx));
2270 mompdg= TMath::Abs(mom->GetPdgCode());
2271
2272 eT= mcpart->E()*TMath::Sin(mcpart->Theta()); //transform to transverse Energy
2273
2274 fphietaPhotons->Fill(eta,phi,eT);
2275
2276
2277 bool foundmatch=kFALSE;
2278 for(int m=0;m<nTracks && foundmatch==kFALSE;m++){
2279 if(m==iTr) continue;
2280
2281 matchingtrack = static_cast<AliAODMCParticle*>(fAODMCParticles->At(m));
2282
2283 if(! matchingtrack->IsPrimary()) continue;
2284 if(! matchingtrack->IsPhysicalPrimary()) continue;
2285 if(matchingtrack->GetStatus()> 10 ) continue;
2286
2287 Double_t etamatching = matchingtrack->Eta();
2288 Double_t phimatching = matchingtrack->Phi();
2289
2290 if(TMath::Abs(eta-etamatching)<=0.02 && TMath::Abs(phi-phimatching)<=0.03){
2291 foundmatch=kTRUE;
2292 fphietaOthers->Fill(matchingtrack->Eta(),matchingtrack->Phi(),eT);
2293 fphietaOthersBis->Fill(matchingtrack->Eta(),matchingtrack->Phi(),matchingtrack->Pt());
2294 }
2295 }
2296
2297 //int grandmaidx = mom->GetMother();
2298
2299 /*if((mcpart->IsPrimary()) || (mompdg==22 && grandmaidx== -1)){
2300 prompt = kTRUE;
2301 }
2302 else{
2303 prompt = kFALSE;
2304 }
2305 cout<<iTr<<"\t";
2306 mcpart->Print();
2307 cout<<"IsPrimary: "<<mcpart->IsPrimary()<<"\t IsSecondaryfromDecay: "<<mcpart->IsSecondaryFromWeakDecay()<<"\t Is SecondaryfromMaterial: "<<mcpart->IsSecondaryFromMaterial()<<endl;
2308 cout<<"NOW THE MOTHER\n"<<momidx<<"\t";
2309 mom->Print();
2310 cout<<"IsPrimary: "<<mom->IsPrimary()<<"\t IsSecondaryfromDecay: "<<mom->IsSecondaryFromWeakDecay()<<"\t Is SecondaryfromMaterial: "<<mom->IsSecondaryFromMaterial()<<endl;
2311 cout<<"Coordinates of the photon (eta,phi)= ("<<eta<<","<<phi<<")"<<endl;
2312 cout<<"\n\n"<<endl;
2313 */
2314
2315 distance=0.;
2316 phip=0., etap=0.;
2317 sumEiso=0,sumUE=0;
2318
2319 for(int iTrack=0;iTrack<nTracks;iTrack++){
2320
2321 if(iTrack==photonlabel)continue;
2322
2323 mcpp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
2324
2325 if(!mcpp) continue;
2326
2327 if((mcpp->GetPdgCode())==22) continue;
2328
2329 if(mcpp->GetStatus()>10) continue;
2330
2331 phip = mcpp->Phi();
2332 etap = mcpp->Eta();
2333
2334 //Depending on which Isolation method and UE method is considered.
2335
2336 distance= TMath::Sqrt((phi-phip)*(phi-phip) + (eta-etap)*(eta-etap));
2337 //cout<<iTrack<<endl;
2338 //cout<<"Coordinates of this particle (eta,phi)= ("<<etap<<","<<phip<<")"<<endl;
2339 //cout<<"distance of this particle from the photon: "<<distance<<endl;
2340
2341 if(distance <= 0.4){ //to be changed with fIsoConeR
2342 sumEiso += mcpp->E()*TMath::Sin(mcpp->Theta());
2343 //cout<<"\n\n Transverse Energy of this particle : "<<mcpp->E()*TMath::Sin(mcpp->Theta())<<endl;
2344 }
2345 else{
2346 if(!fTPC4Iso){
2347 if(TMath::Abs(etap)>=0.7 || (phip<=1.4 || phip>= TMath::Pi()))
2348 continue;
2349 else{
2350 switch(fUEMethod){
2351 case 0: //Phi band
2352 if(TMath::Abs(eta-etap)<0.4) //to be changed with fIsoConeRadius
2353 sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2354 else continue;
2355
2356 break;
2357 case 1: //Eta band
2358 if(TMath::Abs(phi-phip)<0.4)
2359 sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2360 else continue;
2361
2362 break;
2363 }
2364 }
2365 }
2366 else{
2367 if(TMath::Abs(etap)>=1.0)
2368 continue;
2369 else{
2370 switch(fUEMethod){
2371 case 0: //Phi band
2372 {if(TMath::Abs(eta-etap)<0.4) //to be changed with fIsoConeRadius
2373 sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2374 else continue;
2375 break;
2376 }
2377 case 1: //Eta band
2378 { if(TMath::Abs(phi-phip)<0.4)
2379 sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2380 else continue;
2381
2382 break;
2383 }
2384 case 2: //Orthogonal Cones
2385 { double etacone1= eta;
2386 double etacone2= eta;
2387 double phicone1= phi - TMath::PiOver2();
2388 double phicone2= phi + TMath::PiOver2();
2389
2390 if (phicone1 < 0.) phicone1 += 2*TMath::Pi();
2391
2392 if(TMath::Sqrt(TMath::Power(etap-etacone1,2)+TMath::Power(phip-phicone1,2))< 0.4 ||
2393 TMath::Sqrt(TMath::Power(etap-etacone2,2)+TMath::Power(phip-phicone2,2))< 0.4) //to be changed with fIsoConeRadius
2394 {sumUE += mcpp->Pt();}
2395 else continue;
2396
2397 break;
2398 }
2399 case 3: //Full TPC
2400 { // Double_t phiup= phi +TMath::Pi()+fIsoConeRadius;
2401 // Double_t phidown= phi +TMath::Pi()-fIsoConeRadius;
2402 //
2403 // if(phip < phidown || phip > phiup ) //TO BE CHECKED
2404 // continue;
2405 break;
2406 }
2407 }
2408 }
2409 }
2410 }
2411 }
2412 //cout<<"\n\nTotal Energy inside the Isolation Cone : "<<sumEiso<<endl;
2413 //cout<<"Total UE Energy : "<<sumUE<<endl;
2414 if(!fTPC4Iso){
2415 switch (fUEMethod){
2416 case 0:
2417 sumUE = sumUE * (isoConeArea / phiBandArea);
2418 break;
2419 case 1:
2420 sumUE = sumUE * (isoConeArea / etaBandArea);
2421 break;
2422 }
2423 }
2424 else{
2425 switch (fUEMethod){
2426 case 0:
2427 sumUE = sumUE * (isoConeArea / phiBandAreaTr);
2428 break;
2429 case 1:
2430 sumUE = sumUE * (isoConeArea / etaBandAreaTr);
2431 break;
2432 case 2:
2433 sumUE = sumUE * (isoConeArea / perpConesArea);
2434 break;
2435 case 3:
2436 sumUE = sumUE * (isoConeArea / fullTPCArea);
2437 break;
2438 }
2439 }
2440 // cout<<"Total SCALED UE Energy : "<<sumUE<<" calculated with method "<<fUEMethod<<"which brings a normalisation factor: "<<phiBandArea<<endl;
2441
2442 outputValuesMC[0] = nFSParticles;
2443 outputValuesMC[1] = eT;
2444 //outputValuesMC[2] = sumEiso;
2445 //outputValuesMC[3] = sumUE;
2446 outputValuesMC[2] = mompdg;
2447 //outputValuesMC[5] = eta;
2448 //outputValuesMC[6] = phi;
2449 outputValuesMC[3] = mcpart->GetLabel();
2450 // EtaPhiMCPhoton
2451 // EtMC
2452 // EtIsoCone
2453 // EtMother
2454 // UE Et
2455 // Mother PDG
2456 //fill some histograms or a THnSparse or a TTree.
2457 fOutMCTruth -> Fill(outputValuesMC);
2458
2459
2460 }
2461 }
2462 else{
2463 maxE=0.;
2464 int indexmaxE=0;
2465 //getting the index of the particle with the maximum energy.
2466 for(int iTr=0;iTr<nTracks;iTr++){
2467 mcsearch= static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTr));
2468
2469 if(!mcsearch) continue;
2470
2471 if(mcsearch->GetStatus()>10) continue;
2472 if(mcsearch->GetPdgCode()!=22) continue;
2473 if(TMath::Abs(mcsearch->Eta())>0.3) continue;
2474 if(mcsearch->Phi()<= 1.8 ||mcsearch->Phi()>= TMath::Pi()) continue;
2475
2476 mcfirstEnergy= mcsearch->E();
2477 if(mcfirstEnergy>maxE){
2478 maxE=mcfirstEnergy;
2479 indexmaxE=iTr;
2480 }
2481 else continue;
2482 }
2483 mcfirst= static_cast<AliAODMCParticle*>(fAODMCParticles->At(indexmaxE));
2484 mcfirstEnergy=mcfirst->E()*TMath::Sin(mcfirst->Theta());
2485
2486 int momidx= mcfirst->GetMother();
2487 mcfirstmom = static_cast<AliAODMCParticle*>(fAODMCParticles->At(momidx));
2488 mompdg= TMath::Abs(mcfirstmom->GetPdgCode());
2489 mcFirstEta = mcfirst->Eta();
2490 mcFirstPhi = mcfirst->Phi();
2491
2492 phip=0., etap=0.;
2493 sumEiso=0,sumUE=0;
2494
2495 for(Int_t iTrack=1;iTrack<nTracks ;iTrack++){
2496 if(iTrack==indexmaxE) continue;
2497
2498 mcpp= static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
2499 phip = mcpp->Phi();
2500 etap = mcpp->Eta();
2501 if(!mcpp)
2502 continue;
2503
2504 if(mcpp->GetStatus()>10) continue;
2505 if(mcpp->GetPdgCode()==22)continue;
2506 if(TMath::Abs(etap>0.7)) continue;
2507 if(phip<=1.4 || phip>= TMath::Pi()) continue;
2508 distance=0.;
2509 distance= TMath::Sqrt((mcFirstPhi- phip)*(mcFirstPhi- phip) + (mcFirstEta- etap)*(mcFirstEta- etap));
2510
2511 if(distance<=0.4){
2512 sumEiso += mcpp->E()*TMath::Sin(mcpp->Theta());
2513 }
2514 else{
2515 if(!fTPC4Iso){
2516 if(TMath::Abs(etap)>=0.7 || (phip<=1.4 || phip>= TMath::Pi()))
2517 continue;
2518 else{
2519 switch(fUEMethod){
2520 case 0: //Phi band
2521 if(TMath::Abs(mcFirstEta-etap)<0.4) //to be changed with fIsoConeRadius
2522 sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2523 else continue;
2524
2525 break;
2526 case 1: //Eta band
2527 if(TMath::Abs(mcFirstPhi-phip)<0.4)
2528 sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2529 else continue;
2530
2531 break;
2532 }
2533 }
2534 }
2535 else{
2536 if(TMath::Abs(etap)>=1.0)
2537 continue;
2538 else{
2539 switch(fUEMethod){
2540 case 0: //Phi band
2541 { if(TMath::Abs(mcFirstEta-etap)<0.4) //to be changed with fIsoConeRadius
2542 sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2543 else continue;
2544 break;
2545 }
2546 case 1: //Eta band
2547 {if(TMath::Abs(mcFirstPhi-phip)<0.4)
2548 sumUE += mcpp->E()*TMath::Sin(mcpp->Theta());
2549 else continue;
2550
2551 break;
2552 }
2553 case 2: //Orthogonal Cones
2554 { double phicone1= mcFirstPhi - TMath::PiOver2();
2555 double phicone2= mcFirstPhi + TMath::PiOver2();
2556
2557 if (phicone1 < 0.) phicone1 += 2*TMath::Pi();
2558
2559 if(TMath::Sqrt(TMath::Power(etap-mcFirstEta,2)+TMath::Power(phip-phicone1,2))< 0.4 ||
2560 TMath::Sqrt(TMath::Power(etap-mcFirstEta,2)+TMath::Power(phip-phicone2,2))< 0.4) //to be changed with fIsoConeRadius
2561 {sumUE += mcpp->Pt();}
2562 else continue;
2563 break;
2564 }
2565 case 3: //Full TPC
2566 { // Double_t phiup= phi +TMath::Pi()+fIsoConeRadius;
2567 // Double_t phidown= phi +TMath::Pi()-fIsoConeRadius;
2568 //
2569 // if(phip < phidown || phip > phiup ) //TO BE CHECKED
2570 // continue;
2571 break;
2572 }
2573 }
2574 }
2575 }
2576 }
2577 }
2578 // cout<<"\n\nTotal Energy inside the Isolation Cone : "<<sumEiso<<endl;
2579 if(!fTPC4Iso){
2580 switch (fUEMethod){
2581 case 0:
2582 sumUE = sumUE * (isoConeArea / phiBandArea);
2583 break;
2584 case 1:
2585 sumUE = sumUE * (isoConeArea / etaBandArea);
2586 break;
2587 }
2588 }
2589 else{
2590 switch (fUEMethod){
2591 case 0:
2592 sumUE = sumUE * (isoConeArea / phiBandAreaTr);
2593 break;
2594 case 1:
2595 sumUE = sumUE * (isoConeArea / etaBandAreaTr);
2596 break;
2597 case 2:
2598 sumUE = sumUE * (isoConeArea / perpConesArea);
2599 break;
2600 case 3:
2601 sumUE = sumUE * (isoConeArea / fullTPCArea);
2602 break;
2603 }
2604 }
2605 //cout<<"Total UE Energy : "<<sumUE<<" calculated with method "<<fUEMethod<<endl;
2606
2607 //Fill the Output TTree for MC Truth
2608 }
2609
2610 return;
2611}
2612
2613
2614/*
2615 else{
2616
2617 eT = mcpmaxE->E(); //transform to transverse Energy
2618 phi = mcpmaxE->Phi();
2619 eta = mcpmaxE->Eta();
2620 distance=0.;
2621 phip=0., etap=0.;
2622 for(iTrack=0;iTrack<nTracks;iTrack++){
2623
2624 if(iTrack==maxindex)
2625 continue;
2626
2627 mcpp=static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
2628 if(!mcpp)
2629 continue;
2630
2631 phip = mcpp->Phi();
2632 etap = mcpp->Eta();
2633 distance= TMath::Sqrt((phi-phip)*(phi-phip)+(eta-etap)*(eta-etap));
2634
2635 if(distance <= 0.4) //to be changed with fIsoConeR
2636 sum += mcpp->Pt();
2637
2638 else
2639 continue;
2640 }
2641 //fill some histograms (PDG, ET, Eta/phi distributions, sum in pT)
2642 }
2643 */
1c662fe8 2644