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