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