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