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