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