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