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