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