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