]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/CaloTrackCorrelations/AliAnaPhoton.cxx
D_s systematic uncertainties in 0-10% and 20-50% (Anastasia)
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaPhoton.cxx
CommitLineData
85c4406e 1/**************************************************************************
1c5acb87 2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes hereby granted *
cadbb0f3 9 * without fee, provided that the above copyright notice appears in all *
1c5acb87 10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
1c5acb87 15
16//_________________________________________________________________________
17//
18// Class for the photon identification.
19// Clusters from calorimeters are identified as photons
20// and kept in the AOD. Few histograms produced.
85c4406e 21// Produces input for other analysis classes like AliAnaPi0,
22// AliAnaParticleHadronCorrelation ...
1c5acb87 23//
85c4406e 24// -- Author: Gustavo Conesa (LNF-INFN)
1c5acb87 25//////////////////////////////////////////////////////////////////////////////
85c4406e 26
27
28// --- ROOT system ---
1c5acb87 29#include <TH2F.h>
477d6cee 30#include <TClonesArray.h>
0c1383b5 31#include <TObjString.h>
123fc3bd 32#include "TParticle.h"
6175da48 33#include "TDatabasePDG.h"
1c5acb87 34
85c4406e 35// --- Analysis system ---
36#include "AliAnaPhoton.h"
1c5acb87 37#include "AliCaloTrackReader.h"
123fc3bd 38#include "AliStack.h"
1c5acb87 39#include "AliCaloPID.h"
6639984f 40#include "AliMCAnalysisUtils.h"
ff45398a 41#include "AliFiducialCut.h"
0ae57829 42#include "AliVCluster.h"
591cc579 43#include "AliAODMCParticle.h"
c8fe2783 44#include "AliMixedEvent.h"
fc195fd0 45#include "AliAODEvent.h"
2ad19c3d 46#include "AliESDEvent.h"
c8fe2783 47
85c4406e 48// --- Detectors ---
c5693f62 49#include "AliPHOSGeoUtils.h"
50#include "AliEMCALGeometry.h"
1c5acb87 51
52ClassImp(AliAnaPhoton)
521636d2 53
85c4406e 54//____________________________
55AliAnaPhoton::AliAnaPhoton() :
0cea6003 56AliAnaCaloTrackCorrBaseClass(),
85c4406e 57fMinDist(0.), fMinDist2(0.), fMinDist3(0.),
58fRejectTrackMatch(0), fFillTMHisto(kFALSE),
59fTimeCutMin(-10000), fTimeCutMax(10000),
60fNCellsCut(0),
61fNLMCutMin(-1), fNLMCutMax(10),
62fFillSSHistograms(kFALSE), fFillOnlySimpleSSHisto(1),
57cfa545 63fNOriginHistograms(9), fNPrimaryHistograms(5),
1a8c88c1 64fMomentum(), fPrimaryMom(),
85c4406e 65// Histograms
80b5ae80 66
67// Control histograms
68fhNCellsE(0), fhCellsE(0),
69fhMaxCellDiffClusterE(0), fhTimePt(0), fhEtaPhi(0),
70
85c4406e 71fhEPhoton(0), fhPtPhoton(0),
80b5ae80 72fhPhiPhoton(0), fhEtaPhoton(0),
73fhEtaPhiPhoton(0), fhEtaPhi05Photon(0),
85c4406e 74fhPtCentralityPhoton(0), fhPtEventPlanePhoton(0),
75
76// Shower shape histograms
77fhNLocMax(0),
78fhDispE(0), fhLam0E(0), fhLam1E(0),
79fhDispETRD(0), fhLam0ETRD(0), fhLam1ETRD(0),
80fhDispETM(0), fhLam0ETM(0), fhLam1ETM(0),
81fhDispETMTRD(0), fhLam0ETMTRD(0), fhLam1ETMTRD(0),
82
83fhNCellsLam0LowE(0), fhNCellsLam1LowE(0), fhNCellsDispLowE(0),
84fhNCellsLam0HighE(0), fhNCellsLam1HighE(0), fhNCellsDispHighE(0),
85
86fhEtaLam0LowE(0), fhPhiLam0LowE(0),
87fhEtaLam0HighE(0), fhPhiLam0HighE(0),
88fhLam0DispLowE(0), fhLam0DispHighE(0),
89fhLam1Lam0LowE(0), fhLam1Lam0HighE(0),
90fhDispLam1LowE(0), fhDispLam1HighE(0),
91fhDispEtaE(0), fhDispPhiE(0),
92fhSumEtaE(0), fhSumPhiE(0), fhSumEtaPhiE(0),
93fhDispEtaPhiDiffE(0), fhSphericityE(0),
94fhDispSumEtaDiffE(0), fhDispSumPhiDiffE(0),
95
96// MC histograms
97fhMCPhotonELambda0NoOverlap(0), fhMCPhotonELambda0TwoOverlap(0), fhMCPhotonELambda0NOverlap(0),
98// Embedding
99fhEmbeddedSignalFractionEnergy(0),
100fhEmbedPhotonELambda0FullSignal(0), fhEmbedPhotonELambda0MostlySignal(0),
101fhEmbedPhotonELambda0MostlyBkg(0), fhEmbedPhotonELambda0FullBkg(0),
102fhEmbedPi0ELambda0FullSignal(0), fhEmbedPi0ELambda0MostlySignal(0),
103fhEmbedPi0ELambda0MostlyBkg(0), fhEmbedPi0ELambda0FullBkg(0),
bc41680b 104
b2e375c7 105fhTimePtPhotonNoCut(0), fhTimePtPhotonSPD(0),
85c4406e 106fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
85c4406e 107fhPtPhotonNPileUpSPDVtx(0), fhPtPhotonNPileUpTrkVtx(0),
108fhPtPhotonNPileUpSPDVtxTimeCut(0), fhPtPhotonNPileUpTrkVtxTimeCut(0),
6df33fcb 109fhPtPhotonNPileUpSPDVtxTimeCut2(0), fhPtPhotonNPileUpTrkVtxTimeCut2(0),
bc41680b 110
6df33fcb 111fhEClusterSM(0), fhEPhotonSM(0),
4428fe32 112fhPtClusterSM(0), fhPtPhotonSM(0)
85c4406e 113{
1c5acb87 114 //default ctor
115
57cfa545 116 for(Int_t i = 0; i < fgkNmcTypes; i++)
4bfeae64 117 {
4c8f7c2e 118 fhMCPt [i] = 0;
119 fhMCE [i] = 0;
120 fhMCPhi [i] = 0;
121 fhMCEta [i] = 0;
85c4406e 122 fhMCDeltaE [i] = 0;
4c8f7c2e 123 fhMCDeltaPt[i] = 0;
85c4406e 124 fhMC2E [i] = 0;
4c8f7c2e 125 fhMC2Pt [i] = 0;
521636d2 126 }
127
57cfa545 128 for(Int_t i = 0; i < fgkNmcPrimTypes; i++)
4bfeae64 129 {
3d5d5078 130 fhPtPrimMC [i] = 0;
131 fhEPrimMC [i] = 0;
132 fhPhiPrimMC[i] = 0;
4cf13296 133 fhEtaPrimMC[i] = 0;
3d5d5078 134 fhYPrimMC [i] = 0;
135
136 fhPtPrimMCAcc [i] = 0;
137 fhEPrimMCAcc [i] = 0;
138 fhPhiPrimMCAcc[i] = 0;
4cf13296 139 fhEtaPrimMCAcc[i] = 0;
3d5d5078 140 fhYPrimMCAcc [i] = 0;
57cfa545 141 }
142
143 for(Int_t i = 0; i < 7; i++)
144 {
d2655d46 145 fhDispEtaDispPhi[i] = 0;
146 fhLambda0DispPhi[i] = 0;
147 fhLambda0DispEta[i] = 0;
57cfa545 148
bc41680b 149 fhPtPhotonPileUp[i] = 0;
fad96885 150 fhClusterTimeDiffPhotonPileUp [i] = 0;
57cfa545 151
152 for(Int_t j = 0; j < fgkNssTypes; j++)
d2655d46 153 {
154 fhMCDispEtaDispPhi[i][j] = 0;
155 fhMCLambda0DispEta[i][j] = 0;
156 fhMCLambda0DispPhi[i][j] = 0;
157 }
85c4406e 158 }
3d5d5078 159
57cfa545 160 for(Int_t i = 0; i < fgkNssTypes; i++)
4bfeae64 161 {
f66d95af 162 fhMCELambda0 [i] = 0;
163 fhMCELambda1 [i] = 0;
164 fhMCEDispersion [i] = 0;
85c4406e 165 fhMCNCellsE [i] = 0;
166 fhMCMaxCellDiffClusterE[i] = 0;
bfdcf7fb 167 fhLambda0DispEta[i] = 0;
168 fhLambda0DispPhi[i] = 0;
85c4406e 169
f66d95af 170 fhMCLambda0vsClusterMaxCellDiffE0[i] = 0;
171 fhMCLambda0vsClusterMaxCellDiffE2[i] = 0;
172 fhMCLambda0vsClusterMaxCellDiffE6[i] = 0;
173 fhMCNCellsvsClusterMaxCellDiffE0 [i] = 0;
174 fhMCNCellsvsClusterMaxCellDiffE2 [i] = 0;
175 fhMCNCellsvsClusterMaxCellDiffE6 [i] = 0;
34c16486 176
177 fhMCEDispEta [i] = 0;
178 fhMCEDispPhi [i] = 0;
179 fhMCESumEtaPhi [i] = 0;
180 fhMCEDispEtaPhiDiff[i] = 0;
181 fhMCESphericity [i] = 0;
521636d2 182 }
183
85c4406e 184 for(Int_t i = 0; i < 5; i++)
185 {
58ea8ce5 186 fhClusterCutsE [i] = 0;
187 fhClusterCutsPt[i] = 0;
85c4406e 188 }
189
190 // Track matching residuals
191 for(Int_t i = 0; i < 2; i++)
192 {
126b8c62 193 fhTrackMatchedDEta [i] = 0; fhTrackMatchedDPhi [i] = 0; fhTrackMatchedDEtaDPhi [i] = 0;
194 fhTrackMatchedDEtaNeg[i] = 0; fhTrackMatchedDPhiNeg[i] = 0; fhTrackMatchedDEtaDPhiNeg[i] = 0;
195 fhTrackMatchedDEtaPos[i] = 0; fhTrackMatchedDPhiPos[i] = 0; fhTrackMatchedDEtaDPhiPos[i] = 0;
85c4406e 196 fhTrackMatchedDEtaTRD[i] = 0; fhTrackMatchedDPhiTRD[i] = 0;
197 fhTrackMatchedDEtaMCOverlap[i] = 0; fhTrackMatchedDPhiMCOverlap[i] = 0;
198 fhTrackMatchedDEtaMCNoOverlap[i] = 0; fhTrackMatchedDPhiMCNoOverlap[i] = 0;
199 fhTrackMatchedDEtaMCConversion[i] = 0; fhTrackMatchedDPhiMCConversion[i] = 0;
200 fhTrackMatchedMCParticle[i] = 0; fhTrackMatchedMCParticle[i] = 0;
201 fhdEdx[i] = 0; fhEOverP[i] = 0;
202 fhEOverPTRD[i] = 0;
203 }
204
1c5acb87 205 //Initialize parameters
206 InitParameters();
85c4406e 207
1c5acb87 208}
209
b2e375c7 210//_________________________________________________________________________________________
da13b2b8 211Bool_t AliAnaPhoton::ClusterSelected(AliVCluster* calo, Int_t nMaxima)
c4a7d28a 212{
213 //Select clusters if they pass different cuts
fad96885 214
da13b2b8 215 Float_t ptcluster = fMomentum.Pt();
216 Float_t ecluster = fMomentum.E();
217 Float_t etacluster = fMomentum.Eta();
218 Float_t phicluster = fMomentum.Phi();
58ea8ce5 219
6df33fcb 220 if(phicluster < 0) phicluster+=TMath::TwoPi();
c2a62a94 221
b2e375c7 222 Bool_t matched = IsTrackMatched(calo,GetReader()->GetInputEvent());
fad96885 223
2db10729 224 AliDebug(2,Form("Current Event %d; Before selection : E %2.2f, pT %2.2f, phi %2.2f, eta %2.2f",
c4a7d28a 225 GetReader()->GetEventNumber(),
2db10729 226 ecluster,ptcluster, phicluster*TMath::RadToDeg(),etacluster));
fad96885 227
58ea8ce5 228 fhClusterCutsE [1]->Fill( ecluster);
229 fhClusterCutsPt[1]->Fill(ptcluster);
c4a7d28a 230
c2a62a94 231 if(ecluster > 0.5) fhEtaPhi->Fill(etacluster, phicluster);
232
6df33fcb 233 Int_t nSM = GetModuleNumber(calo);
234 if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0)
235 {
236 fhEClusterSM ->Fill(ecluster ,nSM);
237 fhPtClusterSM->Fill(ptcluster,nSM);
238 }
239
c4a7d28a 240 //.......................................
241 //If too small or big energy, skip it
afb3af8a 242 if(ecluster < GetMinEnergy() || ecluster > GetMaxEnergy() ) return kFALSE ;
09273901 243
2db10729 244 AliDebug(2,Form("\t Cluster %d Pass E Cut",calo->GetID()));
09273901 245
58ea8ce5 246 fhClusterCutsE [2]->Fill( ecluster);
247 fhClusterCutsPt[2]->Fill(ptcluster);
85c4406e 248
c4a7d28a 249 //.......................................
250 // TOF cut, BE CAREFUL WITH THIS CUT
251 Double_t tof = calo->GetTOF()*1e9;
252 if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
09273901 253
2db10729 254 AliDebug(2,Form("\t Cluster %d Pass Time Cut",calo->GetID()));
09273901 255
58ea8ce5 256 fhClusterCutsE [3]->Fill( ecluster);
257 fhClusterCutsPt[3]->Fill(ptcluster);
85c4406e 258
c4a7d28a 259 //.......................................
260 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
09273901 261
2db10729 262 AliDebug(2,Form("\t Cluster %d Pass NCell Cut",calo->GetID()));
09273901 263
58ea8ce5 264 fhClusterCutsE [4]->Fill( ecluster);
265 fhClusterCutsPt[4]->Fill(ptcluster);
85c4406e 266
9e51e29a 267 if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
2db10729 268 AliDebug(2,Form("\t Cluster %d pass NLM %d of out of range",calo->GetID(), nMaxima));
85c4406e 269
58ea8ce5 270 fhClusterCutsE [5]->Fill( ecluster);
271 fhClusterCutsPt[5]->Fill(ptcluster);
85c4406e 272
c4a7d28a 273 //.......................................
274 //Check acceptance selection
34c16486 275 if(IsFiducialCutOn())
276 {
da13b2b8 277 Bool_t in = GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),GetCalorimeter()) ;
c4a7d28a 278 if(! in ) return kFALSE ;
279 }
09273901 280
2db10729 281 AliDebug(2,Form("\t Fiducial cut passed"));
09273901 282
58ea8ce5 283 fhClusterCutsE [6]->Fill( ecluster);
284 fhClusterCutsPt[6]->Fill(ptcluster);
85c4406e 285
c4a7d28a 286 //.......................................
287 //Skip matched clusters with tracks
09273901 288
4bfeae64 289 // Fill matching residual histograms before PID cuts
290 if(fFillTMHisto) FillTrackMatchingResidualHistograms(calo,0);
09273901 291
34c16486 292 if(fRejectTrackMatch)
293 {
b2e375c7 294 if(matched)
34c16486 295 {
2db10729 296 AliDebug(2,"\t Reject track-matched clusters");
c4a7d28a 297 return kFALSE ;
298 }
85c4406e 299 else
2db10729 300 AliDebug(2,"\t Track-matching cut passed");
c4a7d28a 301 }// reject matched clusters
09273901 302
58ea8ce5 303 fhClusterCutsE [7]->Fill( ecluster);
304 fhClusterCutsPt[7]->Fill(ptcluster);
85c4406e 305
c4a7d28a 306 //.......................................
307 //Check Distance to Bad channel, set bit.
308 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
309 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
85c4406e 310 if(distBad < fMinDist)
34c16486 311 {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
c4a7d28a 312 return kFALSE ;
313 }
2db10729 314 else AliDebug(2,Form("\t Bad channel cut passed %4.2f > %2.2f",distBad, fMinDist));
fc195fd0 315
58ea8ce5 316 fhClusterCutsE [8]->Fill( ecluster);
317 fhClusterCutsPt[8]->Fill(ptcluster);
09273901 318
2db10729 319 AliDebug(1,Form("Current Event %d; After selection : E %2.2f, pT %2.2f, phi %2.2f, eta %2.2f",
85c4406e 320 GetReader()->GetEventNumber(),
2db10729 321 ecluster, ptcluster,fMomentum.Phi()*TMath::RadToDeg(),fMomentum.Eta()));
c4a7d28a 322
323 //All checks passed, cluster selected
324 return kTRUE;
85c4406e 325
c4a7d28a 326}
327
34c16486 328//___________________________________________
329void AliAnaPhoton::FillAcceptanceHistograms()
330{
3d5d5078 331 //Fill acceptance histograms if MC data is available
332
34c16486 333 Double_t photonY = -100 ;
334 Double_t photonE = -1 ;
335 Double_t photonPt = -1 ;
336 Double_t photonPhi = 100 ;
337 Double_t photonEta = -1 ;
85c4406e 338
34c16486 339 Int_t pdg = 0 ;
340 Int_t tag = 0 ;
f1c9c78f 341 Int_t status = 0 ;
34c16486 342 Int_t mcIndex = 0 ;
f1c9c78f 343 Int_t nprim = 0 ;
344 Bool_t inacceptance = kFALSE ;
85c4406e 345
f1c9c78f 346 TParticle * primStack = 0;
347 AliAODMCParticle * primAOD = 0;
f1c9c78f 348
349 // Get the ESD MC particles container
350 AliStack * stack = 0;
351 if( GetReader()->ReadStack() )
85c4406e 352 {
f1c9c78f 353 stack = GetMCStack();
e7757dd9 354 if( !stack )
355 {
356 AliFatal("Stack not available, is the MC handler called? STOP");
357 return;
358 }
f1c9c78f 359 nprim = stack->GetNtrack();
360 }
361
362 // Get the AOD MC particles container
363 TClonesArray * mcparticles = 0;
364 if( GetReader()->ReadAODMCParticles() )
34c16486 365 {
f1c9c78f 366 mcparticles = GetReader()->GetAODMCParticles();
e7757dd9 367 if( !mcparticles )
368 {
369 AliFatal("Standard MCParticles not available!");
370 return;
371 }
f1c9c78f 372 nprim = mcparticles->GetEntriesFast();
373 }
374
375 for(Int_t i=0 ; i < nprim; i++)
376 {
377 if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
378
379 if(GetReader()->ReadStack())
34c16486 380 {
f1c9c78f 381 primStack = stack->Particle(i) ;
28a99be3 382 if(!primStack)
383 {
2db10729 384 AliWarning("ESD primaries pointer not available!!");
28a99be3 385 continue;
386 }
387
f1c9c78f 388 pdg = primStack->GetPdgCode();
389 status = primStack->GetStatusCode();
3d5d5078 390
f1c9c78f 391 if(primStack->Energy() == TMath::Abs(primStack->Pz())) continue ; //Protection against floating point exception
392
393 //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
394 // prim->GetName(), prim->GetPdgCode());
395
396 //Photon kinematics
1a8c88c1 397 primStack->Momentum(fMomentum);
f1c9c78f 398
1f0cae67 399 photonY = 0.5*TMath::Log((primStack->Energy()+primStack->Pz())/(primStack->Energy()-primStack->Pz())) ;
f1c9c78f 400 }
401 else
402 {
403 primAOD = (AliAODMCParticle *) mcparticles->At(i);
28a99be3 404 if(!primAOD)
405 {
2db10729 406 AliWarning("AOD primaries pointer not available!!");
28a99be3 407 continue;
408 }
409
f1c9c78f 410 pdg = primAOD->GetPdgCode();
411 status = primAOD->GetStatus();
412
413 if(primAOD->E() == TMath::Abs(primAOD->Pz())) continue ; //Protection against floating point exception
414
415 //Photon kinematics
1a8c88c1 416 fMomentum.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
f1c9c78f 417
1f0cae67 418 photonY = 0.5*TMath::Log((primAOD->E()+primAOD->Pz())/(primAOD->E()-primAOD->Pz())) ;
f1c9c78f 419 }
420
421 // Select only photons in the final state
422 if(pdg != 22 ) continue ;
423
424 // If too small or too large pt, skip, same cut as for data analysis
1a8c88c1 425 photonPt = fMomentum.Pt () ;
f1c9c78f 426
427 if(photonPt < GetMinPt() || photonPt > GetMaxPt() ) continue ;
428
1a8c88c1 429 photonE = fMomentum.E () ;
430 photonEta = fMomentum.Eta() ;
431 photonPhi = fMomentum.Phi() ;
f1c9c78f 432
433 if(photonPhi < 0) photonPhi+=TMath::TwoPi();
434
435 // Check if photons hit desired acceptance
3b8377c7 436 inacceptance = kTRUE;
f1c9c78f 437
438 // Check same fidutial borders as in data analysis on top of real acceptance if real was requested.
1a8c88c1 439 if( IsFiducialCutOn() && !GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),GetCalorimeter())) inacceptance = kFALSE ;
f1c9c78f 440
441 // Check if photons hit the Calorimeter acceptance
3b8377c7 442 if(IsRealCaloAcceptanceOn()) // defined on base class
f1c9c78f 443 {
444 if(GetReader()->ReadStack() &&
0cea6003 445 !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(GetCalorimeter(), primStack)) inacceptance = kFALSE ;
f1c9c78f 446 if(GetReader()->ReadAODMCParticles() &&
0cea6003 447 !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(GetCalorimeter(), primAOD )) inacceptance = kFALSE ;
f1c9c78f 448 }
449
450 // Get tag of this particle photon from fragmentation, decay, prompt ...
451 // Set the origin of the photon.
0cea6003 452 tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(),GetCalorimeter());
f1c9c78f 453
454 if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
455 {
456 // A conversion photon from a hadron, skip this kind of photon
457 // printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!\n ");
458 // GetMCAnalysisUtils()->PrintMCTag(tag);
459
460 continue;
461 }
462
463 // Consider only final state particles, but this depends on generator,
464 // status 1 is the usual one, in case of not being ok, leave the possibility
465 // to not consider this.
466 if(status > 1) continue ; // Avoid "partonic" photons
467
468 Bool_t takeIt = kFALSE ;
da13b2b8 469 if(status == 1 && GetMCAnalysisUtils()->GetMCGenerator() != AliMCAnalysisUtils::kBoxLike ) takeIt = kTRUE ;
f1c9c78f 470
57cfa545 471 if (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) continue;
f1c9c78f 472
473 //Origin of photon
760b98f5 474 if (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))
f1c9c78f 475 {
476 mcIndex = kmcPPrompt;
477 }
478 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
479 {
480 mcIndex = kmcPFragmentation ;
481 }
482 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
483 {
484 mcIndex = kmcPISR;
485 }
486 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
487 {
488 mcIndex = kmcPPi0Decay;
489 }
57cfa545 490 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
491 {
492 mcIndex = kmcPEtaDecay;
493 }
494 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
f1c9c78f 495 {
496 mcIndex = kmcPOtherDecay;
497 }
498 else
499 {
500 // Other decay but from non final state particle
501 mcIndex = kmcPOtherDecay;
502 }//Other origin
503
504 if(!takeIt && (mcIndex == kmcPPi0Decay || mcIndex == kmcPOtherDecay)) takeIt = kTRUE ;
505
506 if(!takeIt) continue ;
3d5d5078 507
3b8377c7 508 //Fill histograms for all photons
f1c9c78f 509 fhYPrimMC[kmcPPhoton]->Fill(photonPt, photonY) ;
510 if(TMath::Abs(photonY) < 1.0)
511 {
512 fhEPrimMC [kmcPPhoton]->Fill(photonE ) ;
513 fhPtPrimMC [kmcPPhoton]->Fill(photonPt) ;
514 fhPhiPrimMC[kmcPPhoton]->Fill(photonE , photonPhi) ;
515 fhEtaPrimMC[kmcPPhoton]->Fill(photonE , photonEta) ;
516 }
3b8377c7 517
f1c9c78f 518 if(inacceptance)
519 {
520 fhEPrimMCAcc [kmcPPhoton]->Fill(photonE ) ;
521 fhPtPrimMCAcc [kmcPPhoton]->Fill(photonPt) ;
522 fhPhiPrimMCAcc[kmcPPhoton]->Fill(photonE , photonPhi) ;
523 fhEtaPrimMCAcc[kmcPPhoton]->Fill(photonE , photonEta) ;
524 fhYPrimMCAcc [kmcPPhoton]->Fill(photonE , photonY) ;
525 }//Accepted
526
3b8377c7 527 //Fill histograms for photons origin
760b98f5 528 if(mcIndex < fNPrimaryHistograms)
f1c9c78f 529 {
760b98f5 530 fhYPrimMC[mcIndex]->Fill(photonPt, photonY) ;
531 if(TMath::Abs(photonY) < 1.0)
532 {
533 fhEPrimMC [mcIndex]->Fill(photonE ) ;
534 fhPtPrimMC [mcIndex]->Fill(photonPt) ;
535 fhPhiPrimMC[mcIndex]->Fill(photonE , photonPhi) ;
536 fhEtaPrimMC[mcIndex]->Fill(photonE , photonEta) ;
537 }
538
539 if(inacceptance)
540 {
541 fhEPrimMCAcc [mcIndex]->Fill(photonE ) ;
542 fhPtPrimMCAcc [mcIndex]->Fill(photonPt) ;
543 fhPhiPrimMCAcc[mcIndex]->Fill(photonE , photonPhi) ;
544 fhEtaPrimMCAcc[mcIndex]->Fill(photonE , photonEta) ;
545 fhYPrimMCAcc [mcIndex]->Fill(photonE , photonY) ;
546 }//Accepted
f1c9c78f 547 }
f1c9c78f 548 }//loop on primaries
549
3d5d5078 550}
521636d2 551
bc41680b 552//________________________________________________________________________________
26569ec8 553void AliAnaPhoton::FillPileUpHistograms(AliVCluster* cluster, AliVCaloCells *cells,
554 Int_t absIdMax)
b2e375c7 555{
bc41680b 556 // Fill some histograms to understand pile-up
b2e375c7 557
da13b2b8 558 Float_t pt = fMomentum.Pt();
bc41680b 559 Float_t time = cluster->GetTOF()*1.e9;
b2e375c7 560
bc41680b 561 AliVEvent * event = GetReader()->GetInputEvent();
b2e375c7 562
bc41680b 563 if(GetReader()->IsPileUpFromSPD()) fhPtPhotonPileUp[0]->Fill(pt);
564 if(GetReader()->IsPileUpFromEMCal()) fhPtPhotonPileUp[1]->Fill(pt);
565 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtPhotonPileUp[2]->Fill(pt);
566 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtPhotonPileUp[3]->Fill(pt);
567 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtPhotonPileUp[4]->Fill(pt);
568 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtPhotonPileUp[5]->Fill(pt);
569 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPhotonPileUp[6]->Fill(pt);
b2e375c7 570
bc41680b 571 fhTimePtPhotonNoCut->Fill(pt,time);
572 if(GetReader()->IsPileUpFromSPD()) fhTimePtPhotonSPD->Fill(pt,time);
b2e375c7 573
bc41680b 574 // cells inside the cluster
26569ec8 575
126b8c62 576 //Loop on cells inside cluster, max cell must be over 100 MeV and time in BC=0
bc41680b 577 if(cells->GetCellAmplitude(absIdMax) > 0.1 && TMath::Abs(time) < 30)
b2e375c7 578 {
bc41680b 579 for (Int_t ipos = 0; ipos < cluster->GetNCells(); ipos++)
b2e375c7 580 {
bc41680b 581 Int_t absId = cluster->GetCellsAbsId()[ipos];
126b8c62 582
583 if( absId == absIdMax ) continue ;
584
bc41680b 585 Double_t tcell = cells->GetCellTime(absId);
b2e375c7 586 Float_t amp = cells->GetCellAmplitude(absId);
587 Int_t bc = GetReader()->GetInputEvent()->GetBunchCrossNumber();
b2e375c7 588
bc41680b 589 GetCaloUtils()->GetEMCALRecoUtils()->AcceptCalibrateCell(absId,bc,amp,tcell,cells);
590 tcell*=1e9;
b2e375c7 591
bc41680b 592 Float_t diff = (time-tcell);
b2e375c7 593
36769d30 594 if( cells->GetCellAmplitude(absIdMax) < 0.1 ) continue ;
126b8c62 595
bc41680b 596 if(GetReader()->IsPileUpFromSPD()) fhClusterTimeDiffPhotonPileUp[0]->Fill(pt, diff);
597 if(GetReader()->IsPileUpFromEMCal()) fhClusterTimeDiffPhotonPileUp[1]->Fill(pt, diff);
598 if(GetReader()->IsPileUpFromSPDOrEMCal()) fhClusterTimeDiffPhotonPileUp[2]->Fill(pt, diff);
599 if(GetReader()->IsPileUpFromSPDAndEMCal()) fhClusterTimeDiffPhotonPileUp[3]->Fill(pt, diff);
600 if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhClusterTimeDiffPhotonPileUp[4]->Fill(pt, diff);
601 if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhClusterTimeDiffPhotonPileUp[5]->Fill(pt, diff);
602 if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhClusterTimeDiffPhotonPileUp[6]->Fill(pt, diff);
b2e375c7 603
bc41680b 604 }//loop
85c4406e 605 }
acd56ca4 606
2ad19c3d 607 AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
608 AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
609
610 // N pile up vertices
0f7e7205 611 Int_t nVtxSPD = -1;
612 Int_t nVtxTrk = -1;
2ad19c3d 613
614 if (esdEv)
615 {
0f7e7205 616 nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
617 nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
85c4406e 618
2ad19c3d 619 }//ESD
620 else if (aodEv)
621 {
0f7e7205 622 nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
623 nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
2ad19c3d 624 }//AOD
625
bc41680b 626 if(pt < 8)
627 {
628 fhTimeNPileUpVertSPD ->Fill(time,nVtxSPD);
629 fhTimeNPileUpVertTrack->Fill(time,nVtxTrk);
630 }
2ad19c3d 631
85c4406e 632 fhPtPhotonNPileUpSPDVtx->Fill(pt,nVtxSPD);
0f7e7205 633 fhPtPhotonNPileUpTrkVtx->Fill(pt,nVtxTrk);
634
635 if(TMath::Abs(time) < 25)
85c4406e 636 {
637 fhPtPhotonNPileUpSPDVtxTimeCut->Fill(pt,nVtxSPD);
638 fhPtPhotonNPileUpTrkVtxTimeCut->Fill(pt,nVtxTrk);
0f7e7205 639 }
640
85c4406e 641 if(time < 75 && time > -25)
642 {
643 fhPtPhotonNPileUpSPDVtxTimeCut2->Fill(pt,nVtxSPD);
644 fhPtPhotonNPileUpTrkVtxTimeCut2->Fill(pt,nVtxTrk);
645 }
646
2ad19c3d 647}
648
26569ec8 649//_______________________________________________________________________________
650void AliAnaPhoton::FillShowerShapeHistograms(AliVCluster* cluster,
b8e951b5 651 Int_t mcTag, Float_t maxCellFraction)
34c16486 652{
85c4406e 653 //Fill cluster Shower Shape histograms
521636d2 654
655 if(!fFillSSHistograms || GetMixedEvent()) return;
85c4406e 656
521636d2 657 Float_t energy = cluster->E();
658 Int_t ncells = cluster->GetNCells();
521636d2 659 Float_t lambda0 = cluster->GetM02();
660 Float_t lambda1 = cluster->GetM20();
661 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
662
da13b2b8 663 Float_t eta = fMomentum.Eta();
664 Float_t phi = fMomentum.Phi();
521636d2 665 if(phi < 0) phi+=TMath::TwoPi();
666
667 fhLam0E ->Fill(energy,lambda0);
668 fhLam1E ->Fill(energy,lambda1);
669 fhDispE ->Fill(energy,disp);
85c4406e 670
bdb39dbd 671 if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
4d1d8f00 672 GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD() )
34c16486 673 {
521636d2 674 fhLam0ETRD->Fill(energy,lambda0);
675 fhLam1ETRD->Fill(energy,lambda1);
676 fhDispETRD->Fill(energy,disp);
521636d2 677 }
678
34c16486 679 Float_t l0 = 0., l1 = 0.;
85c4406e 680 Float_t dispp= 0., dEta = 0., dPhi = 0.;
681 Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
bdb39dbd 682 if(GetCalorimeter() == kEMCAL && !fFillOnlySimpleSSHisto)
34c16486 683 {
684 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
685 l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
686 //printf("AliAnaPhoton::FillShowerShapeHistogram - l0 %2.6f, l1 %2.6f, disp %2.6f, dEta %2.6f, dPhi %2.6f, sEta %2.6f, sPhi %2.6f, sEtaPhi %2.6f \n",
687 // l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi );
688 //printf("AliAnaPhoton::FillShowerShapeHistogram - dispersion %f, dispersion eta+phi %f \n",
689 // disp, dPhi+dEta );
690 fhDispEtaE -> Fill(energy,dEta);
691 fhDispPhiE -> Fill(energy,dPhi);
692 fhSumEtaE -> Fill(energy,sEta);
693 fhSumPhiE -> Fill(energy,sPhi);
694 fhSumEtaPhiE -> Fill(energy,sEtaPhi);
695 fhDispEtaPhiDiffE -> Fill(energy,dPhi-dEta);
696 if(dEta+dPhi>0)fhSphericityE -> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
697 if(dEta+sEta>0)fhDispSumEtaDiffE -> Fill(energy,(dEta-sEta)/((dEta+sEta)/2.));
85c4406e 698 if(dPhi+sPhi>0)fhDispSumPhiDiffE -> Fill(energy,(dPhi-sPhi)/((dPhi+sPhi)/2.));
34c16486 699
bfdcf7fb 700 Int_t ebin = -1;
701 if (energy < 2 ) ebin = 0;
702 else if (energy < 4 ) ebin = 1;
703 else if (energy < 6 ) ebin = 2;
704 else if (energy < 10) ebin = 3;
85c4406e 705 else if (energy < 15) ebin = 4;
706 else if (energy < 20) ebin = 5;
707 else ebin = 6;
bfdcf7fb 708
709 fhDispEtaDispPhi[ebin]->Fill(dEta ,dPhi);
710 fhLambda0DispEta[ebin]->Fill(lambda0,dEta);
711 fhLambda0DispPhi[ebin]->Fill(lambda0,dPhi);
712
34c16486 713 }
714
85c4406e 715 // if track-matching was of, check effect of track-matching residual cut
b5dbb99b 716
717 if(!fRejectTrackMatch)
718 {
719 Float_t dZ = cluster->GetTrackDz();
720 Float_t dR = cluster->GetTrackDx();
34c16486 721 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
722 {
b5dbb99b 723 dR = 2000., dZ = 2000.;
724 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
85c4406e 725 }
b5dbb99b 726
727 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
728 {
729 fhLam0ETM ->Fill(energy,lambda0);
730 fhLam1ETM ->Fill(energy,lambda1);
731 fhDispETM ->Fill(energy,disp);
732
bdb39dbd 733 if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
4d1d8f00 734 GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD() )
34c16486 735 {
b5dbb99b 736 fhLam0ETMTRD->Fill(energy,lambda0);
737 fhLam1ETMTRD->Fill(energy,lambda1);
738 fhDispETMTRD->Fill(energy,disp);
739 }
740 }
85c4406e 741 }// if track-matching was of, check effect of matching residual cut
b5dbb99b 742
764ab1f4 743
4301207e 744 if(!fFillOnlySimpleSSHisto)
745 {
764ab1f4 746 if(energy < 2)
747 {
748 fhNCellsLam0LowE ->Fill(ncells,lambda0);
749 fhNCellsLam1LowE ->Fill(ncells,lambda1);
750 fhNCellsDispLowE ->Fill(ncells,disp);
751
752 fhLam1Lam0LowE ->Fill(lambda1,lambda0);
753 fhLam0DispLowE ->Fill(lambda0,disp);
754 fhDispLam1LowE ->Fill(disp,lambda1);
755 fhEtaLam0LowE ->Fill(eta,lambda0);
85c4406e 756 fhPhiLam0LowE ->Fill(phi,lambda0);
764ab1f4 757 }
85c4406e 758 else
764ab1f4 759 {
760 fhNCellsLam0HighE ->Fill(ncells,lambda0);
761 fhNCellsLam1HighE ->Fill(ncells,lambda1);
762 fhNCellsDispHighE ->Fill(ncells,disp);
763
764 fhLam1Lam0HighE ->Fill(lambda1,lambda0);
765 fhLam0DispHighE ->Fill(lambda0,disp);
766 fhDispLam1HighE ->Fill(disp,lambda1);
767 fhEtaLam0HighE ->Fill(eta, lambda0);
768 fhPhiLam0HighE ->Fill(phi, lambda0);
769 }
521636d2 770 }
4301207e 771
34c16486 772 if(IsDataMC())
773 {
f66d95af 774 AliVCaloCells* cells = 0;
bdb39dbd 775 if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
f66d95af 776 else cells = GetPHOSCells();
3d5d5078 777
778 //Fill histograms to check shape of embedded clusters
779 Float_t fraction = 0;
85c4406e 780 // printf("check embedding %i\n",GetReader()->IsEmbeddedClusterSelectionOn());
781
34c16486 782 if(GetReader()->IsEmbeddedClusterSelectionOn())
783 {//Only working for EMCAL
85c4406e 784 // printf("embedded\n");
3d5d5078 785 Float_t clusterE = 0; // recalculate in case corrections applied.
786 Float_t cellE = 0;
34c16486 787 for(Int_t icell = 0; icell < cluster->GetNCells(); icell++)
788 {
3d5d5078 789 cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
85c4406e 790 clusterE+=cellE;
3d5d5078 791 fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
792 }
793
794 //Fraction of total energy due to the embedded signal
795 fraction/=clusterE;
796
2db10729 797 AliDebug(1,Form("Energy fraction of embedded signal %2.3f, Energy %2.3f",fraction, clusterE));
3d5d5078 798
799 fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
800
85c4406e 801 } // embedded fraction
26569ec8 802
f66d95af 803 // Check the origin and fill histograms
34c16486 804
805 Int_t mcIndex = -1;
806
85c4406e 807 if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
3d5d5078 808 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
34c16486 809 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
810 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
811 {
812 mcIndex = kmcssPhoton ;
85c4406e 813
34c16486 814 if(!GetReader()->IsEmbeddedClusterSelectionOn())
815 {
3d5d5078 816 //Check particle overlaps in cluster
817
85c4406e 818 // Compare the primary depositing more energy with the rest,
8d6b7f60 819 // if no photon/electron as comon ancestor (conversions), count as other particle
4914e781 820 const UInt_t nlabels = cluster->GetNLabels();
821 Int_t overpdg[nlabels];
822 Int_t noverlaps = GetMCAnalysisUtils()->GetNOverlaps(cluster->GetLabels(), nlabels,mcTag,-1,GetReader(),overpdg);
823
8d6b7f60 824 //printf("N overlaps %d \n",noverlaps);
3d5d5078 825
f27fe026 826 if(noverlaps == 0)
34c16486 827 {
3d5d5078 828 fhMCPhotonELambda0NoOverlap ->Fill(energy, lambda0);
3d5d5078 829 }
f27fe026 830 else if(noverlaps == 1)
85c4406e 831 {
3d5d5078 832 fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0);
3d5d5078 833 }
f27fe026 834 else if(noverlaps > 1)
85c4406e 835 {
3d5d5078 836 fhMCPhotonELambda0NOverlap ->Fill(energy, lambda0);
3d5d5078 837 }
85c4406e 838 else
34c16486 839 {
2db10729 840 AliWarning(Form("n overlaps = %d!!", noverlaps));
3d5d5078 841 }
842 }//No embedding
843
844 //Fill histograms to check shape of embedded clusters
34c16486 845 if(GetReader()->IsEmbeddedClusterSelectionOn())
846 {
85c4406e 847 if (fraction > 0.9)
3d5d5078 848 {
849 fhEmbedPhotonELambda0FullSignal ->Fill(energy, lambda0);
3d5d5078 850 }
851 else if(fraction > 0.5)
852 {
853 fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0);
3d5d5078 854 }
855 else if(fraction > 0.1)
85c4406e 856 {
3d5d5078 857 fhEmbedPhotonELambda0MostlyBkg ->Fill(energy, lambda0);
3d5d5078 858 }
859 else
860 {
861 fhEmbedPhotonELambda0FullBkg ->Fill(energy, lambda0);
3d5d5078 862 }
863 } // embedded
864
521636d2 865 }//photon no conversion
4431f13a 866 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
867 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
868 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
869 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
870 {
871 mcIndex = kmcssConversion ;
872 }//conversion photon
873
34c16486 874 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))
875 {
876 mcIndex = kmcssElectron ;
521636d2 877 }//electron
34c16486 878 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) )
879 {
880 mcIndex = kmcssPi0 ;
3d5d5078 881
882 //Fill histograms to check shape of embedded clusters
34c16486 883 if(GetReader()->IsEmbeddedClusterSelectionOn())
884 {
85c4406e 885 if (fraction > 0.9)
3d5d5078 886 {
887 fhEmbedPi0ELambda0FullSignal ->Fill(energy, lambda0);
3d5d5078 888 }
889 else if(fraction > 0.5)
890 {
891 fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0);
3d5d5078 892 }
893 else if(fraction > 0.1)
85c4406e 894 {
3d5d5078 895 fhEmbedPi0ELambda0MostlyBkg ->Fill(energy, lambda0);
3d5d5078 896 }
897 else
898 {
899 fhEmbedPi0ELambda0FullBkg ->Fill(energy, lambda0);
3d5d5078 900 }
85c4406e 901 } // embedded
3d5d5078 902
521636d2 903 }//pi0
34c16486 904 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) )
905 {
906 mcIndex = kmcssEta ;
85c4406e 907 }//eta
908 else
34c16486 909 {
85c4406e 910 mcIndex = kmcssOther ;
911 }//other particles
521636d2 912
34c16486 913 fhMCELambda0 [mcIndex]->Fill(energy, lambda0);
914 fhMCELambda1 [mcIndex]->Fill(energy, lambda1);
915 fhMCEDispersion [mcIndex]->Fill(energy, disp);
916 fhMCNCellsE [mcIndex]->Fill(energy, ncells);
917 fhMCMaxCellDiffClusterE[mcIndex]->Fill(energy, maxCellFraction);
918
764ab1f4 919 if(!fFillOnlySimpleSSHisto)
34c16486 920 {
764ab1f4 921 if (energy < 2.)
922 {
923 fhMCLambda0vsClusterMaxCellDiffE0[mcIndex]->Fill(lambda0, maxCellFraction);
924 fhMCNCellsvsClusterMaxCellDiffE0 [mcIndex]->Fill(ncells, maxCellFraction);
925 }
926 else if(energy < 6.)
927 {
928 fhMCLambda0vsClusterMaxCellDiffE2[mcIndex]->Fill(lambda0, maxCellFraction);
929 fhMCNCellsvsClusterMaxCellDiffE2 [mcIndex]->Fill(ncells, maxCellFraction);
930 }
931 else
932 {
933 fhMCLambda0vsClusterMaxCellDiffE6[mcIndex]->Fill(lambda0, maxCellFraction);
934 fhMCNCellsvsClusterMaxCellDiffE6 [mcIndex]->Fill(ncells, maxCellFraction);
935 }
936
bdb39dbd 937 if(GetCalorimeter() == kEMCAL)
764ab1f4 938 {
939 fhMCEDispEta [mcIndex]-> Fill(energy,dEta);
940 fhMCEDispPhi [mcIndex]-> Fill(energy,dPhi);
941 fhMCESumEtaPhi [mcIndex]-> Fill(energy,sEtaPhi);
942 fhMCEDispEtaPhiDiff [mcIndex]-> Fill(energy,dPhi-dEta);
85c4406e 943 if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
764ab1f4 944
945 Int_t ebin = -1;
946 if (energy < 2 ) ebin = 0;
947 else if (energy < 4 ) ebin = 1;
948 else if (energy < 6 ) ebin = 2;
949 else if (energy < 10) ebin = 3;
85c4406e 950 else if (energy < 15) ebin = 4;
951 else if (energy < 20) ebin = 5;
952 else ebin = 6;
764ab1f4 953
954 fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta ,dPhi);
955 fhMCLambda0DispEta[ebin][mcIndex]->Fill(lambda0,dEta);
85c4406e 956 fhMCLambda0DispPhi[ebin][mcIndex]->Fill(lambda0,dPhi);
764ab1f4 957 }
34c16486 958 }
521636d2 959 }//MC data
960
961}
962
4bfeae64 963//__________________________________________________________________________
85c4406e 964void AliAnaPhoton::FillTrackMatchingResidualHistograms(AliVCluster* cluster,
22ad7981 965 Int_t cut)
4bfeae64 966{
967 // If selected, fill histograms with residuals of matched clusters, help to define track matching cut
968 // Residual filled for different cuts 0 (No cut), after 1 PID cut
85c4406e 969
4bfeae64 970 Float_t dZ = cluster->GetTrackDz();
971 Float_t dR = cluster->GetTrackDx();
972
973 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
974 {
975 dR = 2000., dZ = 2000.;
976 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
85c4406e 977 }
978
b2e375c7 979 AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
980
981 Bool_t positive = kFALSE;
982 if(track) positive = (track->Charge()>0);
983
4bfeae64 984 if(fhTrackMatchedDEta[cut] && TMath::Abs(dR) < 999)
985 {
986 fhTrackMatchedDEta[cut]->Fill(cluster->E(),dZ);
987 fhTrackMatchedDPhi[cut]->Fill(cluster->E(),dR);
4bfeae64 988 if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhi[cut]->Fill(dZ,dR);
b2e375c7 989
990 if(track)
991 {
992 if(positive)
993 {
994 fhTrackMatchedDEtaPos[cut]->Fill(cluster->E(),dZ);
995 fhTrackMatchedDPhiPos[cut]->Fill(cluster->E(),dR);
996 if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiPos[cut]->Fill(dZ,dR);
997 }
998 else
999 {
1000 fhTrackMatchedDEtaNeg[cut]->Fill(cluster->E(),dZ);
1001 fhTrackMatchedDPhiNeg[cut]->Fill(cluster->E(),dR);
1002 if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiNeg[cut]->Fill(dZ,dR);
1003 }
1004 }
4bfeae64 1005
1006 Int_t nSMod = GetModuleNumber(cluster);
1007
bdb39dbd 1008 if(GetCalorimeter()==kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
4d1d8f00 1009 nSMod >= GetFirstSMCoveredByTRD() )
4bfeae64 1010 {
1011 fhTrackMatchedDEtaTRD[cut]->Fill(cluster->E(),dZ);
1012 fhTrackMatchedDPhiTRD[cut]->Fill(cluster->E(),dR);
1013 }
1014
1015 // Check dEdx and E/p of matched clusters
b2e375c7 1016
4bfeae64 1017 if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
85c4406e 1018 {
85c4406e 1019 if(track)
4bfeae64 1020 {
4bfeae64 1021 Float_t dEdx = track->GetTPCsignal();
1022 Float_t eOverp = cluster->E()/track->P();
1023
1024 fhdEdx[cut] ->Fill(cluster->E(), dEdx);
1025 fhEOverP[cut]->Fill(cluster->E(), eOverp);
1026
bdb39dbd 1027 if(GetCalorimeter()==kEMCAL && GetFirstSMCoveredByTRD() >= 0 &&
4d1d8f00 1028 nSMod >= GetFirstSMCoveredByTRD() )
4bfeae64 1029 fhEOverPTRD[cut]->Fill(cluster->E(), eOverp);
1030
1031
1032 }
1033 else
2db10729 1034 AliWarning(Form("Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT?", dR,dZ));
4bfeae64 1035
1036 if(IsDataMC())
1037 {
1038
0cea6003 1039 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(),GetReader(),GetCalorimeter());
4bfeae64 1040
1041 if ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
1042 {
1043 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
85c4406e 1044 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 2.5 );
4bfeae64 1045 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 0.5 );
1046 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 1.5 );
1047 else fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 3.5 );
1048
1049 // Check if several particles contributed to cluster and discard overlapped mesons
85c4406e 1050 if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
34c16486 1051 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
1052 {
4bfeae64 1053 if(cluster->GetNLabels()==1)
1054 {
1055 fhTrackMatchedDEtaMCNoOverlap[cut]->Fill(cluster->E(),dZ);
1056 fhTrackMatchedDPhiMCNoOverlap[cut]->Fill(cluster->E(),dR);
1057 }
85c4406e 1058 else
4bfeae64 1059 {
1060 fhTrackMatchedDEtaMCOverlap[cut]->Fill(cluster->E(),dZ);
1061 fhTrackMatchedDPhiMCOverlap[cut]->Fill(cluster->E(),dR);
1062 }
1063
1064 }// Check overlaps
1065
1066 }
1067 else
1068 {
1069 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
85c4406e 1070 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 6.5 );
4bfeae64 1071 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 4.5 );
1072 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) ) fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 5.5 );
1073 else fhTrackMatchedMCParticle[cut]->Fill(cluster->E(), 7.5 );
1074
1075 // Check if several particles contributed to cluster
85c4406e 1076 if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ||
34c16486 1077 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
1078 {
4bfeae64 1079 fhTrackMatchedDEtaMCConversion[cut]->Fill(cluster->E(),dZ);
1080 fhTrackMatchedDPhiMCConversion[cut]->Fill(cluster->E(),dR);
1081
85c4406e 1082 }// Check overlaps
4bfeae64 1083
1084 }
1085
85c4406e 1086 } // MC
4bfeae64 1087
1088 } // residuals window
1089
1090 } // Small residual
1091
1092}
1093
1094//___________________________________________
0c1383b5 1095TObjString * AliAnaPhoton::GetAnalysisCuts()
85c4406e 1096{
0c1383b5 1097 //Save parameters used for analysis
1098 TString parList ; //this will be list of parameters used for this analysis.
5ae09196 1099 const Int_t buffersize = 255;
1100 char onePar[buffersize] ;
0c1383b5 1101
2db10729 1102 snprintf(onePar,buffersize,"--- AliAnaPhoton ---:") ;
85c4406e 1103 parList+=onePar ;
2db10729 1104 snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
0c1383b5 1105 parList+=onePar ;
2db10729 1106 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster);",fMinDist) ;
0c1383b5 1107 parList+=onePar ;
2db10729 1108 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation);",fMinDist2) ;
0c1383b5 1109 parList+=onePar ;
2db10729 1110 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study);",fMinDist3) ;
0c1383b5 1111 parList+=onePar ;
2db10729 1112 snprintf(onePar,buffersize,"fRejectTrackMatch: %d",fRejectTrackMatch) ;
85c4406e 1113 parList+=onePar ;
0c1383b5 1114
1115 //Get parameters set in base class.
1116 parList += GetBaseParametersList() ;
1117
1118 //Get parameters set in PID class.
1119 parList += GetCaloPID()->GetPIDParametersList() ;
1120
1121 //Get parameters set in FiducialCut class (not available yet)
85c4406e 1122 //parlist += GetFidCut()->GetFidCutParametersList()
0c1383b5 1123
1124 return new TObjString(parList) ;
1125}
1126
1c5acb87 1127//________________________________________________________________________
1128TList * AliAnaPhoton::GetCreateOutputObjects()
c2a62a94 1129{
85c4406e 1130 // Create histograms to be saved in output file and
477d6cee 1131 // store them in outputContainer
85c4406e 1132 TList * outputContainer = new TList() ;
1133 outputContainer->SetName("PhotonHistos") ;
4a745797 1134
85c4406e 1135 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
1136 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Float_t phimax = GetHistogramRanges()->GetHistoPhiMax(); Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
1137 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins(); Float_t etamax = GetHistogramRanges()->GetHistoEtaMax(); Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
745913ae 1138 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
85c4406e 1139 Int_t nbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t nmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t nmin = GetHistogramRanges()->GetHistoNClusterCellMin();
1140 Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins(); Float_t timemax = GetHistogramRanges()->GetHistoTimeMax(); Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
1141
1142 Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
1143 Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
09273901 1144 Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
85c4406e 1145 Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
1146 Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
09273901 1147 Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
1148
85c4406e 1149 Int_t ndedxbins = GetHistogramRanges()->GetHistodEdxBins();
1150 Float_t dedxmax = GetHistogramRanges()->GetHistodEdxMax();
31ae6d59 1151 Float_t dedxmin = GetHistogramRanges()->GetHistodEdxMin();
85c4406e 1152 Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
1153 Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
31ae6d59 1154 Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
09273901 1155
d2655d46 1156 Int_t bin[] = {0,2,4,6,10,15,20,100}; // energy bins for SS studies
1157
9e51e29a 1158 TString cut[] = {"Open","Reader","E","Time","NCells","NLM","Fidutial","Matching","Bad","PID"};
85c4406e 1159 for (Int_t i = 0; i < 10 ; i++)
fc195fd0 1160 {
58ea8ce5 1161 fhClusterCutsE[i] = new TH1F(Form("hE_Cut_%d_%s", i, cut[i].Data()),
fc195fd0 1162 Form("Number of clusters that pass cuts <= %d, %s", i, cut[i].Data()),
85c4406e 1163 nptbins,ptmin,ptmax);
58ea8ce5 1164 fhClusterCutsE[i]->SetYTitle("d#it{N}/d#it{E} ");
1165 fhClusterCutsE[i]->SetXTitle("#it{E} (GeV)");
1166 outputContainer->Add(fhClusterCutsE[i]) ;
1167
1168 fhClusterCutsPt[i] = new TH1F(Form("hPt_Cut_%d_%s", i, cut[i].Data()),
1169 Form("Number of clusters that pass cuts <= %d, %s", i, cut[i].Data()),
1170 nptbins,ptmin,ptmax);
1171 fhClusterCutsPt[i]->SetYTitle("d#it{N}/d#it{E} ");
1172 fhClusterCutsPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1173 outputContainer->Add(fhClusterCutsPt[i]) ;
fc195fd0 1174 }
1175
6df33fcb 1176 fhEClusterSM = new TH2F("hEClusterSM","Raw clusters E and super-module number",
1177 nptbins,ptmin,ptmax,
1178 GetCaloUtils()->GetNumberOfSuperModulesUsed(),0,GetCaloUtils()->GetNumberOfSuperModulesUsed());
1179 fhEClusterSM->SetYTitle("SuperModule ");
184ca640 1180 fhEClusterSM->SetXTitle("#it{E} (GeV)");
6df33fcb 1181 outputContainer->Add(fhEClusterSM) ;
1182
11baad66 1183 fhPtClusterSM = new TH2F("hPtClusterSM","Raw clusters #it{p}_{T} and super-module number",
6df33fcb 1184 nptbins,ptmin,ptmax,
1185 GetCaloUtils()->GetNumberOfSuperModulesUsed(),0,GetCaloUtils()->GetNumberOfSuperModulesUsed());
1186 fhPtClusterSM->SetYTitle("SuperModule ");
184ca640 1187 fhPtClusterSM->SetXTitle("#it{E} (GeV)");
6df33fcb 1188 outputContainer->Add(fhPtClusterSM) ;
1189
1190 fhEPhotonSM = new TH2F("hEPhotonSM","Selected clusters E and super-module number",
1191 nptbins,ptmin,ptmax,
1192 GetCaloUtils()->GetNumberOfSuperModulesUsed(),0,GetCaloUtils()->GetNumberOfSuperModulesUsed());
1193 fhEPhotonSM->SetYTitle("SuperModule ");
184ca640 1194 fhEPhotonSM->SetXTitle("#it{E} (GeV)");
6df33fcb 1195 outputContainer->Add(fhEPhotonSM) ;
1196
11baad66 1197 fhPtPhotonSM = new TH2F("hPtPhotonSM","Selected clusters #it{p}_{T} and super-module number",
6df33fcb 1198 nptbins,ptmin,ptmax,
1199 GetCaloUtils()->GetNumberOfSuperModulesUsed(),0,GetCaloUtils()->GetNumberOfSuperModulesUsed());
1200 fhPtPhotonSM->SetYTitle("SuperModule ");
184ca640 1201 fhPtPhotonSM->SetXTitle("#it{E} (GeV)");
6df33fcb 1202 outputContainer->Add(fhPtPhotonSM) ;
1203
85c4406e 1204 fhNCellsE = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax);
184ca640 1205 fhNCellsE->SetXTitle("#it{E} (GeV)");
c4a7d28a 1206 fhNCellsE->SetYTitle("# of cells in cluster");
85c4406e 1207 outputContainer->Add(fhNCellsE);
f15c25da 1208
85c4406e 1209 fhCellsE = new TH2F ("hCellsE","energy of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nptbins*2,ptmin,ptmax);
184ca640 1210 fhCellsE->SetXTitle("#it{E}_{cluster} (GeV)");
1211 fhCellsE->SetYTitle("#it{E}_{cell} (GeV)");
85c4406e 1212 outputContainer->Add(fhCellsE);
5c46c992 1213
b2e375c7 1214 fhTimePt = new TH2F ("hTimePt","time of cluster vs pT of clusters", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
11baad66 1215 fhTimePt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
184ca640 1216 fhTimePt->SetYTitle("#it{time} (ns)");
b2e375c7 1217 outputContainer->Add(fhTimePt);
6175da48 1218
f66d95af 1219 fhMaxCellDiffClusterE = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
85c4406e 1220 nptbins,ptmin,ptmax, 500,0,1.);
184ca640 1221 fhMaxCellDiffClusterE->SetXTitle("#it{E}_{cluster} (GeV) ");
1222 fhMaxCellDiffClusterE->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
85c4406e 1223 outputContainer->Add(fhMaxCellDiffClusterE);
f66d95af 1224
85c4406e 1225 fhEPhoton = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
184ca640 1226 fhEPhoton->SetYTitle("#it{counts}");
1227 fhEPhoton->SetXTitle("#it{E}_{#gamma}(GeV)");
85c4406e 1228 outputContainer->Add(fhEPhoton) ;
20218aea 1229
11baad66 1230 fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs #it{p}_{T}",nptbins,ptmin,ptmax);
184ca640 1231 fhPtPhoton->SetYTitle("#it{counts}");
1232 fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/#it{c})");
85c4406e 1233 outputContainer->Add(fhPtPhoton) ;
1234
0cea6003 1235 if(IsHighMultiplicityAnalysisOn())
1236 {
1237 fhPtCentralityPhoton = new TH2F("hPtCentralityPhoton","centrality vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,100);
1238 fhPtCentralityPhoton->SetYTitle("Centrality");
1239 fhPtCentralityPhoton->SetXTitle("#it{p}_{T}(GeV/#it{c})");
1240 outputContainer->Add(fhPtCentralityPhoton) ;
1241
1242 fhPtEventPlanePhoton = new TH2F("hPtEventPlanePhoton","centrality vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
1243 fhPtEventPlanePhoton->SetYTitle("Event plane angle (rad)");
1244 fhPtEventPlanePhoton->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1245 outputContainer->Add(fhPtEventPlanePhoton) ;
1246 }
85c4406e 1247
c2a62a94 1248 fhEtaPhi = new TH2F
184ca640 1249 ("hEtaPhi","cluster,#it{E} > 0.5 GeV, #eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
c2a62a94 1250 fhEtaPhi->SetYTitle("#phi (rad)");
1251 fhEtaPhi->SetXTitle("#eta");
1252 outputContainer->Add(fhEtaPhi) ;
85c4406e 1253
477d6cee 1254 fhPhiPhoton = new TH2F
11baad66 1255 ("hPhiPhoton","#phi_{#gamma} vs #it{p}_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
6175da48 1256 fhPhiPhoton->SetYTitle("#phi (rad)");
184ca640 1257 fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/#it{c})");
85c4406e 1258 outputContainer->Add(fhPhiPhoton) ;
477d6cee 1259
1260 fhEtaPhoton = new TH2F
11baad66 1261 ("hEtaPhoton","#eta_{#gamma} vs #it{p}_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 1262 fhEtaPhoton->SetYTitle("#eta");
184ca640 1263 fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/#it{c})");
477d6cee 1264 outputContainer->Add(fhEtaPhoton) ;
1265
6175da48 1266 fhEtaPhiPhoton = new TH2F
85c4406e 1267 ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
6175da48 1268 fhEtaPhiPhoton->SetYTitle("#phi (rad)");
1269 fhEtaPhiPhoton->SetXTitle("#eta");
1270 outputContainer->Add(fhEtaPhiPhoton) ;
34c16486 1271 if(GetMinPt() < 0.5)
1272 {
20218aea 1273 fhEtaPhi05Photon = new TH2F
74e3eb22 1274 ("hEtaPhi05Photon","#eta vs #phi, E < 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
20218aea 1275 fhEtaPhi05Photon->SetYTitle("#phi (rad)");
1276 fhEtaPhi05Photon->SetXTitle("#eta");
1277 outputContainer->Add(fhEtaPhi05Photon) ;
1278 }
85c4406e 1279
9e51e29a 1280 fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
85c4406e 1281 nptbins,ptmin,ptmax,10,0,10);
9e51e29a 1282 fhNLocMax ->SetYTitle("N maxima");
184ca640 1283 fhNLocMax ->SetXTitle("#it{E} (GeV)");
85c4406e 1284 outputContainer->Add(fhNLocMax) ;
9e51e29a 1285
521636d2 1286 //Shower shape
34c16486 1287 if(fFillSSHistograms)
1288 {
85c4406e 1289 fhLam0E = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
521636d2 1290 fhLam0E->SetYTitle("#lambda_{0}^{2}");
184ca640 1291 fhLam0E->SetXTitle("#it{E} (GeV)");
85c4406e 1292 outputContainer->Add(fhLam0E);
521636d2 1293
85c4406e 1294 fhLam1E = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
521636d2 1295 fhLam1E->SetYTitle("#lambda_{1}^{2}");
184ca640 1296 fhLam1E->SetXTitle("#it{E} (GeV)");
85c4406e 1297 outputContainer->Add(fhLam1E);
521636d2 1298
85c4406e 1299 fhDispE = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
521636d2 1300 fhDispE->SetYTitle("D^{2}");
184ca640 1301 fhDispE->SetXTitle("#it{E} (GeV) ");
521636d2 1302 outputContainer->Add(fhDispE);
85c4406e 1303
b5dbb99b 1304 if(!fRejectTrackMatch)
1305 {
85c4406e 1306 fhLam0ETM = new TH2F ("hLam0ETM","#lambda_{0}^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
b5dbb99b 1307 fhLam0ETM->SetYTitle("#lambda_{0}^{2}");
184ca640 1308 fhLam0ETM->SetXTitle("#it{E} (GeV)");
85c4406e 1309 outputContainer->Add(fhLam0ETM);
b5dbb99b 1310
85c4406e 1311 fhLam1ETM = new TH2F ("hLam1ETM","#lambda_{1}^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
b5dbb99b 1312 fhLam1ETM->SetYTitle("#lambda_{1}^{2}");
184ca640 1313 fhLam1ETM->SetXTitle("#it{E} (GeV)");
85c4406e 1314 outputContainer->Add(fhLam1ETM);
b5dbb99b 1315
85c4406e 1316 fhDispETM = new TH2F ("hDispETM"," dispersion^{2} vs E, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
b5dbb99b 1317 fhDispETM->SetYTitle("D^{2}");
184ca640 1318 fhDispETM->SetXTitle("#it{E} (GeV) ");
b5dbb99b 1319 outputContainer->Add(fhDispETM);
1320 }
521636d2 1321
bdb39dbd 1322 if(GetCalorimeter() == kEMCAL && GetFirstSMCoveredByTRD() >= 0)
b5dbb99b 1323 {
85c4406e 1324 fhLam0ETRD = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
521636d2 1325 fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
184ca640 1326 fhLam0ETRD->SetXTitle("#it{E} (GeV)");
85c4406e 1327 outputContainer->Add(fhLam0ETRD);
521636d2 1328
85c4406e 1329 fhLam1ETRD = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
521636d2 1330 fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
184ca640 1331 fhLam1ETRD->SetXTitle("#it{E} (GeV)");
85c4406e 1332 outputContainer->Add(fhLam1ETRD);
521636d2 1333
85c4406e 1334 fhDispETRD = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
521636d2 1335 fhDispETRD->SetYTitle("Dispersion^{2}");
184ca640 1336 fhDispETRD->SetXTitle("#it{E} (GeV) ");
b5dbb99b 1337 outputContainer->Add(fhDispETRD);
1338
4d1d8f00 1339 if(!fRejectTrackMatch && GetFirstSMCoveredByTRD() >=0 )
b5dbb99b 1340 {
85c4406e 1341 fhLam0ETMTRD = new TH2F ("hLam0ETMTRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
b5dbb99b 1342 fhLam0ETMTRD->SetYTitle("#lambda_{0}^{2}");
184ca640 1343 fhLam0ETMTRD->SetXTitle("#it{E} (GeV)");
85c4406e 1344 outputContainer->Add(fhLam0ETMTRD);
b5dbb99b 1345
85c4406e 1346 fhLam1ETMTRD = new TH2F ("hLam1ETMTRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
b5dbb99b 1347 fhLam1ETMTRD->SetYTitle("#lambda_{1}^{2}");
184ca640 1348 fhLam1ETMTRD->SetXTitle("#it{E} (GeV)");
85c4406e 1349 outputContainer->Add(fhLam1ETMTRD);
b5dbb99b 1350
85c4406e 1351 fhDispETMTRD = new TH2F ("hDispETMTRD"," dispersion^{2} vs E, EMCAL SM covered by TRD, cut on track-matching residual |#Delta #eta| < 0.05, |#Delta #phi| < 0.05", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
b5dbb99b 1352 fhDispETMTRD->SetYTitle("Dispersion^{2}");
184ca640 1353 fhDispETMTRD->SetXTitle("#it{E} (GeV) ");
85c4406e 1354 outputContainer->Add(fhDispETMTRD);
1355 }
1356 }
521636d2 1357
764ab1f4 1358 if(!fFillOnlySimpleSSHisto)
34c16486 1359 {
85c4406e 1360 fhNCellsLam0LowE = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
764ab1f4 1361 fhNCellsLam0LowE->SetXTitle("N_{cells}");
1362 fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
85c4406e 1363 outputContainer->Add(fhNCellsLam0LowE);
764ab1f4 1364
184ca640 1365 fhNCellsLam0HighE = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, #it{E} > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
764ab1f4 1366 fhNCellsLam0HighE->SetXTitle("N_{cells}");
1367 fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
85c4406e 1368 outputContainer->Add(fhNCellsLam0HighE);
764ab1f4 1369
85c4406e 1370 fhNCellsLam1LowE = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
764ab1f4 1371 fhNCellsLam1LowE->SetXTitle("N_{cells}");
1372 fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
85c4406e 1373 outputContainer->Add(fhNCellsLam1LowE);
764ab1f4 1374
184ca640 1375 fhNCellsLam1HighE = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, #it{E} > 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
764ab1f4 1376 fhNCellsLam1HighE->SetXTitle("N_{cells}");
1377 fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
85c4406e 1378 outputContainer->Add(fhNCellsLam1HighE);
764ab1f4 1379
85c4406e 1380 fhNCellsDispLowE = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
764ab1f4 1381 fhNCellsDispLowE->SetXTitle("N_{cells}");
1382 fhNCellsDispLowE->SetYTitle("D^{2}");
85c4406e 1383 outputContainer->Add(fhNCellsDispLowE);
764ab1f4 1384
85c4406e 1385 fhNCellsDispHighE = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", nbins,nmin, nmax, ssbins,ssmin,ssmax);
764ab1f4 1386 fhNCellsDispHighE->SetXTitle("N_{cells}");
1387 fhNCellsDispHighE->SetYTitle("D^{2}");
85c4406e 1388 outputContainer->Add(fhNCellsDispHighE);
764ab1f4 1389
85c4406e 1390 fhEtaLam0LowE = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
764ab1f4 1391 fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
1392 fhEtaLam0LowE->SetXTitle("#eta");
85c4406e 1393 outputContainer->Add(fhEtaLam0LowE);
764ab1f4 1394
85c4406e 1395 fhPhiLam0LowE = new TH2F ("hPhiLam0LowE","#phi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
764ab1f4 1396 fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
1397 fhPhiLam0LowE->SetXTitle("#phi");
85c4406e 1398 outputContainer->Add(fhPhiLam0LowE);
764ab1f4 1399
184ca640 1400 fhEtaLam0HighE = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, #it{E} > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
764ab1f4 1401 fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
1402 fhEtaLam0HighE->SetXTitle("#eta");
85c4406e 1403 outputContainer->Add(fhEtaLam0HighE);
764ab1f4 1404
184ca640 1405 fhPhiLam0HighE = new TH2F ("hPhiLam0HighE","#phi vs #lambda_{0}^{2}, #it{E} > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
764ab1f4 1406 fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
1407 fhPhiLam0HighE->SetXTitle("#phi");
85c4406e 1408 outputContainer->Add(fhPhiLam0HighE);
764ab1f4 1409
85c4406e 1410 fhLam1Lam0LowE = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
764ab1f4 1411 fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
1412 fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
85c4406e 1413 outputContainer->Add(fhLam1Lam0LowE);
764ab1f4 1414
184ca640 1415 fhLam1Lam0HighE = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of #it{E} > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
764ab1f4 1416 fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
1417 fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
85c4406e 1418 outputContainer->Add(fhLam1Lam0HighE);
764ab1f4 1419
85c4406e 1420 fhLam0DispLowE = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
764ab1f4 1421 fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
1422 fhLam0DispLowE->SetYTitle("D^{2}");
85c4406e 1423 outputContainer->Add(fhLam0DispLowE);
764ab1f4 1424
184ca640 1425 fhLam0DispHighE = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of #it{E} > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
764ab1f4 1426 fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
1427 fhLam0DispHighE->SetYTitle("D^{2}");
85c4406e 1428 outputContainer->Add(fhLam0DispHighE);
764ab1f4 1429
85c4406e 1430 fhDispLam1LowE = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
764ab1f4 1431 fhDispLam1LowE->SetXTitle("D^{2}");
1432 fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
85c4406e 1433 outputContainer->Add(fhDispLam1LowE);
764ab1f4 1434
184ca640 1435 fhDispLam1HighE = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of #it{E} > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
764ab1f4 1436 fhDispLam1HighE->SetXTitle("D^{2}");
1437 fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
85c4406e 1438 outputContainer->Add(fhDispLam1HighE);
764ab1f4 1439
bdb39dbd 1440 if(GetCalorimeter() == kEMCAL)
34c16486 1441 {
85c4406e 1442 fhDispEtaE = new TH2F ("hDispEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
184ca640 1443 fhDispEtaE->SetXTitle("#it{E} (GeV)");
764ab1f4 1444 fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
85c4406e 1445 outputContainer->Add(fhDispEtaE);
764ab1f4 1446
85c4406e 1447 fhDispPhiE = new TH2F ("hDispPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
184ca640 1448 fhDispPhiE->SetXTitle("#it{E} (GeV)");
764ab1f4 1449 fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
85c4406e 1450 outputContainer->Add(fhDispPhiE);
764ab1f4 1451
85c4406e 1452 fhSumEtaE = new TH2F ("hSumEtaE","#delta^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
184ca640 1453 fhSumEtaE->SetXTitle("#it{E} (GeV)");
764ab1f4 1454 fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
85c4406e 1455 outputContainer->Add(fhSumEtaE);
764ab1f4 1456
85c4406e 1457 fhSumPhiE = new TH2F ("hSumPhiE","#delta^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
1458 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
184ca640 1459 fhSumPhiE->SetXTitle("#it{E} (GeV)");
764ab1f4 1460 fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
85c4406e 1461 outputContainer->Add(fhSumPhiE);
764ab1f4 1462
85c4406e 1463 fhSumEtaPhiE = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
1464 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
184ca640 1465 fhSumEtaPhiE->SetXTitle("#it{E} (GeV)");
764ab1f4 1466 fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
1467 outputContainer->Add(fhSumEtaPhiE);
1468
85c4406e 1469 fhDispEtaPhiDiffE = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
1470 nptbins,ptmin,ptmax,200, -10,10);
184ca640 1471 fhDispEtaPhiDiffE->SetXTitle("#it{E} (GeV)");
764ab1f4 1472 fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
85c4406e 1473 outputContainer->Add(fhDispEtaPhiDiffE);
bfdcf7fb 1474
85c4406e 1475 fhSphericityE = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
1476 nptbins,ptmin,ptmax, 200, -1,1);
184ca640 1477 fhSphericityE->SetXTitle("#it{E} (GeV)");
764ab1f4 1478 fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1479 outputContainer->Add(fhSphericityE);
bfdcf7fb 1480
85c4406e 1481 fhDispSumEtaDiffE = new TH2F ("hDispSumEtaDiffE","#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average vs E", nptbins,ptmin,ptmax, 200,-0.01,0.01);
184ca640 1482 fhDispSumEtaDiffE->SetXTitle("#it{E} (GeV)");
764ab1f4 1483 fhDispSumEtaDiffE->SetYTitle("#sigma^{2}_{#eta #eta} - #delta^{2}_{#eta #eta} / average");
85c4406e 1484 outputContainer->Add(fhDispSumEtaDiffE);
764ab1f4 1485
85c4406e 1486 fhDispSumPhiDiffE = new TH2F ("hDispSumPhiDiffE","#sigma^{2}_{#phi #phi} - #delta^{2}_{#phi #phi} / average vs E", nptbins,ptmin,ptmax, 200,-0.01,0.01);
184ca640 1487 fhDispSumPhiDiffE->SetXTitle("#it{E} (GeV)");
764ab1f4 1488 fhDispSumPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi} - #delta^{2}_{#phi #phi} / average");
85c4406e 1489 outputContainer->Add(fhDispSumPhiDiffE);
764ab1f4 1490
1491 for(Int_t i = 0; i < 7; i++)
1492 {
85c4406e 1493 fhDispEtaDispPhi[i] = new TH2F (Form("hDispEtaDispPhi_EBin%d",i),Form("#sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
1494 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
764ab1f4 1495 fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1496 fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
85c4406e 1497 outputContainer->Add(fhDispEtaDispPhi[i]);
764ab1f4 1498
85c4406e 1499 fhLambda0DispEta[i] = new TH2F (Form("hLambda0DispEta_EBin%d",i),Form("#lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
1500 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
764ab1f4 1501 fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
1502 fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
85c4406e 1503 outputContainer->Add(fhLambda0DispEta[i]);
764ab1f4 1504
85c4406e 1505 fhLambda0DispPhi[i] = new TH2F (Form("hLambda0DispPhi_EBin%d",i),Form("#lambda^{2}_{0}} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",bin[i],bin[i+1]),
1506 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
764ab1f4 1507 fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
1508 fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
85c4406e 1509 outputContainer->Add(fhLambda0DispPhi[i]);
764ab1f4 1510 }
34c16486 1511 }
1512 }
521636d2 1513 } // Shower shape
1514
09273901 1515 // Track Matching
1516
b5dbb99b 1517 if(fFillTMHisto)
1518 {
b2e375c7 1519 TString cutTM [] = {"NoCut",""};
b5dbb99b 1520
b2e375c7 1521 for(Int_t i = 0; i < 2; i++)
31ae6d59 1522 {
b2e375c7 1523 fhTrackMatchedDEta[i] = new TH2F
1524 (Form("hTrackMatchedDEta%s",cutTM[i].Data()),
1525 Form("d#eta of cluster-track vs cluster energy, %s",cutTM[i].Data()),
85c4406e 1526 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
b2e375c7 1527 fhTrackMatchedDEta[i]->SetYTitle("d#eta");
184ca640 1528 fhTrackMatchedDEta[i]->SetXTitle("#it{E}_{cluster} (GeV)");
4bfeae64 1529
b2e375c7 1530 fhTrackMatchedDPhi[i] = new TH2F
1531 (Form("hTrackMatchedDPhi%s",cutTM[i].Data()),
1532 Form("d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
85c4406e 1533 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
b2e375c7 1534 fhTrackMatchedDPhi[i]->SetYTitle("d#phi (rad)");
184ca640 1535 fhTrackMatchedDPhi[i]->SetXTitle("#it{E}_{cluster} (GeV)");
b2e375c7 1536
1537 fhTrackMatchedDEtaDPhi[i] = new TH2F
1538 (Form("hTrackMatchedDEtaDPhi%s",cutTM[i].Data()),
1539 Form("d#eta vs d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
1540 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1541 fhTrackMatchedDEtaDPhi[i]->SetYTitle("d#phi (rad)");
1542 fhTrackMatchedDEtaDPhi[i]->SetXTitle("d#eta");
1543
1544 fhTrackMatchedDEtaPos[i] = new TH2F
1545 (Form("hTrackMatchedDEtaPos%s",cutTM[i].Data()),
1546 Form("d#eta of cluster-track vs cluster energy, %s",cutTM[i].Data()),
85c4406e 1547 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
b2e375c7 1548 fhTrackMatchedDEtaPos[i]->SetYTitle("d#eta");
184ca640 1549 fhTrackMatchedDEtaPos[i]->SetXTitle("#it{E}_{cluster} (GeV)");
8d6b7f60 1550
b2e375c7 1551 fhTrackMatchedDPhiPos[i] = new TH2F
1552 (Form("hTrackMatchedDPhiPos%s",cutTM[i].Data()),
1553 Form("d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
85c4406e 1554 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
b2e375c7 1555 fhTrackMatchedDPhiPos[i]->SetYTitle("d#phi (rad)");
184ca640 1556 fhTrackMatchedDPhiPos[i]->SetXTitle("#it{E}_{cluster} (GeV)");
b2e375c7 1557
1558 fhTrackMatchedDEtaDPhiPos[i] = new TH2F
1559 (Form("hTrackMatchedDEtaDPhiPos%s",cutTM[i].Data()),
1560 Form("d#eta vs d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
1561 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1562 fhTrackMatchedDEtaDPhiPos[i]->SetYTitle("d#phi (rad)");
1563 fhTrackMatchedDEtaDPhiPos[i]->SetXTitle("d#eta");
1564
1565 fhTrackMatchedDEtaNeg[i] = new TH2F
1566 (Form("hTrackMatchedDEtaNeg%s",cutTM[i].Data()),
1567 Form("d#eta of cluster-track vs cluster energy, %s",cutTM[i].Data()),
85c4406e 1568 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
b2e375c7 1569 fhTrackMatchedDEtaNeg[i]->SetYTitle("d#eta");
184ca640 1570 fhTrackMatchedDEtaNeg[i]->SetXTitle("#it{E}_{cluster} (GeV)");
4bfeae64 1571
b2e375c7 1572 fhTrackMatchedDPhiNeg[i] = new TH2F
1573 (Form("hTrackMatchedDPhiNeg%s",cutTM[i].Data()),
1574 Form("d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
85c4406e 1575 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
b2e375c7 1576 fhTrackMatchedDPhiNeg[i]->SetYTitle("d#phi (rad)");
184ca640 1577 fhTrackMatchedDPhiNeg[i]->SetXTitle("#it{E}_{cluster} (GeV)");
b2e375c7 1578
1579 fhTrackMatchedDEtaDPhiNeg[i] = new TH2F
1580 (Form("hTrackMatchedDEtaDPhiNeg%s",cutTM[i].Data()),
1581 Form("d#eta vs d#phi of cluster-track vs cluster energy, %s",cutTM[i].Data()),
1582 nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1583 fhTrackMatchedDEtaDPhiNeg[i]->SetYTitle("d#phi (rad)");
1584 fhTrackMatchedDEtaDPhiNeg[i]->SetXTitle("d#eta");
1585
1586 fhdEdx[i] = new TH2F (Form("hdEdx%s",cutTM[i].Data()),Form("matched track <dE/dx> vs cluster E, %s",cutTM[i].Data()),
1587 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
184ca640 1588 fhdEdx[i]->SetXTitle("#it{E} (GeV)");
b2e375c7 1589 fhdEdx[i]->SetYTitle("<dE/dx>");
1590
1591 fhEOverP[i] = new TH2F (Form("hEOverP%s",cutTM[i].Data()),Form("matched track E/p vs cluster E, %s",cutTM[i].Data()),
1592 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
184ca640 1593 fhEOverP[i]->SetXTitle("#it{E} (GeV)");
b2e375c7 1594 fhEOverP[i]->SetYTitle("E/p");
1595
1596 outputContainer->Add(fhTrackMatchedDEta[i]) ;
1597 outputContainer->Add(fhTrackMatchedDPhi[i]) ;
1598 outputContainer->Add(fhTrackMatchedDEtaDPhi[i]) ;
1599 outputContainer->Add(fhTrackMatchedDEtaPos[i]) ;
1600 outputContainer->Add(fhTrackMatchedDPhiPos[i]) ;
1601 outputContainer->Add(fhTrackMatchedDEtaDPhiPos[i]) ;
1602 outputContainer->Add(fhTrackMatchedDEtaNeg[i]) ;
1603 outputContainer->Add(fhTrackMatchedDPhiNeg[i]) ;
1604 outputContainer->Add(fhTrackMatchedDEtaDPhiNeg[i]) ;
1605 outputContainer->Add(fhdEdx[i]);
1606 outputContainer->Add(fhEOverP[i]);
1607
bdb39dbd 1608 if(GetCalorimeter()==kEMCAL && GetFirstSMCoveredByTRD() >=0 )
b2e375c7 1609 {
1610 fhTrackMatchedDEtaTRD[i] = new TH2F
1611 (Form("hTrackMatchedDEtaTRD%s",cutTM[i].Data()),
1612 Form("d#eta of cluster-track vs cluster energy, SM behind TRD, %s",cutTM[i].Data()),
1613 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1614 fhTrackMatchedDEtaTRD[i]->SetYTitle("d#eta");
184ca640 1615 fhTrackMatchedDEtaTRD[i]->SetXTitle("#it{E}_{cluster} (GeV)");
b2e375c7 1616
1617 fhTrackMatchedDPhiTRD[i] = new TH2F
1618 (Form("hTrackMatchedDPhiTRD%s",cutTM[i].Data()),
1619 Form("d#phi of cluster-track vs cluster energy, SM behing TRD, %s",cutTM[i].Data()),
1620 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1621 fhTrackMatchedDPhiTRD[i]->SetYTitle("d#phi (rad)");
184ca640 1622 fhTrackMatchedDPhiTRD[i]->SetXTitle("#it{E}_{cluster} (GeV)");
b2e375c7 1623
1624 fhEOverPTRD[i] = new TH2F
1625 (Form("hEOverPTRD%s",cutTM[i].Data()),
1626 Form("matched track E/p vs cluster E, behind TRD, %s",cutTM[i].Data()),
1627 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
184ca640 1628 fhEOverPTRD[i]->SetXTitle("#it{E} (GeV)");
b2e375c7 1629 fhEOverPTRD[i]->SetYTitle("E/p");
1630
1631 outputContainer->Add(fhTrackMatchedDEtaTRD[i]) ;
1632 outputContainer->Add(fhTrackMatchedDPhiTRD[i]) ;
1633 outputContainer->Add(fhEOverPTRD[i]);
1634 }
8d6b7f60 1635
b2e375c7 1636 if(IsDataMC())
1637 {
1638 fhTrackMatchedDEtaMCNoOverlap[i] = new TH2F
1639 (Form("hTrackMatchedDEtaMCNoOverlap%s",cutTM[i].Data()),
1640 Form("d#eta of cluster-track vs cluster energy, no other MC particles overlap %s",cutTM[i].Data()),
1641 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1642 fhTrackMatchedDEtaMCNoOverlap[i]->SetYTitle("d#eta");
184ca640 1643 fhTrackMatchedDEtaMCNoOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
b2e375c7 1644
1645 fhTrackMatchedDPhiMCNoOverlap[i] = new TH2F
1646 (Form("hTrackMatchedDPhiMCNoOverlap%s",cutTM[i].Data()),
1647 Form("d#phi of cluster-track vs cluster energy, no other MC particles overlap %s",cutTM[i].Data()),
1648 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1649 fhTrackMatchedDPhiMCNoOverlap[i]->SetYTitle("d#phi (rad)");
184ca640 1650 fhTrackMatchedDPhiMCNoOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
b2e375c7 1651
1652 outputContainer->Add(fhTrackMatchedDEtaMCNoOverlap[i]) ;
1653 outputContainer->Add(fhTrackMatchedDPhiMCNoOverlap[i]) ;
1654 fhTrackMatchedDEtaMCOverlap[i] = new TH2F
1655 (Form("hTrackMatchedDEtaMCOverlap%s",cutTM[i].Data()),
1656 Form("d#eta of cluster-track vs cluster energy, several MC particles overlap %s",cutTM[i].Data()),
1657 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1658 fhTrackMatchedDEtaMCOverlap[i]->SetYTitle("d#eta");
184ca640 1659 fhTrackMatchedDEtaMCOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
b2e375c7 1660
1661 fhTrackMatchedDPhiMCOverlap[i] = new TH2F
1662 (Form("hTrackMatchedDPhiMCOverlap%s",cutTM[i].Data()),
1663 Form("d#phi of cluster-track vs cluster energy, several MC particles overlap %s",cutTM[i].Data()),
1664 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1665 fhTrackMatchedDPhiMCOverlap[i]->SetYTitle("d#phi (rad)");
184ca640 1666 fhTrackMatchedDPhiMCOverlap[i]->SetXTitle("#it{E}_{cluster} (GeV)");
b2e375c7 1667
1668 outputContainer->Add(fhTrackMatchedDEtaMCOverlap[i]) ;
1669 outputContainer->Add(fhTrackMatchedDPhiMCOverlap[i]) ;
1670
1671 fhTrackMatchedDEtaMCConversion[i] = new TH2F
1672 (Form("hTrackMatchedDEtaMCConversion%s",cutTM[i].Data()),
1673 Form("d#eta of cluster-track vs cluster energy, no other MC particles overlap appart from conversions %s",cutTM[i].Data()),
1674 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1675 fhTrackMatchedDEtaMCConversion[i]->SetYTitle("d#eta");
184ca640 1676 fhTrackMatchedDEtaMCConversion[i]->SetXTitle("#it{E}_{cluster} (GeV)");
b2e375c7 1677
1678 fhTrackMatchedDPhiMCConversion[i] = new TH2F
1679 (Form("hTrackMatchedDPhiMCConversion%s",cutTM[i].Data()),
1680 Form("d#phi of cluster-track vs cluster energy, no other MC particles overlap appart from conversions %s",cutTM[i].Data()),
1681 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1682 fhTrackMatchedDPhiMCConversion[i]->SetYTitle("d#phi (rad)");
184ca640 1683 fhTrackMatchedDPhiMCConversion[i]->SetXTitle("#it{E}_{cluster} (GeV)");
b2e375c7 1684
1685 outputContainer->Add(fhTrackMatchedDEtaMCConversion[i]) ;
1686 outputContainer->Add(fhTrackMatchedDPhiMCConversion[i]) ;
1687
1688 fhTrackMatchedMCParticle[i] = new TH2F
1689 (Form("hTrackMatchedMCParticle%s",cutTM[i].Data()),
1690 Form("Origin of particle vs energy %s",cutTM[i].Data()),
1691 nptbins,ptmin,ptmax,8,0,8);
184ca640 1692 fhTrackMatchedMCParticle[i]->SetXTitle("#it{E} (GeV)");
b2e375c7 1693 //fhTrackMatchedMCParticle[i]->SetYTitle("Particle type");
1694
1695 fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(1 ,"Photon");
1696 fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(2 ,"Electron");
1697 fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1698 fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(4 ,"Rest");
1699 fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1700 fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1701 fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1702 fhTrackMatchedMCParticle[i]->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1703
1704 outputContainer->Add(fhTrackMatchedMCParticle[i]);
1705 }
31ae6d59 1706 }
85c4406e 1707 }
09273901 1708
0cea6003 1709 if(IsPileUpAnalysisOn())
2ad19c3d 1710 {
5e5e056f 1711 TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1712
1713 for(Int_t i = 0 ; i < 7 ; i++)
1714 {
1715 fhPtPhotonPileUp[i] = new TH1F(Form("hPtPhotonPileUp%s",pileUpName[i].Data()),
11baad66 1716 Form("Selected photon #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1717 fhPtPhotonPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
5e5e056f 1718 outputContainer->Add(fhPtPhotonPileUp[i]);
fad96885 1719
fad96885 1720 fhClusterTimeDiffPhotonPileUp[i] = new TH2F(Form("hClusterTimeDiffPhotonPileUp%s",pileUpName[i].Data()),
bc41680b 1721 Form("Photon cluster E vs #it{t}_{max}-#it{t}_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
1722 nptbins,ptmin,ptmax,400,-200,200);
184ca640 1723 fhClusterTimeDiffPhotonPileUp[i]->SetXTitle("#it{E} (GeV)");
bc41680b 1724 fhClusterTimeDiffPhotonPileUp[i]->SetYTitle("#it{t}_{max}-#it{t}_{cell} (ns)");
fad96885 1725 outputContainer->Add(fhClusterTimeDiffPhotonPileUp[i]);
5e5e056f 1726 }
1727
b2e375c7 1728 fhTimePtPhotonNoCut = new TH2F ("hTimePtPhoton_NoCut","time of photon cluster vs pT of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
11baad66 1729 fhTimePtPhotonNoCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
184ca640 1730 fhTimePtPhotonNoCut->SetYTitle("#it{time} (ns)");
b2e375c7 1731 outputContainer->Add(fhTimePtPhotonNoCut);
1732
1733 fhTimePtPhotonSPD = new TH2F ("hTimePtPhoton_SPD","time of photon cluster vs pT of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
11baad66 1734 fhTimePtPhotonSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
184ca640 1735 fhTimePtPhotonSPD->SetYTitle("#it{time} (ns)");
b2e375c7 1736 outputContainer->Add(fhTimePtPhotonSPD);
bc41680b 1737
85c4406e 1738 fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,20,0,20);
2ad19c3d 1739 fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
184ca640 1740 fhTimeNPileUpVertSPD->SetXTitle("#it{time} (ns)");
fad96885 1741 outputContainer->Add(fhTimeNPileUpVertSPD);
2ad19c3d 1742
85c4406e 1743 fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 20,0,20 );
2ad19c3d 1744 fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
184ca640 1745 fhTimeNPileUpVertTrack->SetXTitle("#it{time} (ns)");
85c4406e 1746 outputContainer->Add(fhTimeNPileUpVertTrack);
2ad19c3d 1747
85c4406e 1748 fhPtPhotonNPileUpSPDVtx = new TH2F ("hPtPhoton_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
1749 nptbins,ptmin,ptmax,20,0,20);
1750 fhPtPhotonNPileUpSPDVtx->SetYTitle("# vertex ");
11baad66 1751 fhPtPhotonNPileUpSPDVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
85c4406e 1752 outputContainer->Add(fhPtPhotonNPileUpSPDVtx);
0f7e7205 1753
85c4406e 1754 fhPtPhotonNPileUpTrkVtx = new TH2F ("hPtPhoton_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
1755 nptbins,ptmin,ptmax, 20,0,20 );
1756 fhPtPhotonNPileUpTrkVtx->SetYTitle("# vertex ");
11baad66 1757 fhPtPhotonNPileUpTrkVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
85c4406e 1758 outputContainer->Add(fhPtPhotonNPileUpTrkVtx);
0f7e7205 1759
85c4406e 1760 fhPtPhotonNPileUpSPDVtxTimeCut = new TH2F ("hPtPhoton_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
1761 nptbins,ptmin,ptmax,20,0,20);
1762 fhPtPhotonNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
11baad66 1763 fhPtPhotonNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
85c4406e 1764 outputContainer->Add(fhPtPhotonNPileUpSPDVtxTimeCut);
0f7e7205 1765
85c4406e 1766 fhPtPhotonNPileUpTrkVtxTimeCut = new TH2F ("hPtPhoton_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
1767 nptbins,ptmin,ptmax, 20,0,20 );
1768 fhPtPhotonNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
11baad66 1769 fhPtPhotonNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
85c4406e 1770 outputContainer->Add(fhPtPhotonNPileUpTrkVtxTimeCut);
0f7e7205 1771
85c4406e 1772 fhPtPhotonNPileUpSPDVtxTimeCut2 = new TH2F ("hPtPhoton_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
1773 nptbins,ptmin,ptmax,20,0,20);
1774 fhPtPhotonNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
11baad66 1775 fhPtPhotonNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
85c4406e 1776 outputContainer->Add(fhPtPhotonNPileUpSPDVtxTimeCut2);
0f7e7205 1777
85c4406e 1778 fhPtPhotonNPileUpTrkVtxTimeCut2 = new TH2F ("hPtPhoton_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
1779 nptbins,ptmin,ptmax, 20,0,20 );
1780 fhPtPhotonNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
11baad66 1781 fhPtPhotonNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
85c4406e 1782 outputContainer->Add(fhPtPhotonNPileUpTrkVtxTimeCut2);
1783
2ad19c3d 1784 }
bc41680b 1785
1786
2ad19c3d 1787
34c16486 1788 if(IsDataMC())
1789 {
57cfa545 1790 TString ptype[] = { "#gamma" , "#gamma_{#pi decay}" , "#gamma_{#eta decay}", "#gamma_{other decay}",
1791 "#pi^{0}" , "#eta" , "e^{#pm}" , "#gamma->e^{#pm}" ,
1792 "hadron?" , "Anti-N" , "Anti-P" ,
1793 "#gamma_{prompt}", "#gamma_{fragmentation}", "#gamma_{ISR}" , "String" } ;
3d5d5078 1794
57cfa545 1795 TString pname[] = { "Photon" , "PhotonPi0Decay" , "PhotonEtaDecay", "PhotonOtherDecay",
1796 "Pi0" , "Eta" , "Electron" , "Conversion" ,
1797 "Hadron" , "AntiNeutron" , "AntiProton" ,
1798 "PhotonPrompt", "PhotonFragmentation", "PhotonISR" , "String" } ;
521636d2 1799
34c16486 1800 for(Int_t i = 0; i < fNOriginHistograms; i++)
85c4406e 1801 {
3d5d5078 1802 fhMCE[i] = new TH1F(Form("hE_MC%s",pname[i].Data()),
85c4406e 1803 Form("cluster from %s : E ",ptype[i].Data()),
1804 nptbins,ptmin,ptmax);
184ca640 1805 fhMCE[i]->SetXTitle("#it{E} (GeV)");
85c4406e 1806 outputContainer->Add(fhMCE[i]) ;
521636d2 1807
4c8f7c2e 1808 fhMCPt[i] = new TH1F(Form("hPt_MC%s",pname[i].Data()),
11baad66 1809 Form("cluster from %s : #it{p}_{T} ",ptype[i].Data()),
85c4406e 1810 nptbins,ptmin,ptmax);
11baad66 1811 fhMCPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
4c8f7c2e 1812 outputContainer->Add(fhMCPt[i]) ;
521636d2 1813
4c8f7c2e 1814 fhMCEta[i] = new TH2F(Form("hEta_MC%s",pname[i].Data()),
85c4406e 1815 Form("cluster from %s : #eta ",ptype[i].Data()),
1816 nptbins,ptmin,ptmax,netabins,etamin,etamax);
4c8f7c2e 1817 fhMCEta[i]->SetYTitle("#eta");
184ca640 1818 fhMCEta[i]->SetXTitle("#it{E} (GeV)");
4c8f7c2e 1819 outputContainer->Add(fhMCEta[i]) ;
521636d2 1820
4c8f7c2e 1821 fhMCPhi[i] = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
85c4406e 1822 Form("cluster from %s : #phi ",ptype[i].Data()),
1823 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
4c8f7c2e 1824 fhMCPhi[i]->SetYTitle("#phi (rad)");
184ca640 1825 fhMCPhi[i]->SetXTitle("#it{E} (GeV)");
4c8f7c2e 1826 outputContainer->Add(fhMCPhi[i]) ;
1827
1828
d9105d92 1829 fhMCDeltaE[i] = new TH2F (Form("hDeltaE_MC%s",pname[i].Data()),
85c4406e 1830 Form("MC - Reco E from %s",pname[i].Data()),
1831 nptbins,ptmin,ptmax, 200,-50,50);
184ca640 1832 fhMCDeltaE[i]->SetYTitle("#Delta #it{E} (GeV)");
1833 fhMCDeltaE[i]->SetXTitle("#it{E} (GeV)");
4c8f7c2e 1834 outputContainer->Add(fhMCDeltaE[i]);
1835
d9105d92 1836 fhMCDeltaPt[i] = new TH2F (Form("hDeltaPt_MC%s",pname[i].Data()),
11baad66 1837 Form("MC - Reco #it{p}_{T} from %s",pname[i].Data()),
85c4406e 1838 nptbins,ptmin,ptmax, 200,-50,50);
184ca640 1839 fhMCDeltaPt[i]->SetXTitle("p_{T,rec} (GeV/#it{c})");
11baad66 1840 fhMCDeltaPt[i]->SetYTitle("#Delta #it{p}_{T} (GeV/#it{c})");
4c8f7c2e 1841 outputContainer->Add(fhMCDeltaPt[i]);
85c4406e 1842
4c8f7c2e 1843 fhMC2E[i] = new TH2F (Form("h2E_MC%s",pname[i].Data()),
85c4406e 1844 Form("E distribution, reconstructed vs generated from %s",pname[i].Data()),
1845 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
184ca640 1846 fhMC2E[i]->SetXTitle("#it{E}_{rec} (GeV)");
1847 fhMC2E[i]->SetYTitle("#it{E}_{gen} (GeV)");
85c4406e 1848 outputContainer->Add(fhMC2E[i]);
4c8f7c2e 1849
1850 fhMC2Pt[i] = new TH2F (Form("h2Pt_MC%s",pname[i].Data()),
85c4406e 1851 Form("p_T distribution, reconstructed vs generated from %s",pname[i].Data()),
1852 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
184ca640 1853 fhMC2Pt[i]->SetXTitle("p_{T,rec} (GeV/#it{c})");
1854 fhMC2Pt[i]->SetYTitle("p_{T,gen} (GeV/#it{c})");
4c8f7c2e 1855 outputContainer->Add(fhMC2Pt[i]);
1856
521636d2 1857
1858 }
3d5d5078 1859
57cfa545 1860 TString pptype[] = { "#gamma" , "#gamma_{#pi decay}" ,
1861 "#gamma_{#eta decay}", "#gamma_{other decay}" ,
1862 "#gamma_{prompt}" , "#gamma_{fragmentation}", "#gamma_{ISR}" } ;
f66d95af 1863
57cfa545 1864 TString ppname[] = { "Photon" , "PhotonPi0Decay" ,
1865 "PhotonEtaDecay", "PhotonOtherDecay" ,
1866 "PhotonPrompt" , "PhotonFragmentation", "PhotonISR" } ;
f66d95af 1867
34c16486 1868 for(Int_t i = 0; i < fNPrimaryHistograms; i++)
85c4406e 1869 {
f66d95af 1870 fhEPrimMC[i] = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
85c4406e 1871 Form("primary photon %s : E ",pptype[i].Data()),
1872 nptbins,ptmin,ptmax);
184ca640 1873 fhEPrimMC[i]->SetXTitle("#it{E} (GeV)");
85c4406e 1874 outputContainer->Add(fhEPrimMC[i]) ;
3d5d5078 1875
f66d95af 1876 fhPtPrimMC[i] = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
11baad66 1877 Form("primary photon %s : #it{p}_{T} ",pptype[i].Data()),
85c4406e 1878 nptbins,ptmin,ptmax);
11baad66 1879 fhPtPrimMC[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3d5d5078 1880 outputContainer->Add(fhPtPrimMC[i]) ;
1881
f66d95af 1882 fhYPrimMC[i] = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
85c4406e 1883 Form("primary photon %s : Rapidity ",pptype[i].Data()),
a421d7a3 1884 nptbins,ptmin,ptmax,200,-2,2);
3d5d5078 1885 fhYPrimMC[i]->SetYTitle("Rapidity");
184ca640 1886 fhYPrimMC[i]->SetXTitle("#it{E} (GeV)");
3d5d5078 1887 outputContainer->Add(fhYPrimMC[i]) ;
1888
4cf13296 1889 fhEtaPrimMC[i] = new TH2F(Form("hEtaPrim_MC%s",ppname[i].Data()),
1890 Form("primary photon %s : #eta",pptype[i].Data()),
a421d7a3 1891 nptbins,ptmin,ptmax,200,-2,2);
4cf13296 1892 fhEtaPrimMC[i]->SetYTitle("#eta");
184ca640 1893 fhEtaPrimMC[i]->SetXTitle("#it{E} (GeV)");
4cf13296 1894 outputContainer->Add(fhEtaPrimMC[i]) ;
1895
f66d95af 1896 fhPhiPrimMC[i] = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
85c4406e 1897 Form("primary photon %s : #phi ",pptype[i].Data()),
a421d7a3 1898 nptbins,ptmin,ptmax,nphibins,0,TMath::TwoPi());
3d5d5078 1899 fhPhiPrimMC[i]->SetYTitle("#phi (rad)");
184ca640 1900 fhPhiPrimMC[i]->SetXTitle("#it{E} (GeV)");
3d5d5078 1901 outputContainer->Add(fhPhiPrimMC[i]) ;
85c4406e 1902
3d5d5078 1903
f66d95af 1904 fhEPrimMCAcc[i] = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
85c4406e 1905 Form("primary photon %s in acceptance: E ",pptype[i].Data()),
1906 nptbins,ptmin,ptmax);
184ca640 1907 fhEPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
85c4406e 1908 outputContainer->Add(fhEPrimMCAcc[i]) ;
3d5d5078 1909
f66d95af 1910 fhPtPrimMCAcc[i] = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
11baad66 1911 Form("primary photon %s in acceptance: #it{p}_{T} ",pptype[i].Data()),
85c4406e 1912 nptbins,ptmin,ptmax);
11baad66 1913 fhPtPrimMCAcc[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3d5d5078 1914 outputContainer->Add(fhPtPrimMCAcc[i]) ;
1915
f66d95af 1916 fhYPrimMCAcc[i] = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
85c4406e 1917 Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
1918 nptbins,ptmin,ptmax,100,-1,1);
3d5d5078 1919 fhYPrimMCAcc[i]->SetYTitle("Rapidity");
184ca640 1920 fhYPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3d5d5078 1921 outputContainer->Add(fhYPrimMCAcc[i]) ;
4cf13296 1922
1923 fhEtaPrimMCAcc[i] = new TH2F(Form("hEtaPrimAcc_MC%s",ppname[i].Data()),
1924 Form("primary photon %s in acceptance: #eta ",pptype[i].Data()),
1925 nptbins,ptmin,ptmax,netabins,etamin,etamax);
1926 fhEtaPrimMCAcc[i]->SetYTitle("#eta");
184ca640 1927 fhEtaPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
667a3592 1928 outputContainer->Add(fhEtaPrimMCAcc[i]) ;
3d5d5078 1929
f66d95af 1930 fhPhiPrimMCAcc[i] = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
85c4406e 1931 Form("primary photon %s in acceptance: #phi ",pptype[i].Data()),
1932 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
3d5d5078 1933 fhPhiPrimMCAcc[i]->SetYTitle("#phi (rad)");
184ca640 1934 fhPhiPrimMCAcc[i]->SetXTitle("#it{E} (GeV)");
3d5d5078 1935 outputContainer->Add(fhPhiPrimMCAcc[i]) ;
1936
1937 }
85c4406e 1938
34c16486 1939 if(fFillSSHistograms)
1940 {
85c4406e 1941 TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
3d5d5078 1942
1943 TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
1944
34c16486 1945 for(Int_t i = 0; i < 6; i++)
85c4406e 1946 {
3d5d5078 1947 fhMCELambda0[i] = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
1948 Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
85c4406e 1949 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1950 fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
184ca640 1951 fhMCELambda0[i]->SetXTitle("#it{E} (GeV)");
85c4406e 1952 outputContainer->Add(fhMCELambda0[i]) ;
521636d2 1953
3d5d5078 1954 fhMCELambda1[i] = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
1955 Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
85c4406e 1956 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1957 fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
184ca640 1958 fhMCELambda1[i]->SetXTitle("#it{E} (GeV)");
85c4406e 1959 outputContainer->Add(fhMCELambda1[i]) ;
34c16486 1960
3d5d5078 1961 fhMCEDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
1962 Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
85c4406e 1963 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1964 fhMCEDispersion[i]->SetYTitle("D^{2}");
184ca640 1965 fhMCEDispersion[i]->SetXTitle("#it{E} (GeV)");
85c4406e 1966 outputContainer->Add(fhMCEDispersion[i]) ;
34c16486 1967
f66d95af 1968 fhMCNCellsE[i] = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
85c4406e 1969 Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()),
1970 nptbins,ptmin,ptmax, nbins,nmin,nmax);
184ca640 1971 fhMCNCellsE[i]->SetXTitle("#it{E} (GeV)");
f66d95af 1972 fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
85c4406e 1973 outputContainer->Add(fhMCNCellsE[i]);
f66d95af 1974
1975 fhMCMaxCellDiffClusterE[i] = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
34c16486 1976 Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
85c4406e 1977 nptbins,ptmin,ptmax, 500,0,1.);
184ca640 1978 fhMCMaxCellDiffClusterE[i]->SetXTitle("#it{E}_{cluster} (GeV) ");
1979 fhMCMaxCellDiffClusterE[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
85c4406e 1980 outputContainer->Add(fhMCMaxCellDiffClusterE[i]);
f66d95af 1981
764ab1f4 1982 if(!fFillOnlySimpleSSHisto)
34c16486 1983 {
764ab1f4 1984 fhMCLambda0vsClusterMaxCellDiffE0[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
1985 Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
85c4406e 1986 ssbins,ssmin,ssmax,500,0,1.);
764ab1f4 1987 fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
184ca640 1988 fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
85c4406e 1989 outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ;
764ab1f4 1990
1991 fhMCLambda0vsClusterMaxCellDiffE2[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
1992 Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
85c4406e 1993 ssbins,ssmin,ssmax,500,0,1.);
764ab1f4 1994 fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
184ca640 1995 fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
85c4406e 1996 outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ;
34c16486 1997
764ab1f4 1998 fhMCLambda0vsClusterMaxCellDiffE6[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
184ca640 1999 Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, #it{E} > 6 GeV",ptypess[i].Data()),
85c4406e 2000 ssbins,ssmin,ssmax,500,0,1.);
764ab1f4 2001 fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
184ca640 2002 fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
85c4406e 2003 outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ;
34c16486 2004
764ab1f4 2005 fhMCNCellsvsClusterMaxCellDiffE0[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
2006 Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
85c4406e 2007 nbins/5,nmin,nmax/5,500,0,1.);
764ab1f4 2008 fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
184ca640 2009 fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
85c4406e 2010 outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ;
34c16486 2011
764ab1f4 2012 fhMCNCellsvsClusterMaxCellDiffE2[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
2013 Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
85c4406e 2014 nbins/5,nmin,nmax/5,500,0,1.);
764ab1f4 2015 fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
184ca640 2016 fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
85c4406e 2017 outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ;
34c16486 2018
764ab1f4 2019 fhMCNCellsvsClusterMaxCellDiffE6[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
184ca640 2020 Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, #it{E} > 6 GeV",ptypess[i].Data()),
85c4406e 2021 nbins/5,nmin,nmax/5,500,0,1.);
764ab1f4 2022 fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
184ca640 2023 fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("#it{E} (GeV)");
85c4406e 2024 outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ;
34c16486 2025
bdb39dbd 2026 if(GetCalorimeter()==kEMCAL)
34c16486 2027 {
764ab1f4 2028 fhMCEDispEta[i] = new TH2F (Form("hEDispEtaE_MC%s",pnamess[i].Data()),
2029 Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data()),
85c4406e 2030 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
184ca640 2031 fhMCEDispEta[i]->SetXTitle("#it{E} (GeV)");
764ab1f4 2032 fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
85c4406e 2033 outputContainer->Add(fhMCEDispEta[i]);
764ab1f4 2034
2035 fhMCEDispPhi[i] = new TH2F (Form("hEDispPhiE_MC%s",pnamess[i].Data()),
2036 Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data()),
85c4406e 2037 nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
184ca640 2038 fhMCEDispPhi[i]->SetXTitle("#it{E} (GeV)");
764ab1f4 2039 fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
85c4406e 2040 outputContainer->Add(fhMCEDispPhi[i]);
764ab1f4 2041
2042 fhMCESumEtaPhi[i] = new TH2F (Form("hESumEtaPhiE_MC%s",pnamess[i].Data()),
85c4406e 2043 Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptypess[i].Data()),
2044 nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
184ca640 2045 fhMCESumEtaPhi[i]->SetXTitle("#it{E} (GeV)");
764ab1f4 2046 fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
2047 outputContainer->Add(fhMCESumEtaPhi[i]);
2048
2049 fhMCEDispEtaPhiDiff[i] = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pnamess[i].Data()),
85c4406e 2050 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data()),
2051 nptbins,ptmin,ptmax,200,-10,10);
184ca640 2052 fhMCEDispEtaPhiDiff[i]->SetXTitle("#it{E} (GeV)");
764ab1f4 2053 fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
85c4406e 2054 outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
764ab1f4 2055
2056 fhMCESphericity[i] = new TH2F (Form("hESphericity_MC%s",pnamess[i].Data()),
85c4406e 2057 Form("cluster from %s : (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptypess[i].Data()),
2058 nptbins,ptmin,ptmax, 200,-1,1);
184ca640 2059 fhMCESphericity[i]->SetXTitle("#it{E} (GeV)");
764ab1f4 2060 fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
2061 outputContainer->Add(fhMCESphericity[i]);
2062
2063 for(Int_t ie = 0; ie < 7; ie++)
2064 {
2065 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
85c4406e 2066 Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pnamess[i].Data(),bin[ie],bin[ie+1]),
2067 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
764ab1f4 2068 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
2069 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
85c4406e 2070 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
764ab1f4 2071
2072 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pnamess[i].Data()),
85c4406e 2073 Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pnamess[i].Data(),bin[ie],bin[ie+1]),
2074 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
764ab1f4 2075 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
2076 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
85c4406e 2077 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
764ab1f4 2078
2079 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pnamess[i].Data()),
85c4406e 2080 Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",pnamess[i].Data(),bin[ie],bin[ie+1]),
2081 ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
764ab1f4 2082 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
2083 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
85c4406e 2084 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
764ab1f4 2085 }
34c16486 2086 }
34c16486 2087 }
85c4406e 2088 }// loop
3d5d5078 2089
2090 if(!GetReader()->IsEmbeddedClusterSelectionOn())
2091 {
2092 fhMCPhotonELambda0NoOverlap = new TH2F("hELambda0_MCPhoton_NoOverlap",
2093 "cluster from Photon : E vs #lambda_{0}^{2}",
85c4406e 2094 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2095 fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
184ca640 2096 fhMCPhotonELambda0NoOverlap->SetXTitle("#it{E} (GeV)");
85c4406e 2097 outputContainer->Add(fhMCPhotonELambda0NoOverlap) ;
3d5d5078 2098
3d5d5078 2099 fhMCPhotonELambda0TwoOverlap = new TH2F("hELambda0_MCPhoton_TwoOverlap",
2100 "cluster from Photon : E vs #lambda_{0}^{2}",
85c4406e 2101 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2102 fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
184ca640 2103 fhMCPhotonELambda0TwoOverlap->SetXTitle("#it{E} (GeV)");
85c4406e 2104 outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ;
3d5d5078 2105
3d5d5078 2106 fhMCPhotonELambda0NOverlap = new TH2F("hELambda0_MCPhoton_NOverlap",
2107 "cluster from Photon : E vs #lambda_{0}^{2}",
85c4406e 2108 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2109 fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
184ca640 2110 fhMCPhotonELambda0NOverlap->SetXTitle("#it{E} (GeV)");
85c4406e 2111 outputContainer->Add(fhMCPhotonELambda0NOverlap) ;
521636d2 2112
85c4406e 2113 } //No embedding
3d5d5078 2114
3d5d5078 2115 if(GetReader()->IsEmbeddedClusterSelectionOn())
2116 {
2117
2118 fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
34c16486 2119 "Energy Fraction of embedded signal versus cluster energy",
85c4406e 2120 nptbins,ptmin,ptmax,100,0.,1.);
3d5d5078 2121 fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
184ca640 2122 fhEmbeddedSignalFractionEnergy->SetXTitle("#it{E} (GeV)");
85c4406e 2123 outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
3d5d5078 2124
2125 fhEmbedPhotonELambda0FullSignal = new TH2F("hELambda0_EmbedPhoton_FullSignal",
34c16486 2126 "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
85c4406e 2127 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2128 fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
184ca640 2129 fhEmbedPhotonELambda0FullSignal->SetXTitle("#it{E} (GeV)");
85c4406e 2130 outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ;
34c16486 2131
3d5d5078 2132 fhEmbedPhotonELambda0MostlySignal = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
34c16486 2133 "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
85c4406e 2134 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2135 fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
184ca640 2136 fhEmbedPhotonELambda0MostlySignal->SetXTitle("#it{E} (GeV)");
85c4406e 2137 outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ;
3d5d5078 2138
3d5d5078 2139 fhEmbedPhotonELambda0MostlyBkg = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
34c16486 2140 "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
85c4406e 2141 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2142 fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
184ca640 2143 fhEmbedPhotonELambda0MostlyBkg->SetXTitle("#it{E} (GeV)");
85c4406e 2144 outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ;
34c16486 2145
3d5d5078 2146 fhEmbedPhotonELambda0FullBkg = new TH2F("hELambda0_EmbedPhoton_FullBkg",
34c16486 2147 "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
85c4406e 2148 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2149 fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
184ca640 2150 fhEmbedPhotonELambda0FullBkg->SetXTitle("#it{E} (GeV)");
85c4406e 2151 outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ;
3d5d5078 2152
3d5d5078 2153 fhEmbedPi0ELambda0FullSignal = new TH2F("hELambda0_EmbedPi0_FullSignal",
34c16486 2154 "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
85c4406e 2155 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2156 fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
184ca640 2157 fhEmbedPi0ELambda0FullSignal->SetXTitle("#it{E} (GeV)");
85c4406e 2158 outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ;
34c16486 2159
3d5d5078 2160 fhEmbedPi0ELambda0MostlySignal = new TH2F("hELambda0_EmbedPi0_MostlySignal",
34c16486 2161 "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
85c4406e 2162 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2163 fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
184ca640 2164 fhEmbedPi0ELambda0MostlySignal->SetXTitle("#it{E} (GeV)");
85c4406e 2165 outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ;
3d5d5078 2166
3d5d5078 2167 fhEmbedPi0ELambda0MostlyBkg = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
34c16486 2168 "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
85c4406e 2169 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2170 fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
184ca640 2171 fhEmbedPi0ELambda0MostlyBkg->SetXTitle("#it{E} (GeV)");
85c4406e 2172 outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ;
3d5d5078 2173
3d5d5078 2174 fhEmbedPi0ELambda0FullBkg = new TH2F("hELambda0_EmbedPi0_FullBkg",
34c16486 2175 "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
85c4406e 2176 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 2177 fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
184ca640 2178 fhEmbedPi0ELambda0FullBkg->SetXTitle("#it{E} (GeV)");
85c4406e 2179 outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ;
34c16486 2180
3d5d5078 2181 }// embedded histograms
2182
521636d2 2183
2184 }// Fill SS MC histograms
2185
477d6cee 2186 }//Histos with MC
1035a8d9 2187
477d6cee 2188 return outputContainer ;
2189
1c5acb87 2190}
2191
34c16486 2192//_______________________
6639984f 2193void AliAnaPhoton::Init()
2194{
6639984f 2195 //Init
2db10729 2196
6639984f 2197 //Do some checks
2db10729 2198
2199 if ( GetCalorimeter() == kPHOS && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD() )
2200 AliFatal("!!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!");
2201 else if( GetCalorimeter() == kEMCAL && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD() )
2202 AliFatal("!!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!");
6639984f 2203
49b5c49b 2204 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC) GetCaloPID()->SwitchOnBayesian();
2205
6639984f 2206}
2207
1c5acb87 2208//____________________________________________________________________________
2209void AliAnaPhoton::InitParameters()
2210{
2211
2212 //Initialize the parameters of the analysis.
a3aebfff 2213 AddToHistogramsName("AnaPhoton_");
521636d2 2214
6175da48 2215 fMinDist = 2.;
2216 fMinDist2 = 4.;
2217 fMinDist3 = 5.;
1e86c71e 2218
caa8a222 2219 fTimeCutMin =-1000000;
2220 fTimeCutMax = 1000000;
6175da48 2221 fNCellsCut = 0;
2ac125bf 2222
1e86c71e 2223 fRejectTrackMatch = kTRUE ;
2e3e315d 2224
1c5acb87 2225}
2226
2227//__________________________________________________________________
85c4406e 2228void AliAnaPhoton::MakeAnalysisFillAOD()
1c5acb87 2229{
f8006433 2230 //Do photon analysis and fill aods
f37fa8d2 2231
85c4406e 2232 //Get the vertex
5025c139 2233 Double_t v[3] = {0,0,0}; //vertex ;
2234 GetReader()->GetVertex(v);
f8006433 2235
f37fa8d2 2236 //Select the Calorimeter of the photon
85c4406e 2237 TObjArray * pl = 0x0;
2238 AliVCaloCells* cells = 0;
bdb39dbd 2239 if (GetCalorimeter() == kPHOS )
71e3889f 2240 {
2241 pl = GetPHOSClusters();
2242 cells = GetPHOSCells();
2243 }
bdb39dbd 2244 else if (GetCalorimeter() == kEMCAL)
71e3889f 2245 {
2246 pl = GetEMCALClusters();
2247 cells = GetEMCALCells();
2248 }
5ae09196 2249
85c4406e 2250 if(!pl)
34c16486 2251 {
2db10729 2252 AliWarning(Form("TObjArray with %s clusters is NULL!",GetCalorimeterString().Data()));
5ae09196 2253 return;
2254 }
521636d2 2255
fc195fd0 2256 // Loop on raw clusters before filtering in the reader and fill control histogram
bdb39dbd 2257 if((GetReader()->GetEMCALClusterListName()=="" && GetCalorimeter()==kEMCAL) || GetCalorimeter()==kPHOS)
34c16486 2258 {
2259 for(Int_t iclus = 0; iclus < GetReader()->GetInputEvent()->GetNumberOfCaloClusters(); iclus++ )
2260 {
fc195fd0 2261 AliVCluster * clus = GetReader()->GetInputEvent()->GetCaloCluster(iclus);
bdb39dbd 2262 if (GetCalorimeter() == kPHOS && clus->IsPHOS() && clus->E() > GetReader()->GetPHOSPtMin() )
58ea8ce5 2263 {
2264 fhClusterCutsE [0]->Fill(clus->E());
2265
da13b2b8 2266 clus->GetMomentum(fMomentum,GetVertex(0)) ;
2267 fhClusterCutsPt[0]->Fill(fMomentum.Pt());
58ea8ce5 2268 }
bdb39dbd 2269 else if(GetCalorimeter() == kEMCAL && clus->IsEMCAL() && clus->E() > GetReader()->GetEMCALPtMin())
58ea8ce5 2270 {
2271 fhClusterCutsE [0]->Fill(clus->E());
2272
da13b2b8 2273 clus->GetMomentum(fMomentum,GetVertex(0)) ;
2274 fhClusterCutsPt[0]->Fill(fMomentum.Pt());
58ea8ce5 2275 }
fc195fd0 2276 }
2277 }
85c4406e 2278 else
34c16486 2279 { // reclusterized
fc195fd0 2280 TClonesArray * clusterList = 0;
85c4406e 2281
7d650cb7 2282 if(GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()))
2283 clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetInputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
2284 else if(GetReader()->GetOutputEvent())
4a9e1073 2285 clusterList = dynamic_cast<TClonesArray*> (GetReader()->GetOutputEvent()->FindListObject(GetReader()->GetEMCALClusterListName()));
7d650cb7 2286
34c16486 2287 if(clusterList)
2288 {
fc195fd0 2289 Int_t nclusters = clusterList->GetEntriesFast();
85c4406e 2290 for (Int_t iclus = 0; iclus < nclusters; iclus++)
34c16486 2291 {
85c4406e 2292 AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
58ea8ce5 2293 if(clus)
2294 {
2295 fhClusterCutsE [0]->Fill(clus->E());
2296
da13b2b8 2297 clus->GetMomentum(fMomentum,GetVertex(0)) ;
2298 fhClusterCutsPt[0]->Fill(fMomentum.Pt());
58ea8ce5 2299 }
6265ad55 2300 }
fc195fd0 2301 }
2302 }
fc195fd0 2303
da13b2b8 2304 // Init arrays, variables, get number of clusters
1e86c71e 2305 Int_t nCaloClusters = pl->GetEntriesFast();
20218aea 2306
2db10729 2307 AliDebug(1,Form("Input %s cluster entries %d", GetCalorimeterString().Data(), nCaloClusters));
521636d2 2308
6175da48 2309 //----------------------------------------------------
2310 // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
2311 //----------------------------------------------------
2312 // Loop on clusters
34c16486 2313 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
85c4406e 2314 {
2315 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
0ae57829 2316 //printf("calo %d, %f\n",icalo,calo->E());
521636d2 2317
f8006433 2318 //Get the index where the cluster comes, to retrieve the corresponding vertex
85c4406e 2319 Int_t evtIndex = 0 ;
2320 if (GetMixedEvent())
34c16486 2321 {
85c4406e 2322 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
5025c139 2323 //Get the vertex and check it is not too large in z
96539743 2324 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
c8fe2783 2325 }
521636d2 2326
85c4406e 2327 //Cluster selection, not charged, with photon id and in fiducial cut
34c16486 2328 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2329 {
da13b2b8 2330 calo->GetMomentum(fMomentum,GetVertex(evtIndex)) ;
afb3af8a 2331 }//Assume that come from vertex in straight line
34c16486 2332 else
2333 {
f8006433 2334 Double_t vertex[]={0,0,0};
da13b2b8 2335 calo->GetMomentum(fMomentum,vertex) ;
f8006433 2336 }
c2a62a94 2337
bc41680b 2338 //-----------------------------
6175da48 2339 // Cluster selection
bc41680b 2340 //-----------------------------
9e51e29a 2341 Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
da13b2b8 2342 if(!ClusterSelected(calo,nMaxima)) continue;
85c4406e 2343
6175da48 2344 //----------------------------
26569ec8 2345 // Create AOD for analysis
6175da48 2346 //----------------------------
da13b2b8 2347 AliAODPWG4Particle aodph = AliAODPWG4Particle(fMomentum);
6175da48 2348
2349 //...............................................
85c4406e 2350 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
6175da48 2351 Int_t label = calo->GetLabel();
2352 aodph.SetLabel(label);
6175da48 2353 aodph.SetCaloLabel(calo->GetID(),-1);
bdb39dbd 2354 aodph.SetDetectorTag(GetCalorimeter());
c4a7d28a 2355 //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
521636d2 2356
6175da48 2357 //...............................................
2358 //Set bad channel distance bit
c4a7d28a 2359 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
f37fa8d2 2360 if (distBad > fMinDist3) aodph.SetDistToBad(2) ;
85c4406e 2361 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
f37fa8d2 2362 else aodph.SetDistToBad(0) ;
af7b3903 2363 //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
c8fe2783 2364
bc41680b 2365 //-------------------------------------
8d6b7f60 2366 // Play with the MC stack if available
bc41680b 2367 //-------------------------------------
8d6b7f60 2368
2369 //Check origin of the candidates
2370 Int_t tag = -1;
2371
34c16486 2372 if(IsDataMC())
2373 {
0cea6003 2374 tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(),GetCalorimeter());
8d6b7f60 2375 aodph.SetTag(tag);
2376
2db10729 2377 AliDebug(1,Form("Origin of candidate, bit map %d",aodph.GetTag()));
85c4406e 2378 }//Work with stack also
2379
bc41680b 2380 //--------------------------------------------------------
26569ec8 2381 // Fill some shower shape histograms before PID is applied
bc41680b 2382 //--------------------------------------------------------
521636d2 2383
26569ec8 2384 Float_t maxCellFraction = 0;
2385 Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, calo, maxCellFraction);
2386 if( absIdMax < 0 ) AliFatal("Wrong absID");
2387
2388 FillShowerShapeHistograms(calo,tag,maxCellFraction);
2389
4428fe32 2390 aodph.SetM02(calo->GetM02());
26569ec8 2391 aodph.SetNLM(nMaxima);
5cd6b1a6 2392 aodph.SetTime(calo->GetTOF()*1e9);
2393 aodph.SetNCells(calo->GetNCells());
2394 Int_t nSM = GetModuleNumber(calo);
2395 aodph.SetSModNumber(nSM);
2396
6175da48 2397 //-------------------------------------
26569ec8 2398 // PID selection or bit setting
6175da48 2399 //-------------------------------------
49b5c49b 2400
6175da48 2401 //...............................................
2402 // Data, PID check on
3c1d9afb 2403 if(IsCaloPIDOn())
2404 {
49b5c49b 2405 // Get most probable PID, 2 options check bayesian PID weights or redo PID
2406 // By default, redo PID
09273901 2407
3c1d9afb 2408 aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(calo));
477d6cee 2409
2db10729 2410 AliDebug(1,Form("PDG of identified particle %d",aodph.GetIdentifiedParticleType()));
477d6cee 2411
f37fa8d2 2412 //If cluster does not pass pid, not photon, skip it.
85c4406e 2413 if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
477d6cee 2414
2415 }
85c4406e 2416
6175da48 2417 //...............................................
2418 // Data, PID check off
3c1d9afb 2419 else
2420 {
f37fa8d2 2421 //Set PID bits for later selection (AliAnaPi0 for example)
49b5c49b 2422 //GetIdentifiedParticleType already called in SetPIDBits.
2423
3c1d9afb 2424 GetCaloPID()->SetPIDBits(calo,&aodph, GetCaloUtils(),GetReader()->GetInputEvent());
49b5c49b 2425
2db10729 2426 AliDebug(1,"PID Bits set");
477d6cee 2427 }
2428
2db10729 2429 AliDebug(1,Form("Photon selection cuts passed: pT %3.2f, pdg %d",aodph.Pt(),aodph.GetIdentifiedParticleType()));
09273901 2430
58ea8ce5 2431 fhClusterCutsE [9]->Fill(calo->E());
da13b2b8 2432 fhClusterCutsPt[9]->Fill(fMomentum.Pt());
85c4406e 2433
6df33fcb 2434 if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0)
2435 {
da13b2b8 2436 fhEPhotonSM ->Fill(fMomentum.E (),nSM);
2437 fhPtPhotonSM->Fill(fMomentum.Pt(),nSM);
6df33fcb 2438 }
2439
9e51e29a 2440 fhNLocMax->Fill(calo->E(),nMaxima);
85c4406e 2441
26569ec8 2442 // Few more control histograms for selected clusters
2443 fhMaxCellDiffClusterE->Fill(calo->E() ,maxCellFraction);
2444 fhNCellsE ->Fill(calo->E() ,calo->GetNCells());
da13b2b8 2445 fhTimePt ->Fill(fMomentum.Pt() ,calo->GetTOF()*1.e9);
26569ec8 2446
2447 if(cells)
2448 {
2449 for(Int_t icell = 0; icell < calo->GetNCells(); icell++)
2450 fhCellsE->Fill(calo->E(),cells->GetCellAmplitude(calo->GetCellsAbsId()[icell]));
2451 }
2452
09273901 2453 // Matching after cuts
bc41680b 2454 if( fFillTMHisto ) FillTrackMatchingResidualHistograms(calo,1);
09273901 2455
2ad19c3d 2456 // Fill histograms to undertand pile-up before other cuts applied
2457 // Remember to relax time cuts in the reader
26569ec8 2458 if( IsPileUpAnalysisOn() ) FillPileUpHistograms(calo,cells, absIdMax);
2ad19c3d 2459
26569ec8 2460 // Add AOD with photon object to aod branch
477d6cee 2461 AddAODParticle(aodph);
2462
2463 }//loop
85c4406e 2464
2db10729 2465 AliDebug(1,Form("End fill AODs, with %d entries",GetOutputAODBranch()->GetEntriesFast()));
477d6cee 2466
1c5acb87 2467}
2468
bc41680b 2469//______________________________________________
85c4406e 2470void AliAnaPhoton::MakeAnalysisFillHistograms()
1c5acb87 2471{
6175da48 2472 //Fill histograms
85c4406e 2473
f27fe026 2474 //In case of simulated data, fill acceptance histograms
2475 if(IsDataMC()) FillAcceptanceHistograms();
2476
6175da48 2477 // Get vertex
2244659d 2478 Double_t v[3] = {0,0,0}; //vertex ;
2479 GetReader()->GetVertex(v);
85c4406e 2480 //fhVertex->Fill(v[0],v[1],v[2]);
6175da48 2481 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
2482
2483 //----------------------------------
577d9801 2484 //Loop on stored AOD photons
2485 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2db10729 2486 AliDebug(1,Form("AOD branch entries %d", naod));
521636d2 2487
c8710850 2488 Float_t cen = GetEventCentrality();
85c4406e 2489 // printf("++++++++++ GetEventCentrality() %f\n",cen);
2490
c8710850 2491 Float_t ep = GetEventPlaneAngle();
2492
3c1d9afb 2493 for(Int_t iaod = 0; iaod < naod ; iaod++)
2494 {
577d9801 2495 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2496 Int_t pdg = ph->GetIdentifiedParticleType();
521636d2 2497
2db10729 2498 AliDebug(2,Form("PDG %d, MC TAG %d, Calorimeter <%d>",ph->GetIdentifiedParticleType(),ph->GetTag(), ph->GetDetectorTag())) ;
521636d2 2499
0cea6003 2500 //If PID used, fill histos with photons in Calorimeter GetCalorimeter()
85c4406e 2501 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
bc41680b 2502
bdb39dbd 2503 if(((Int_t) ph->GetDetectorTag()) != GetCalorimeter()) continue;
521636d2 2504
2db10729 2505 AliDebug(2,Form("ID Photon: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta())) ;
521636d2 2506
6175da48 2507 //................................
85c4406e 2508 //Fill photon histograms
577d9801 2509 Float_t ptcluster = ph->Pt();
2510 Float_t phicluster = ph->Phi();
2511 Float_t etacluster = ph->Eta();
2512 Float_t ecluster = ph->E();
521636d2 2513
20218aea 2514 fhEPhoton ->Fill(ecluster);
577d9801 2515 fhPtPhoton ->Fill(ptcluster);
2516 fhPhiPhoton ->Fill(ptcluster,phicluster);
85c4406e 2517 fhEtaPhoton ->Fill(ptcluster,etacluster);
fad96885 2518 if (ecluster > 0.5) fhEtaPhiPhoton ->Fill(etacluster, phicluster);
20218aea 2519 else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
85c4406e 2520
0cea6003 2521 if(IsHighMultiplicityAnalysisOn())
2522 {
2523 fhPtCentralityPhoton ->Fill(ptcluster,cen) ;
2524 fhPtEventPlanePhoton ->Fill(ptcluster,ep ) ;
2525 }
26569ec8 2526
2527// Comment this part, not needed but in case to know how to do it in the future
2528// //Get original cluster, to recover some information
2529// AliVCaloCells* cells = 0;
2530// TObjArray * clusters = 0;
bdb39dbd 2531// if(GetCalorimeter() == kEMCAL)
26569ec8 2532// {
2533// cells = GetEMCALCells();
2534// clusters = GetEMCALClusters();
2535// }
2536// else
2537// {
2538// cells = GetPHOSCells();
2539// clusters = GetPHOSClusters();
2540// }
2541
2542// Int_t iclus = -1;
2543// AliVCluster *cluster = FindCluster(clusters,ph->GetCaloLabel(0),iclus);
2544// if(cluster)
5812a064 2545
6175da48 2546 //.......................................
577d9801 2547 //Play with the MC data if available
34c16486 2548 if(IsDataMC())
2549 {
4c8f7c2e 2550 //....................................................................
2551 // Access MC information in stack if requested, check that it exists.
2db10729 2552 Int_t label = ph->GetLabel();
51a0ace5 2553
85c4406e 2554 if(label < 0)
34c16486 2555 {
2db10729 2556 AliDebug(1,Form("*** bad label ***: label %d", label));
4c8f7c2e 2557 continue;
2558 }
2559
2560 Float_t eprim = 0;
2561 Float_t ptprim = 0;
51a0ace5 2562 Bool_t ok = kFALSE;
1a8c88c1 2563 fPrimaryMom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
2db10729 2564
51a0ace5 2565 if(ok)
34c16486 2566 {
1a8c88c1 2567 eprim = fPrimaryMom.Energy();
2568 ptprim = fPrimaryMom.Pt();
4c8f7c2e 2569 }
2570
577d9801 2571 Int_t tag =ph->GetTag();
51a0ace5 2572 Int_t mcParticleTag = -1;
c5693f62 2573 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[kmcPhoton])
3d5d5078 2574 {
8e81c2cf 2575 fhMCE [kmcPhoton] ->Fill(ecluster);
2576 fhMCPt [kmcPhoton] ->Fill(ptcluster);
2577 fhMCPhi[kmcPhoton] ->Fill(ecluster,phicluster);
2578 fhMCEta[kmcPhoton] ->Fill(ecluster,etacluster);
2579
2580 fhMC2E [kmcPhoton] ->Fill(ecluster, eprim);
85c4406e 2581 fhMC2Pt [kmcPhoton] ->Fill(ptcluster, ptprim);
8e81c2cf 2582 fhMCDeltaE [kmcPhoton] ->Fill(ecluster,eprim-ecluster);
85c4406e 2583 fhMCDeltaPt[kmcPhoton] ->Fill(ptcluster,ptprim-ptcluster);
8e81c2cf 2584
85c4406e 2585 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) &&
764ab1f4 2586 fhMCE[kmcConversion])
3d5d5078 2587 {
8e81c2cf 2588 fhMCE [kmcConversion] ->Fill(ecluster);
2589 fhMCPt [kmcConversion] ->Fill(ptcluster);
2590 fhMCPhi[kmcConversion] ->Fill(ecluster,phicluster);
2591 fhMCEta[kmcConversion] ->Fill(ecluster,etacluster);
2592
2593 fhMC2E [kmcConversion] ->Fill(ecluster, eprim);
85c4406e 2594 fhMC2Pt [kmcConversion] ->Fill(ptcluster, ptprim);
8e81c2cf 2595 fhMCDeltaE [kmcConversion] ->Fill(ecluster,eprim-ecluster);
85c4406e 2596 fhMCDeltaPt[kmcConversion] ->Fill(ptcluster,ptprim-ptcluster);
2597 }
3d5d5078 2598
57cfa545 2599 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) )
34c16486 2600 {
51a0ace5 2601 mcParticleTag = kmcPrompt;
3d5d5078 2602 }
57cfa545 2603 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) )
3d5d5078 2604 {
51a0ace5 2605 mcParticleTag = kmcFragmentation;
3d5d5078 2606 }
57cfa545 2607 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) )
3d5d5078 2608 {
85c4406e 2609 mcParticleTag = kmcISR;
3d5d5078 2610 }
57cfa545 2611 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) )
3d5d5078 2612 {
85c4406e 2613 mcParticleTag = kmcPi0;
2614 }
57cfa545 2615 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
f66d95af 2616 {
51a0ace5 2617 mcParticleTag = kmcEta;
85c4406e 2618 }
57cfa545 2619 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) )
2620 {
2621 mcParticleTag = kmcPi0Decay;
2622 }
2623 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) )
2624 {
2625 mcParticleTag = kmcEtaDecay;
2626 }
760b98f5 2627 else
2628 {
2629 mcParticleTag = kmcOtherDecay;
2630 }
3d5d5078 2631 }
57cfa545 2632 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) )
3d5d5078 2633 {
51a0ace5 2634 mcParticleTag = kmcAntiNeutron;
3d5d5078 2635 }
57cfa545 2636 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) )
3d5d5078 2637 {
85c4406e 2638 mcParticleTag = kmcAntiProton;
2639 }
57cfa545 2640 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) )
3d5d5078 2641 {
85c4406e 2642 mcParticleTag = kmcElectron;
2643 }
57cfa545 2644 else if( fhMCE[kmcOther] )
34c16486 2645 {
51a0ace5 2646 mcParticleTag = kmcOther;
4c8f7c2e 2647
f8006433 2648 // printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
2649 // ph->GetLabel(),ph->Pt());
2650 // for(Int_t i = 0; i < 20; i++) {
2651 // if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
2652 // }
2653 // printf("\n");
2654
577d9801 2655 }
521636d2 2656
760b98f5 2657 if(mcParticleTag >= 0 && fhMCE[mcParticleTag])
fde324ab 2658 {
2659 fhMCE [mcParticleTag]->Fill(ecluster);
2660 fhMCPt [mcParticleTag]->Fill(ptcluster);
2661 fhMCPhi[mcParticleTag]->Fill(ecluster,phicluster);
2662 fhMCEta[mcParticleTag]->Fill(ecluster,etacluster);
2663
2664 fhMC2E [mcParticleTag]->Fill(ecluster, eprim);
2665 fhMC2Pt [mcParticleTag]->Fill(ptcluster, ptprim);
2666 fhMCDeltaE [mcParticleTag]->Fill(ecluster,eprim-ecluster);
2667 fhMCDeltaPt[mcParticleTag]->Fill(ptcluster,ptprim-ptcluster);
2668 }
577d9801 2669 }//Histograms with MC
521636d2 2670
577d9801 2671 }// aod loop
521636d2 2672
1c5acb87 2673}
2674
2675
2676//__________________________________________________________________
2677void AliAnaPhoton::Print(const Option_t * opt) const
2678{
477d6cee 2679 //Print some relevant parameters set for the analysis
2680
2681 if(! opt)
2682 return;
2683
2684 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 2685 AliAnaCaloTrackCorrBaseClass::Print(" ");
85c4406e 2686
bdb39dbd 2687 printf("Calorimeter = %s\n", GetCalorimeterString().Data()) ;
477d6cee 2688 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
2689 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2690 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 2691 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
4cf55759 2692 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2ac125bf 2693 printf("Number of cells in cluster is > %d \n", fNCellsCut);
477d6cee 2694 printf(" \n") ;
1c5acb87 2695
85c4406e 2696}