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