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