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