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