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