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