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