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