]> 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
Be sure to load mapping when needed
[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 **************************************************************************/
15/* $Id: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
16
17//_________________________________________________________________________
18//
19// Class for the photon identification.
20// Clusters from calorimeters are identified as photons
21// and kept in the AOD. Few histograms produced.
6175da48 22// Produces input for other analysis classes like AliAnaPi0,
23// AliAnaParticleHadronCorrelation ...
1c5acb87 24//
25// -- Author: Gustavo Conesa (LNF-INFN)
26//////////////////////////////////////////////////////////////////////////////
27
28
29// --- ROOT system ---
30#include <TH2F.h>
2244659d 31#include <TH3D.h>
477d6cee 32#include <TClonesArray.h>
0c1383b5 33#include <TObjString.h>
477d6cee 34//#include <Riostream.h>
123fc3bd 35#include "TParticle.h"
6175da48 36#include "TDatabasePDG.h"
1c5acb87 37
38// --- Analysis system ---
39#include "AliAnaPhoton.h"
40#include "AliCaloTrackReader.h"
123fc3bd 41#include "AliStack.h"
1c5acb87 42#include "AliCaloPID.h"
6639984f 43#include "AliMCAnalysisUtils.h"
ff45398a 44#include "AliFiducialCut.h"
0ae57829 45#include "AliVCluster.h"
591cc579 46#include "AliAODMCParticle.h"
c8fe2783 47#include "AliMixedEvent.h"
48
1c5acb87 49
50ClassImp(AliAnaPhoton)
51
52//____________________________________________________________________________
521636d2 53AliAnaPhoton::AliAnaPhoton() :
54 AliAnaPartCorrBaseClass(), fCalorimeter(""),
55 fMinDist(0.), fMinDist2(0.), fMinDist3(0.),
56 fRejectTrackMatch(0), fTimeCutMin(-1), fTimeCutMax(999999),
f66d95af 57 fNCellsCut(0), fFillSSHistograms(kFALSE),
58 fNOriginHistograms(8), fNPrimaryHistograms(4),
521636d2 59 fCheckConversion(kFALSE), fRemoveConvertedPair(kFALSE),
60 fAddConvertedPairsToAOD(kFALSE),
61 fMassCut(0), fConvAsymCut(1.), fConvDEtaCut(2.),
62 fConvDPhiMinCut(-1.), fConvDPhiMaxCut(7.),
63
c4a7d28a 64 // Histograms
f66d95af 65 fhNCellsE(0), fhMaxCellDiffClusterE(0),
66 fhEPhoton(0), fhPtPhoton(0),
521636d2 67 fhPhiPhoton(0), fhEtaPhoton(0),
68 fhEtaPhiPhoton(0), fhEtaPhi05Photon(0),
69
c4a7d28a 70 // Conversion histograms
521636d2 71 fhPtPhotonConv(0), fhEtaPhiPhotonConv(0), fhEtaPhi05PhotonConv(0),
72 fhConvDeltaEta(0), fhConvDeltaPhi(0), fhConvDeltaEtaPhi(0),
73 fhConvAsym(0), fhConvPt(0),
74 fhConvDistEta(0), fhConvDistEn(0), fhConvDistMass(0),
75 fhConvDistEtaCutEta(0), fhConvDistEnCutEta(0), fhConvDistMassCutEta(0),
76 fhConvDistEtaCutMass(0), fhConvDistEnCutMass(0),
77 fhConvDistEtaCutAsy(0), fhConvDistEnCutAsy(0),
78
79 //Shower shape histograms
80 fhDispE(0), fhLam0E(0), fhLam1E(0),
81 fhdDispE(0), fhdLam0E(0), fhdLam1E(0),
82 fhDispETRD(0), fhLam0ETRD(0), fhLam1ETRD(0),
83 fhdDispETRD(0), fhdLam0ETRD(0), fhdLam1ETRD(0),
84
85 fhNCellsLam0LowE(0), fhNCellsLam1LowE(0), fhNCellsDispLowE(0),
86 fhNCellsLam0HighE(0), fhNCellsLam1HighE(0), fhNCellsDispHighE(0),
87 fhNCellsdLam0LowE(0), fhNCellsdLam1LowE(0), fhNCellsdDispLowE(0),
88 fhNCellsdLam0HighE(0), fhNCellsdLam1HighE(0), fhNCellsdDispHighE(0),
89
90 fhEtaLam0LowE(0), fhPhiLam0LowE(0),
91 fhEtaLam0HighE(0), fhPhiLam0HighE(0),
92 fhLam0DispLowE(0), fhLam0DispHighE(0),
93 fhLam1Lam0LowE(0), fhLam1Lam0HighE(0),
94 fhDispLam1LowE(0), fhDispLam1HighE(0),
95 fhEtadLam0LowE(0), fhPhidLam0LowE(0),
96 fhEtadLam0HighE(0), fhPhidLam0HighE(0),
97 fhdLam0dDispLowE(0), fhdLam0dDispHighE(0),
98 fhdLam1dLam0LowE(0), fhdLam1dLam0HighE(0),
99 fhdDispdLam1LowE(0), fhdDispdLam1HighE(0),
100
101 //MC histograms
102 fhDeltaE(0), fhDeltaPt(0),
103 fhRatioE(0), fhRatioPt(0),
104 fh2E(0), fh2Pt(0),
105
106 // Conversion MC histograms
107 fhPtConversionTagged(0), fhPtAntiNeutronTagged(0),
108 fhPtAntiProtonTagged(0), fhPtUnknownTagged(0),
109 fhEtaPhiConversion(0), fhEtaPhi05Conversion(0),
110
111 fhConvDeltaEtaMCConversion(0), fhConvDeltaPhiMCConversion(0), fhConvDeltaEtaPhiMCConversion(0),
112 fhConvAsymMCConversion(0), fhConvPtMCConversion(0),
113 fhConvDispersionMCConversion(0), fhConvM02MCConversion(0),
114
115 fhConvDeltaEtaMCAntiNeutron(0), fhConvDeltaPhiMCAntiNeutron(0), fhConvDeltaEtaPhiMCAntiNeutron(0),
116 fhConvAsymMCAntiNeutron(0), fhConvPtMCAntiNeutron(0),
117 fhConvDispersionMCAntiNeutron(0), fhConvM02MCAntiNeutron(0),
118 fhConvDeltaEtaMCAntiProton(0), fhConvDeltaPhiMCAntiProton(0), fhConvDeltaEtaPhiMCAntiProton(0),
119 fhConvAsymMCAntiProton(0), fhConvPtMCAntiProton(0),
120 fhConvDispersionMCAntiProton(0), fhConvM02MCAntiProton(0),
121 fhConvDeltaEtaMCString(0), fhConvDeltaPhiMCString(0), fhConvDeltaEtaPhiMCString(0),
122 fhConvAsymMCString(0), fhConvPtMCString(0),
123 fhConvDispersionMCString(0), fhConvM02MCString(0),
3d5d5078 124 fhConvDistMCConversion(0), fhConvDistMCConversionCuts(0),
125
126 // Photon SS MC histograms
127 fhMCPhotonELambda0NoOverlap(0), fhMCPhotonELambda0TwoOverlap(0), fhMCPhotonELambda0NOverlap(0),
128 fhMCPhotonEdLambda0NoOverlap(0), fhMCPhotonEdLambda0TwoOverlap(0), fhMCPhotonEdLambda0NOverlap(0),
129
130 //Embedding
131 fhEmbeddedSignalFractionEnergy(0),
132 fhEmbedPhotonELambda0FullSignal(0), fhEmbedPhotonEdLambda0FullSignal(0),
133 fhEmbedPhotonELambda0MostlySignal(0), fhEmbedPhotonEdLambda0MostlySignal(0),
134 fhEmbedPhotonELambda0MostlyBkg(0), fhEmbedPhotonEdLambda0MostlyBkg(0),
135 fhEmbedPhotonELambda0FullBkg(0), fhEmbedPhotonEdLambda0FullBkg(0),
136 fhEmbedPi0ELambda0FullSignal(0), fhEmbedPi0EdLambda0FullSignal(0),
137 fhEmbedPi0ELambda0MostlySignal(0), fhEmbedPi0EdLambda0MostlySignal(0),
138 fhEmbedPi0ELambda0MostlyBkg(0), fhEmbedPi0EdLambda0MostlyBkg(0),
139 fhEmbedPi0ELambda0FullBkg(0), fhEmbedPi0EdLambda0FullBkg(0)
1c5acb87 140{
141 //default ctor
142
f66d95af 143 for(Int_t i = 0; i < 14; i++){
3d5d5078 144 fhPtMC [i] = 0;
145 fhMCE [i] = 0;
521636d2 146 fhPhiMC[i] = 0;
147 fhEtaMC[i] = 0;
148 }
149
3d5d5078 150 for(Int_t i = 0; i < 7; i++){
151 fhPtPrimMC [i] = 0;
152 fhEPrimMC [i] = 0;
153 fhPhiPrimMC[i] = 0;
154 fhYPrimMC [i] = 0;
155
156 fhPtPrimMCAcc [i] = 0;
157 fhEPrimMCAcc [i] = 0;
158 fhPhiPrimMCAcc[i] = 0;
159 fhYPrimMCAcc [i] = 0;
160 }
161
162 for(Int_t i = 0; i < 6; i++){
f66d95af 163 fhMCELambda0 [i] = 0;
164 fhMCELambda1 [i] = 0;
165 fhMCEDispersion [i] = 0;
166 fhMCEdLambda0 [i] = 0;
167 fhMCEdLambda1 [i] = 0;
168 fhMCEdDispersion[i] = 0;
169 fhMCNCellsE [i] = 0;
170 fhMCMaxCellDiffClusterE[i] = 0;
171 fhMCLambda0vsClusterMaxCellDiffE0[i] = 0;
172 fhMCLambda0vsClusterMaxCellDiffE2[i] = 0;
173 fhMCLambda0vsClusterMaxCellDiffE6[i] = 0;
174 fhMCNCellsvsClusterMaxCellDiffE0 [i] = 0;
175 fhMCNCellsvsClusterMaxCellDiffE2 [i] = 0;
176 fhMCNCellsvsClusterMaxCellDiffE6 [i] = 0;
521636d2 177 }
178
1c5acb87 179 //Initialize parameters
180 InitParameters();
181
5ae09196 182}//____________________________________________________________________________
1c5acb87 183AliAnaPhoton::~AliAnaPhoton()
184{
185 //dtor
186
187}
188
c4a7d28a 189//__________________________________________________________________
190Bool_t AliAnaPhoton::ClusterSelected(AliVCluster* calo, TLorentzVector mom)
191{
192 //Select clusters if they pass different cuts
193 if(GetDebug() > 2)
194 printf("AliAnaPhoton::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
195 GetReader()->GetEventNumber(),
196 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
197
198 //.......................................
199 //If too small or big energy, skip it
200 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ;
201 if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
202
203 //.......................................
204 // TOF cut, BE CAREFUL WITH THIS CUT
205 Double_t tof = calo->GetTOF()*1e9;
206 if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
207 if(GetDebug() > 2) printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
208
209 //.......................................
210 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
211 if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
212
213 //.......................................
214 //Check acceptance selection
215 if(IsFiducialCutOn()){
216 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
217 if(! in ) return kFALSE ;
218 }
219 if(GetDebug() > 2) printf("Fiducial cut passed \n");
220
221 //.......................................
222 //Skip matched clusters with tracks
223 if(fRejectTrackMatch){
224 if(IsTrackMatched(calo)) {
225 if(GetDebug() > 2) printf("\t Reject track-matched clusters\n");
226 return kFALSE ;
227 }
228 else
229 if(GetDebug() > 2) printf(" Track-matching cut passed \n");
230 }// reject matched clusters
231
232 //.......................................
233 //Check Distance to Bad channel, set bit.
234 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
235 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
236 if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
237 return kFALSE ;
238 }
239 else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
240 //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
241
242 if(GetDebug() > 0)
243 printf("AliAnaPhoton::ClusterSelected() Current Event %d; After selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
244 GetReader()->GetEventNumber(),
245 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
246
247 //All checks passed, cluster selected
248 return kTRUE;
249
250}
251
3d5d5078 252//_____________________________________________________________
253void AliAnaPhoton::FillAcceptanceHistograms(){
254 //Fill acceptance histograms if MC data is available
255
256 if(GetReader()->ReadStack()){
257 AliStack * stack = GetMCStack();
258 if(stack){
259 for(Int_t i=0 ; i<stack->GetNtrack(); i++){
260 TParticle * prim = stack->Particle(i) ;
261 Int_t pdg = prim->GetPdgCode();
262 //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
263 // prim->GetName(), prim->GetPdgCode());
264
265 if(pdg == 22){
266
267 // Get tag of this particle photon from fragmentation, decay, prompt ...
268 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
269 if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)){
270 //A conversion photon from a hadron, skip this kind of photon
271 // 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,
272 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
273 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
274 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
275 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
276 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
277 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon),
278 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
279 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton),
280 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
281 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon),
282 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton),
283 // GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron));
284
285 return;
286 }
287
288 //Get photon kinematics
289 if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
290
291 Double_t photonY = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
292 Double_t photonE = prim->Energy() ;
293 Double_t photonPt = prim->Pt() ;
294 Double_t photonPhi = TMath::RadToDeg()*prim->Phi() ;
295 if(photonPhi < 0) photonPhi+=TMath::TwoPi();
296 Double_t photonEta = prim->Eta() ;
297
298
299 //Check if photons hit the Calorimeter
300 TLorentzVector lv;
301 prim->Momentum(lv);
302 Bool_t inacceptance = kFALSE;
303 if (fCalorimeter == "PHOS"){
304 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
305 Int_t mod ;
306 Double_t x,z ;
307 if(GetPHOSGeometry()->ImpactOnEmc(prim,mod,z,x))
308 inacceptance = kTRUE;
309 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
310 }
311 else{
312 if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
313 inacceptance = kTRUE ;
314 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
315 }
316 }
317 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
318 if(GetEMCALGeometry()){
319
320 Int_t absID=0;
321
322 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
323
324 if( absID >= 0)
325 inacceptance = kTRUE;
326
327 // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
328 // inacceptance = kTRUE;
329 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
330 }
331 else{
332 if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
333 inacceptance = kTRUE ;
334 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
335 }
336 } //In EMCAL
337
338 //Fill histograms
339
f66d95af 340 fhYPrimMC[mcPPhoton]->Fill(photonPt, photonY) ;
3d5d5078 341 if(TMath::Abs(photonY) < 1.0)
342 {
f66d95af 343 fhEPrimMC [mcPPhoton]->Fill(photonE ) ;
344 fhPtPrimMC [mcPPhoton]->Fill(photonPt) ;
345 fhPhiPrimMC[mcPPhoton]->Fill(photonE , photonPhi) ;
346 fhYPrimMC[mcPPhoton] ->Fill(photonE , photonEta) ;
3d5d5078 347 }
348 if(inacceptance){
f66d95af 349 fhEPrimMCAcc[mcPPhoton] ->Fill(photonE ) ;
350 fhPtPrimMCAcc[mcPPhoton] ->Fill(photonPt) ;
351 fhPhiPrimMCAcc[mcPPhoton]->Fill(photonE , photonPhi) ;
352 fhYPrimMCAcc[mcPPhoton] ->Fill(photonE , photonY) ;
3d5d5078 353 }//Accepted
354
355 //Origin of photon
f66d95af 356 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[mcPPrompt])
3d5d5078 357 {
f66d95af 358 fhYPrimMC[mcPPrompt]->Fill(photonPt, photonY) ;
3d5d5078 359 if(TMath::Abs(photonY) < 1.0){
f66d95af 360 fhEPrimMC [mcPPrompt]->Fill(photonE ) ;
361 fhPtPrimMC [mcPPrompt]->Fill(photonPt) ;
362 fhPhiPrimMC[mcPPrompt]->Fill(photonE , photonPhi) ;
363 fhYPrimMC[mcPPrompt] ->Fill(photonE , photonEta) ;
3d5d5078 364 }
365 if(inacceptance){
f66d95af 366 fhEPrimMCAcc[mcPPrompt] ->Fill(photonE ) ;
367 fhPtPrimMCAcc[mcPPrompt] ->Fill(photonPt) ;
368 fhPhiPrimMCAcc[mcPPrompt]->Fill(photonE , photonPhi) ;
369 fhYPrimMCAcc[mcPPrompt] ->Fill(photonE , photonY) ;
3d5d5078 370 }//Accepted
371 }
f66d95af 372 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[mcPFragmentation])
3d5d5078 373 {
f66d95af 374 fhYPrimMC[mcPFragmentation]->Fill(photonPt, photonY) ;
3d5d5078 375 if(TMath::Abs(photonY) < 1.0){
f66d95af 376 fhEPrimMC [mcPFragmentation]->Fill(photonE ) ;
377 fhPtPrimMC [mcPFragmentation]->Fill(photonPt) ;
378 fhPhiPrimMC[mcPFragmentation]->Fill(photonE , photonPhi) ;
379 fhYPrimMC[mcPFragmentation] ->Fill(photonE , photonEta) ;
3d5d5078 380 }
381 if(inacceptance){
f66d95af 382 fhEPrimMCAcc[mcPFragmentation] ->Fill(photonE ) ;
383 fhPtPrimMCAcc[mcPFragmentation] ->Fill(photonPt) ;
384 fhPhiPrimMCAcc[mcPFragmentation]->Fill(photonE , photonPhi) ;
385 fhYPrimMCAcc[mcPFragmentation] ->Fill(photonE , photonY) ;
3d5d5078 386 }//Accepted
387 }
f66d95af 388 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[mcPISR])
3d5d5078 389 {
f66d95af 390 fhYPrimMC[mcPISR]->Fill(photonPt, photonY) ;
3d5d5078 391 if(TMath::Abs(photonY) < 1.0){
f66d95af 392 fhEPrimMC [mcPISR]->Fill(photonE ) ;
393 fhPtPrimMC [mcPISR]->Fill(photonPt) ;
394 fhPhiPrimMC[mcPISR]->Fill(photonE , photonPhi) ;
395 fhYPrimMC[mcPISR]->Fill(photonE , photonEta) ;
3d5d5078 396 }
397 if(inacceptance){
f66d95af 398 fhEPrimMCAcc[mcPISR] ->Fill(photonE ) ;
399 fhPtPrimMCAcc[mcPISR] ->Fill(photonPt) ;
400 fhPhiPrimMCAcc[mcPISR]->Fill(photonE , photonPhi) ;
401 fhYPrimMCAcc[mcPISR] ->Fill(photonE , photonY) ;
3d5d5078 402 }//Accepted
403 }
f66d95af 404 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[mcPPi0Decay])
3d5d5078 405 {
f66d95af 406 fhYPrimMC[mcPPi0Decay]->Fill(photonPt, photonY) ;
3d5d5078 407 if(TMath::Abs(photonY) < 1.0){
f66d95af 408 fhEPrimMC [mcPPi0Decay]->Fill(photonE ) ;
409 fhPtPrimMC [mcPPi0Decay]->Fill(photonPt) ;
410 fhPhiPrimMC[mcPPi0Decay]->Fill(photonE , photonPhi) ;
411 fhYPrimMC[mcPPi0Decay] ->Fill(photonE , photonEta) ;
3d5d5078 412 }
413 if(inacceptance){
f66d95af 414 fhEPrimMCAcc[mcPPi0Decay] ->Fill(photonE ) ;
415 fhPtPrimMCAcc[mcPPi0Decay] ->Fill(photonPt) ;
416 fhPhiPrimMCAcc[mcPPi0Decay]->Fill(photonE , photonPhi) ;
417 fhYPrimMCAcc[mcPPi0Decay] ->Fill(photonE , photonY) ;
3d5d5078 418 }//Accepted
419 }
420 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
f66d95af 421 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) && fhEPrimMC[mcPOtherDecay])
3d5d5078 422 {
f66d95af 423 fhYPrimMC[mcPOtherDecay]->Fill(photonPt, photonY) ;
3d5d5078 424 if(TMath::Abs(photonY) < 1.0){
f66d95af 425 fhEPrimMC [mcPOtherDecay]->Fill(photonE ) ;
426 fhPtPrimMC [mcPOtherDecay]->Fill(photonPt) ;
427 fhPhiPrimMC[mcPOtherDecay]->Fill(photonE , photonPhi) ;
428 fhYPrimMC[mcPOtherDecay] ->Fill(photonE , photonEta) ;
3d5d5078 429 }
430 if(inacceptance){
f66d95af 431 fhEPrimMCAcc[mcPOtherDecay] ->Fill(photonE ) ;
432 fhPtPrimMCAcc[mcPOtherDecay] ->Fill(photonPt) ;
433 fhPhiPrimMCAcc[mcPOtherDecay]->Fill(photonE , photonPhi) ;
434 fhYPrimMCAcc[mcPOtherDecay] ->Fill(photonE , photonY) ;
3d5d5078 435 }//Accepted
436 }
f66d95af 437 else if(fhEPrimMC[mcPOther])
3d5d5078 438 {
f66d95af 439 fhYPrimMC[mcPOther]->Fill(photonPt, photonY) ;
3d5d5078 440 if(TMath::Abs(photonY) < 1.0){
f66d95af 441 fhEPrimMC [mcPOther]->Fill(photonE ) ;
442 fhPtPrimMC [mcPOther]->Fill(photonPt) ;
443 fhPhiPrimMC[mcPOther]->Fill(photonE , photonPhi) ;
444 fhYPrimMC[mcPOther] ->Fill(photonE , photonEta) ;
3d5d5078 445 }
446 if(inacceptance){
f66d95af 447 fhEPrimMCAcc[mcPOther] ->Fill(photonE ) ;
448 fhPtPrimMCAcc[mcPOther] ->Fill(photonPt) ;
449 fhPhiPrimMCAcc[mcPOther]->Fill(photonE , photonPhi) ;
450 fhYPrimMCAcc[mcPOther] ->Fill(photonE , photonY) ;
3d5d5078 451 }//Accepted
452 }//Other origin
453 }// Primary photon
454 }//loop on primaries
455 }//stack exists and data is MC
456 }//read stack
457 else if(GetReader()->ReadAODMCParticles()){
458 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
459 if(mcparticles){
460 Int_t nprim = mcparticles->GetEntriesFast();
461
462 for(Int_t i=0; i < nprim; i++)
463 {
464 AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
465
466 Int_t pdg = prim->GetPdgCode();
467
468 if(pdg == 22){
469
470 // Get tag of this particle photon from fragmentation, decay, prompt ...
471 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
472 if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)){
473 //A conversion photon from a hadron, skip this kind of photon
474// 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,
475// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
476// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
477// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
478// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
479// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
480// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon),
481// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
482// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton),
483// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
484// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon),
485// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton),
486// GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron));
487
488 return;
489 }
490
491 //Get photon kinematics
492 if(prim->E() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
493
494 Double_t photonY = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
495 Double_t photonE = prim->E() ;
496 Double_t photonPt = prim->Pt() ;
497 Double_t photonPhi = TMath::RadToDeg()*prim->Phi() ;
498 if(photonPhi < 0) photonPhi+=TMath::TwoPi();
499 Double_t photonEta = prim->Eta() ;
500
501 //Check if photons hit the Calorimeter
502 TLorentzVector lv;
503 lv.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
504 Bool_t inacceptance = kFALSE;
505 if (fCalorimeter == "PHOS"){
506 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
507 Int_t mod ;
508 Double_t x,z ;
509 Double_t vtx[]={prim->Xv(),prim->Yv(),prim->Zv()};
510 if(GetPHOSGeometry()->ImpactOnEmc(vtx, prim->Theta(),prim->Phi(),mod,z,x))
511 inacceptance = kTRUE;
512 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
513 }
514 else{
515 if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
516 inacceptance = kTRUE ;
517 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
518 }
519 }
520 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
521 if(GetEMCALGeometry()){
522
523 Int_t absID=0;
524
525 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
526
527 if( absID >= 0)
528 inacceptance = kTRUE;
529
530 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
531 }
532 else{
533 if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
534 inacceptance = kTRUE ;
535 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
536 }
537 } //In EMCAL
538
539 //Fill histograms
540
f66d95af 541 fhYPrimMC[mcPPhoton]->Fill(photonPt, photonY) ;
3d5d5078 542 if(TMath::Abs(photonY) < 1.0)
543 {
f66d95af 544 fhEPrimMC [mcPPhoton]->Fill(photonE ) ;
545 fhPtPrimMC [mcPPhoton]->Fill(photonPt) ;
546 fhPhiPrimMC[mcPPhoton]->Fill(photonE , photonPhi) ;
547 fhYPrimMC[mcPPhoton]->Fill(photonE , photonEta) ;
3d5d5078 548 }
549 if(inacceptance){
f66d95af 550 fhEPrimMCAcc[mcPPhoton] ->Fill(photonE ) ;
551 fhPtPrimMCAcc[mcPPhoton] ->Fill(photonPt) ;
552 fhPhiPrimMCAcc[mcPPhoton]->Fill(photonE , photonPhi) ;
553 fhYPrimMCAcc[mcPPhoton] ->Fill(photonE , photonY) ;
3d5d5078 554 }//Accepted
555
556
557 //Origin of photon
f66d95af 558 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhEPrimMC[mcPPrompt])
3d5d5078 559 {
f66d95af 560 fhYPrimMC[mcPPrompt]->Fill(photonPt, photonY) ;
3d5d5078 561 if(TMath::Abs(photonY) < 1.0){
f66d95af 562 fhEPrimMC [mcPPrompt]->Fill(photonE ) ;
563 fhPtPrimMC [mcPPrompt]->Fill(photonPt) ;
564 fhPhiPrimMC[mcPPrompt]->Fill(photonE , photonPhi) ;
565 fhYPrimMC[mcPPrompt]->Fill(photonE , photonEta) ;
3d5d5078 566 }
567 if(inacceptance){
f66d95af 568 fhEPrimMCAcc[mcPPrompt] ->Fill(photonE ) ;
569 fhPtPrimMCAcc[mcPPrompt] ->Fill(photonPt) ;
570 fhPhiPrimMCAcc[mcPPrompt]->Fill(photonE , photonPhi) ;
571 fhYPrimMCAcc[mcPPrompt] ->Fill(photonE , photonY) ;
3d5d5078 572 }//Accepted
573 }
f66d95af 574 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation) && fhEPrimMC[mcPFragmentation] )
3d5d5078 575 {
f66d95af 576 fhYPrimMC[mcPFragmentation]->Fill(photonPt, photonY) ;
3d5d5078 577 if(TMath::Abs(photonY) < 1.0){
f66d95af 578 fhEPrimMC [mcPFragmentation]->Fill(photonE ) ;
579 fhPtPrimMC [mcPFragmentation]->Fill(photonPt) ;
580 fhPhiPrimMC[mcPFragmentation]->Fill(photonE , photonPhi) ;
581 fhYPrimMC[mcPFragmentation]->Fill(photonE , photonEta) ;
3d5d5078 582 }
583 if(inacceptance){
f66d95af 584 fhEPrimMCAcc[mcPFragmentation] ->Fill(photonE ) ;
585 fhPtPrimMCAcc[mcPFragmentation] ->Fill(photonPt) ;
586 fhPhiPrimMCAcc[mcPFragmentation]->Fill(photonE , photonPhi) ;
587 fhYPrimMCAcc[mcPFragmentation] ->Fill(photonE , photonY) ;
3d5d5078 588 }//Accepted
589 }
f66d95af 590 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR) && fhEPrimMC[mcPISR])
3d5d5078 591 {
f66d95af 592 fhYPrimMC[mcPISR]->Fill(photonPt, photonY) ;
3d5d5078 593 if(TMath::Abs(photonY) < 1.0){
f66d95af 594 fhEPrimMC [mcPISR]->Fill(photonE ) ;
595 fhPtPrimMC [mcPISR]->Fill(photonPt) ;
596 fhPhiPrimMC[mcPISR]->Fill(photonE , photonPhi) ;
597 fhYPrimMC[mcPISR]->Fill(photonE , photonEta) ;
3d5d5078 598 }
599 if(inacceptance){
f66d95af 600 fhEPrimMCAcc[mcPISR] ->Fill(photonE ) ;
601 fhPtPrimMCAcc[mcPISR] ->Fill(photonPt) ;
602 fhPhiPrimMCAcc[mcPISR]->Fill(photonE , photonPhi) ;
603 fhYPrimMCAcc[mcPISR] ->Fill(photonE , photonY) ;
3d5d5078 604 }//Accepted
605 }
f66d95af 606 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay)&& fhEPrimMC[mcPPi0Decay])
3d5d5078 607 {
f66d95af 608 fhYPrimMC[mcPPi0Decay]->Fill(photonPt, photonY) ;
3d5d5078 609 if(TMath::Abs(photonY) < 1.0){
f66d95af 610 fhEPrimMC [mcPPi0Decay]->Fill(photonE ) ;
611 fhPtPrimMC [mcPPi0Decay]->Fill(photonPt) ;
612 fhPhiPrimMC[mcPPi0Decay]->Fill(photonE , photonPhi) ;
613 fhYPrimMC[mcPPi0Decay]->Fill(photonE , photonEta) ;
3d5d5078 614 }
615 if(inacceptance){
f66d95af 616 fhEPrimMCAcc[mcPPi0Decay] ->Fill(photonE ) ;
617 fhPtPrimMCAcc[mcPPi0Decay] ->Fill(photonPt) ;
618 fhPhiPrimMCAcc[mcPPi0Decay]->Fill(photonE , photonPhi) ;
619 fhYPrimMCAcc[mcPPi0Decay] ->Fill(photonE , photonY) ;
3d5d5078 620 }//Accepted
621 }
622 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
f66d95af 623 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) && fhEPrimMC[mcPOtherDecay])
3d5d5078 624 {
f66d95af 625 fhYPrimMC[mcPOtherDecay]->Fill(photonPt, photonY) ;
3d5d5078 626 if(TMath::Abs(photonY) < 1.0){
f66d95af 627 fhEPrimMC [mcPOtherDecay]->Fill(photonE ) ;
628 fhPtPrimMC [mcPOtherDecay]->Fill(photonPt) ;
629 fhPhiPrimMC[mcPOtherDecay]->Fill(photonE , photonPhi) ;
630 fhYPrimMC[mcPOtherDecay]->Fill(photonE , photonEta) ;
3d5d5078 631 }
632 if(inacceptance){
f66d95af 633 fhEPrimMCAcc[mcPOtherDecay] ->Fill(photonE ) ;
634 fhPtPrimMCAcc[mcPOtherDecay] ->Fill(photonPt) ;
635 fhPhiPrimMCAcc[mcPOtherDecay]->Fill(photonE , photonPhi) ;
636 fhYPrimMCAcc[mcPOtherDecay] ->Fill(photonE , photonY) ;
3d5d5078 637 }//Accepted
638 }
f66d95af 639 else if(fhEPrimMC[mcPOther])
3d5d5078 640 {
f66d95af 641 fhYPrimMC[mcPOther]->Fill(photonPt, photonY) ;
3d5d5078 642 if(TMath::Abs(photonY) < 1.0){
f66d95af 643 fhEPrimMC [mcPOther]->Fill(photonE ) ;
644 fhPtPrimMC [mcPOther]->Fill(photonPt) ;
645 fhPhiPrimMC[mcPOther]->Fill(photonE , photonPhi) ;
646 fhYPrimMC[mcPOther]->Fill(photonE , photonEta) ;
3d5d5078 647 }
648 if(inacceptance){
f66d95af 649 fhEPrimMCAcc[mcPOther] ->Fill(photonE ) ;
650 fhPtPrimMCAcc[mcPOther] ->Fill(photonPt) ;
651 fhPhiPrimMCAcc[mcPOther]->Fill(photonE , photonPhi) ;
652 fhYPrimMCAcc[mcPOther] ->Fill(photonE , photonY) ;
3d5d5078 653 }//Accepted
654 }//Other origin
655 }// Primary photon
656 }//loop on primaries
657
658 }//mc array exists and data is MC
659 } // read AOD MC
660}
521636d2 661
662//__________________________________________________________________
663void AliAnaPhoton::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag){
664
665 //Fill cluster Shower Shape histograms
666
667 if(!fFillSSHistograms || GetMixedEvent()) return;
668
669 Float_t energy = cluster->E();
670 Int_t ncells = cluster->GetNCells();
671 Int_t ncells2 = ncells*ncells;
672 Float_t lambda0 = cluster->GetM02();
673 Float_t lambda1 = cluster->GetM20();
674 Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
675
676 TLorentzVector mom;
677 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
678 cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
679 else{
680 Double_t vertex[]={0,0,0};
681 cluster->GetMomentum(mom,vertex) ;
682 }
683
684 Float_t eta = mom.Eta();
685 Float_t phi = mom.Phi();
686 if(phi < 0) phi+=TMath::TwoPi();
687
688 fhLam0E ->Fill(energy,lambda0);
689 fhLam1E ->Fill(energy,lambda1);
690 fhDispE ->Fill(energy,disp);
691 fhdLam0E->Fill(energy,lambda0/ncells2);
692 fhdLam1E->Fill(energy,lambda1/ncells2);
693 fhdDispE->Fill(energy,disp/ncells2);
694
695 if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5){
696 fhLam0ETRD->Fill(energy,lambda0);
697 fhLam1ETRD->Fill(energy,lambda1);
698 fhDispETRD->Fill(energy,disp);
699 fhdLam0ETRD->Fill(energy,lambda0/ncells2);
700 fhdLam1ETRD->Fill(energy,lambda1/ncells2);
701 fhdDispETRD->Fill(energy,disp/ncells2);
702 }
703
704 if(energy < 2){
705 fhNCellsLam0LowE ->Fill(ncells,lambda0);
706 fhNCellsLam1LowE ->Fill(ncells,lambda1);
707 fhNCellsDispLowE ->Fill(ncells,disp);
708 fhNCellsdLam0LowE->Fill(ncells,lambda0/ncells2);
709 fhNCellsdLam1LowE->Fill(ncells,lambda1/ncells2);
710 fhNCellsdDispLowE->Fill(ncells,disp/ncells2);
711
712 fhLam1Lam0LowE ->Fill(lambda1,lambda0);
713 fhLam0DispLowE ->Fill(lambda0,disp);
714 fhDispLam1LowE ->Fill(disp,lambda1);
715 fhEtaLam0LowE ->Fill(eta,lambda0);
716 fhPhiLam0LowE ->Fill(phi,lambda0);
717
718 fhdLam1dLam0LowE->Fill(lambda1/ncells2,lambda0/ncells2);
719 fhdLam0dDispLowE->Fill(lambda0/ncells2,disp/ncells2);
720 fhdDispdLam1LowE->Fill(disp/ncells2,lambda1/ncells2);
721 fhEtadLam0LowE ->Fill(eta,lambda0/ncells2);
722 fhPhidLam0LowE ->Fill(phi,lambda0/ncells2);
723 }
724 else {
725 fhNCellsLam0HighE ->Fill(ncells,lambda0);
726 fhNCellsLam1HighE ->Fill(ncells,lambda1);
727 fhNCellsDispHighE ->Fill(ncells,disp);
728 fhNCellsdLam0HighE->Fill(ncells,lambda0/ncells2);
729 fhNCellsdLam1HighE->Fill(ncells,lambda1/ncells2);
730 fhNCellsdDispHighE->Fill(ncells,disp/ncells2);
731
732 fhLam1Lam0HighE ->Fill(lambda1,lambda0);
733 fhLam0DispHighE ->Fill(lambda0,disp);
734 fhDispLam1HighE ->Fill(disp,lambda1);
735 fhEtaLam0HighE ->Fill(eta, lambda0);
736 fhPhiLam0HighE ->Fill(phi, lambda0);
737
738 fhdLam1dLam0HighE->Fill(lambda1/ncells2,lambda0/ncells2);
739 fhdLam0dDispHighE->Fill(lambda0/ncells2,disp/ncells2);
740 fhdDispdLam1HighE->Fill(disp/ncells2,lambda1/ncells2);
741 fhEtadLam0HighE ->Fill(eta, lambda0/ncells2);
742 fhPhidLam0HighE ->Fill(phi, lambda0/ncells2);
743 }
744
745 if(IsDataMC()){
3d5d5078 746
f66d95af 747 AliVCaloCells* cells = 0;
748 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
749 else cells = GetPHOSCells();
3d5d5078 750
751 //Fill histograms to check shape of embedded clusters
752 Float_t fraction = 0;
f66d95af 753 if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
754
3d5d5078 755 Float_t clusterE = 0; // recalculate in case corrections applied.
756 Float_t cellE = 0;
757 for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
758 cellE = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
759 clusterE+=cellE;
760 fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
761 }
762
763 //Fraction of total energy due to the embedded signal
764 fraction/=clusterE;
765
766 if(GetDebug() > 1 ) printf("AliAnaPhoton::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
767
768 fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
769
770 } // embedded fraction
771
f66d95af 772 // Get the fraction of the cluster energy that carries the cell with highest energy
773 Int_t absID =-1 ;
774 Float_t maxCellFraction = 0.;
775
776 absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
777
778 // Check the origin and fill histograms
521636d2 779 if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
3d5d5078 780 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
781 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
782 !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)){
783 fhMCELambda0[mcssPhoton] ->Fill(energy, lambda0);
784 fhMCEdLambda0[mcssPhoton] ->Fill(energy, lambda0/ncells2);
785 fhMCELambda1[mcssPhoton] ->Fill(energy, lambda1);
786 fhMCEdLambda1[mcssPhoton] ->Fill(energy, lambda1/ncells2);
787 fhMCEDispersion[mcssPhoton] ->Fill(energy, disp);
788 fhMCEdDispersion[mcssPhoton]->Fill(energy, disp/ncells2);
f66d95af 789 fhMCNCellsE[mcssPhoton] ->Fill(energy, ncells);
790 fhMCMaxCellDiffClusterE[mcssPhoton]->Fill(energy,maxCellFraction);
3d5d5078 791
f66d95af 792 if (energy < 2.){
793 fhMCLambda0vsClusterMaxCellDiffE0[mcssPhoton]->Fill(lambda0, maxCellFraction);
794 fhMCNCellsvsClusterMaxCellDiffE0[mcssPhoton] ->Fill(ncells, maxCellFraction);
795 }
796 else if(energy < 6.){
797 fhMCLambda0vsClusterMaxCellDiffE2[mcssPhoton]->Fill(lambda0, maxCellFraction);
798 fhMCNCellsvsClusterMaxCellDiffE2[mcssPhoton] ->Fill(ncells, maxCellFraction);
799 }
800 else{
801 fhMCLambda0vsClusterMaxCellDiffE6[mcssPhoton]->Fill(lambda0, maxCellFraction);
802 fhMCNCellsvsClusterMaxCellDiffE6[mcssPhoton] ->Fill(ncells, maxCellFraction);
803 }
3d5d5078 804
805 if(!GetReader()->IsEmbeddedClusterSelectionOn()){
806 //Check particle overlaps in cluster
807
808 //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
809 Int_t ancPDG = 0, ancStatus = -1;
810 TLorentzVector momentum; TVector3 prodVertex;
811 Int_t ancLabel = 0;
812 Int_t noverlaps = 1;
813 for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
814 ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
815 if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
816 }
817
818 if(noverlaps == 1){
819 fhMCPhotonELambda0NoOverlap ->Fill(energy, lambda0);
820 fhMCPhotonEdLambda0NoOverlap ->Fill(energy, lambda0/ncells2);
821 }
822 else if(noverlaps == 2){
823 fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0);
824 fhMCPhotonEdLambda0TwoOverlap->Fill(energy, lambda0/ncells2);
825 }
826 else if(noverlaps > 2){
827 fhMCPhotonELambda0NOverlap ->Fill(energy, lambda0);
828 fhMCPhotonEdLambda0NOverlap ->Fill(energy, lambda0/ncells2);
829 }
830 else {
831 printf("AliAnaPhoton::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
832 }
833 }//No embedding
834
835 //Fill histograms to check shape of embedded clusters
836 if(GetReader()->IsEmbeddedClusterSelectionOn()){
837
838 if (fraction > 0.9)
839 {
840 fhEmbedPhotonELambda0FullSignal ->Fill(energy, lambda0);
841 fhEmbedPhotonEdLambda0FullSignal ->Fill(energy, lambda0/ncells2);
842 }
843 else if(fraction > 0.5)
844 {
845 fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0);
846 fhEmbedPhotonEdLambda0MostlySignal->Fill(energy, lambda0/ncells2);
847 }
848 else if(fraction > 0.1)
849 {
850 fhEmbedPhotonELambda0MostlyBkg ->Fill(energy, lambda0);
851 fhEmbedPhotonEdLambda0MostlyBkg ->Fill(energy, lambda0/ncells2);
852 }
853 else
854 {
855 fhEmbedPhotonELambda0FullBkg ->Fill(energy, lambda0);
856 fhEmbedPhotonEdLambda0FullBkg ->Fill(energy, lambda0/ncells2);
857 }
858 } // embedded
859
521636d2 860 }//photon no conversion
861 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)){
3d5d5078 862 fhMCELambda0[mcssElectron] ->Fill(energy, lambda0);
863 fhMCEdLambda0[mcssElectron] ->Fill(energy, lambda0/ncells2);
864 fhMCELambda1[mcssElectron] ->Fill(energy, lambda1);
865 fhMCEdLambda1[mcssElectron] ->Fill(energy, lambda1/ncells2);
866 fhMCEDispersion[mcssElectron] ->Fill(energy, disp);
867 fhMCEdDispersion[mcssElectron]->Fill(energy, disp/ncells2);
f66d95af 868 fhMCNCellsE[mcssElectron] ->Fill(energy, ncells);
869 fhMCMaxCellDiffClusterE[mcssElectron]->Fill(energy,maxCellFraction);
870
871 if (energy < 2.){
872 fhMCLambda0vsClusterMaxCellDiffE0[mcssElectron]->Fill(lambda0, maxCellFraction);
873 fhMCNCellsvsClusterMaxCellDiffE0[mcssElectron] ->Fill(ncells, maxCellFraction);
874 }
875 else if(energy < 6.){
876 fhMCLambda0vsClusterMaxCellDiffE2[mcssElectron]->Fill(lambda0, maxCellFraction);
877 fhMCNCellsvsClusterMaxCellDiffE2[mcssElectron] ->Fill(ncells, maxCellFraction);
878 }
879 else{
880 fhMCLambda0vsClusterMaxCellDiffE6[mcssElectron]->Fill(lambda0, maxCellFraction);
881 fhMCNCellsvsClusterMaxCellDiffE6[mcssElectron] ->Fill(ncells, maxCellFraction);
882 }
521636d2 883 }//electron
3d5d5078 884 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
521636d2 885 GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) ){
3d5d5078 886 fhMCELambda0[mcssConversion] ->Fill(energy, lambda0);
887 fhMCEdLambda0[mcssConversion] ->Fill(energy, lambda0/ncells2);
888 fhMCELambda1[mcssConversion] ->Fill(energy, lambda1);
889 fhMCEdLambda1[mcssConversion] ->Fill(energy, lambda1/ncells2);
890 fhMCEDispersion[mcssConversion] ->Fill(energy, disp);
f66d95af 891 fhMCEdDispersion[mcssConversion]->Fill(energy, disp/ncells2);
892 fhMCNCellsE[mcssConversion] ->Fill(energy, ncells);
893 fhMCMaxCellDiffClusterE[mcssConversion]->Fill(energy,maxCellFraction);
894
895 if (energy < 2.){
896 fhMCLambda0vsClusterMaxCellDiffE0[mcssConversion]->Fill(lambda0, maxCellFraction);
897 fhMCNCellsvsClusterMaxCellDiffE0[mcssConversion] ->Fill(ncells, maxCellFraction);
898 }
899 else if(energy < 6.){
900 fhMCLambda0vsClusterMaxCellDiffE2[mcssConversion]->Fill(lambda0, maxCellFraction);
901 fhMCNCellsvsClusterMaxCellDiffE2[mcssConversion] ->Fill(ncells, maxCellFraction);
902 }
903 else{
904 fhMCLambda0vsClusterMaxCellDiffE6[mcssConversion]->Fill(lambda0, maxCellFraction);
905 fhMCNCellsvsClusterMaxCellDiffE6[mcssConversion] ->Fill(ncells, maxCellFraction);
906 }
3d5d5078 907
521636d2 908 }//conversion photon
909 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) ){
3d5d5078 910 fhMCELambda0[mcssPi0] ->Fill(energy, lambda0);
911 fhMCEdLambda0[mcssPi0] ->Fill(energy, lambda0/ncells2);
912 fhMCELambda1[mcssPi0] ->Fill(energy, lambda1);
913 fhMCEdLambda1[mcssPi0] ->Fill(energy, lambda1/ncells2);
914 fhMCEDispersion[mcssPi0] ->Fill(energy, disp);
915 fhMCEdDispersion[mcssPi0]->Fill(energy, disp/ncells2);
f66d95af 916 fhMCNCellsE[mcssPi0] ->Fill(energy, ncells);
917 fhMCMaxCellDiffClusterE[mcssPi0]->Fill(energy,maxCellFraction);
918
919 if (energy < 2.){
920 fhMCLambda0vsClusterMaxCellDiffE0[mcssPi0]->Fill(lambda0, maxCellFraction);
921 fhMCNCellsvsClusterMaxCellDiffE0[mcssPi0] ->Fill(ncells, maxCellFraction);
922 }
923 else if(energy < 6.){
924 fhMCLambda0vsClusterMaxCellDiffE2[mcssPi0]->Fill(lambda0, maxCellFraction);
925 fhMCNCellsvsClusterMaxCellDiffE2[mcssPi0] ->Fill(ncells, maxCellFraction);
926 }
927 else{
928 fhMCLambda0vsClusterMaxCellDiffE6[mcssPi0]->Fill(lambda0, maxCellFraction);
929 fhMCNCellsvsClusterMaxCellDiffE6[mcssPi0] ->Fill(ncells, maxCellFraction);
930 }
3d5d5078 931
932 //Fill histograms to check shape of embedded clusters
933 if(GetReader()->IsEmbeddedClusterSelectionOn()){
934
935 if (fraction > 0.9)
936 {
937 fhEmbedPi0ELambda0FullSignal ->Fill(energy, lambda0);
938 fhEmbedPi0EdLambda0FullSignal ->Fill(energy, lambda0/ncells2);
939 }
940 else if(fraction > 0.5)
941 {
942 fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0);
943 fhEmbedPi0EdLambda0MostlySignal->Fill(energy, lambda0/ncells2);
944 }
945 else if(fraction > 0.1)
946 {
947 fhEmbedPi0ELambda0MostlyBkg ->Fill(energy, lambda0);
948 fhEmbedPi0EdLambda0MostlyBkg ->Fill(energy, lambda0/ncells2);
949 }
950 else
951 {
952 fhEmbedPi0ELambda0FullBkg ->Fill(energy, lambda0);
953 fhEmbedPi0EdLambda0FullBkg ->Fill(energy, lambda0/ncells2);
954 }
955 } // embedded
956
521636d2 957 }//pi0
3d5d5078 958 else if ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta) ){
959 fhMCELambda0[mcssEta] ->Fill(energy, lambda0);
960 fhMCEdLambda0[mcssEta] ->Fill(energy, lambda0/ncells2);
961 fhMCELambda1[mcssEta] ->Fill(energy, lambda1);
962 fhMCEdLambda1[mcssEta] ->Fill(energy, lambda1/ncells2);
963 fhMCEDispersion[mcssEta] ->Fill(energy, disp);
f66d95af 964 fhMCEdDispersion[mcssEta]->Fill(energy, disp/ncells2);
965 fhMCNCellsE[mcssEta] ->Fill(energy, ncells);
966 fhMCMaxCellDiffClusterE[mcssEta]->Fill(energy,maxCellFraction);
967
968 if (energy < 2.){
969 fhMCLambda0vsClusterMaxCellDiffE0[mcssEta]->Fill(lambda0, maxCellFraction);
970 fhMCNCellsvsClusterMaxCellDiffE0[mcssEta] ->Fill(ncells, maxCellFraction);
971 }
972 else if(energy < 6.){
973 fhMCLambda0vsClusterMaxCellDiffE2[mcssEta]->Fill(lambda0, maxCellFraction);
974 fhMCNCellsvsClusterMaxCellDiffE2[mcssEta] ->Fill(ncells, maxCellFraction);
975 }
976 else{
977 fhMCLambda0vsClusterMaxCellDiffE6[mcssEta]->Fill(lambda0, maxCellFraction);
978 fhMCNCellsvsClusterMaxCellDiffE6[mcssEta] ->Fill(ncells, maxCellFraction);
979 }
980
3d5d5078 981 }//eta
521636d2 982 else {
3d5d5078 983 fhMCELambda0[mcssOther] ->Fill(energy, lambda0);
984 fhMCEdLambda0[mcssOther] ->Fill(energy, lambda0/ncells2);
985 fhMCELambda1[mcssOther] ->Fill(energy, lambda1);
986 fhMCEdLambda1[mcssOther] ->Fill(energy, lambda1/ncells2);
987 fhMCEDispersion[mcssOther] ->Fill(energy, disp);
f66d95af 988 fhMCEdDispersion[mcssOther]->Fill(energy, disp/ncells2);
989 fhMCNCellsE[mcssOther] ->Fill(energy, ncells);
990 fhMCMaxCellDiffClusterE[mcssOther]->Fill(energy,maxCellFraction);
991
992 if (energy < 2.){
993 fhMCLambda0vsClusterMaxCellDiffE0[mcssOther]->Fill(lambda0, maxCellFraction);
994 fhMCNCellsvsClusterMaxCellDiffE0[mcssOther] ->Fill(ncells, maxCellFraction);
995 }
996 else if(energy < 6.){
997 fhMCLambda0vsClusterMaxCellDiffE2[mcssOther]->Fill(lambda0, maxCellFraction);
998 fhMCNCellsvsClusterMaxCellDiffE2[mcssOther] ->Fill(ncells, maxCellFraction);
999 }
1000 else{
1001 fhMCLambda0vsClusterMaxCellDiffE6[mcssOther]->Fill(lambda0, maxCellFraction);
1002 fhMCNCellsvsClusterMaxCellDiffE6[mcssOther] ->Fill(ncells, maxCellFraction);
1003 }
1004
521636d2 1005 }//other particles
1006
1007 }//MC data
1008
1009}
1010
0c1383b5 1011//________________________________________________________________________
1012TObjString * AliAnaPhoton::GetAnalysisCuts()
1013{
1014 //Save parameters used for analysis
1015 TString parList ; //this will be list of parameters used for this analysis.
5ae09196 1016 const Int_t buffersize = 255;
1017 char onePar[buffersize] ;
0c1383b5 1018
5ae09196 1019 snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
0c1383b5 1020 parList+=onePar ;
5ae09196 1021 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
0c1383b5 1022 parList+=onePar ;
5ae09196 1023 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
0c1383b5 1024 parList+=onePar ;
5ae09196 1025 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
0c1383b5 1026 parList+=onePar ;
5ae09196 1027 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
0c1383b5 1028 parList+=onePar ;
5ae09196 1029 snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
0c1383b5 1030 parList+=onePar ;
41121cfe 1031 snprintf(onePar,buffersize,"Conversion Selection: fConvAsymCut %1.2f, fConvDEtaCut %1.2f fConvDPhiCut (%1.2f,%1.2f)\n",
1032 fConvAsymCut, fConvDEtaCut, fConvDPhiMinCut, fConvDPhiMaxCut) ;
6175da48 1033 parList+=onePar ;
0c1383b5 1034
1035 //Get parameters set in base class.
1036 parList += GetBaseParametersList() ;
1037
1038 //Get parameters set in PID class.
1039 parList += GetCaloPID()->GetPIDParametersList() ;
1040
1041 //Get parameters set in FiducialCut class (not available yet)
1042 //parlist += GetFidCut()->GetFidCutParametersList()
1043
1044 return new TObjString(parList) ;
1045}
1046
1c5acb87 1047//________________________________________________________________________
1048TList * AliAnaPhoton::GetCreateOutputObjects()
1049{
477d6cee 1050 // Create histograms to be saved in output file and
1051 // store them in outputContainer
1052 TList * outputContainer = new TList() ;
1053 outputContainer->SetName("PhotonHistos") ;
4a745797 1054
521636d2 1055 Int_t nptbins = GetHistoPtBins(); Float_t ptmax = GetHistoPtMax(); Float_t ptmin = GetHistoPtMin();
1056 Int_t nphibins = GetHistoPhiBins(); Float_t phimax = GetHistoPhiMax(); Float_t phimin = GetHistoPhiMin();
1057 Int_t netabins = GetHistoEtaBins(); Float_t etamax = GetHistoEtaMax(); Float_t etamin = GetHistoEtaMin();
1058 Int_t ssbins = GetHistoShowerShapeBins(); Float_t ssmax = GetHistoShowerShapeMax(); Float_t ssmin = GetHistoShowerShapeMin();
e1e62b89 1059 Int_t nbins = GetHistoNClusterCellBins(); Int_t nmax = GetHistoNClusterCellMax(); Int_t nmin = GetHistoNClusterCellMin();
521636d2 1060
e1e62b89 1061 fhNCellsE = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", nptbins,ptmin,ptmax, nbins,nmin,nmax);
c4a7d28a 1062 fhNCellsE->SetXTitle("E (GeV)");
1063 fhNCellsE->SetYTitle("# of cells in cluster");
1064 outputContainer->Add(fhNCellsE);
6175da48 1065
f66d95af 1066 fhMaxCellDiffClusterE = new TH2F ("hMaxCellDiffClusterE","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1067 nptbins,ptmin,ptmax, 500,0,1.);
1068 fhMaxCellDiffClusterE->SetXTitle("E_{cluster} (GeV) ");
1069 fhMaxCellDiffClusterE->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1070 outputContainer->Add(fhMaxCellDiffClusterE);
1071
20218aea 1072 fhEPhoton = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
1073 fhEPhoton->SetYTitle("N");
1074 fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
1075 outputContainer->Add(fhEPhoton) ;
1076
1077 fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax);
477d6cee 1078 fhPtPhoton->SetYTitle("N");
1079 fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
1080 outputContainer->Add(fhPtPhoton) ;
1081
1082 fhPhiPhoton = new TH2F
20218aea 1083 ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
6175da48 1084 fhPhiPhoton->SetYTitle("#phi (rad)");
477d6cee 1085 fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
1086 outputContainer->Add(fhPhiPhoton) ;
1087
1088 fhEtaPhoton = new TH2F
20218aea 1089 ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 1090 fhEtaPhoton->SetYTitle("#eta");
1091 fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
1092 outputContainer->Add(fhEtaPhoton) ;
1093
6175da48 1094 fhEtaPhiPhoton = new TH2F
1095 ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
1096 fhEtaPhiPhoton->SetYTitle("#phi (rad)");
1097 fhEtaPhiPhoton->SetXTitle("#eta");
1098 outputContainer->Add(fhEtaPhiPhoton) ;
20218aea 1099 if(GetMinPt() < 0.5){
1100 fhEtaPhi05Photon = new TH2F
1101 ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
1102 fhEtaPhi05Photon->SetYTitle("#phi (rad)");
1103 fhEtaPhi05Photon->SetXTitle("#eta");
1104 outputContainer->Add(fhEtaPhi05Photon) ;
1105 }
6175da48 1106
1107 //Conversion
20218aea 1108 if(fCheckConversion){
521636d2 1109
1110 fhEtaPhiConversion = new TH2F
1111 ("hEtaPhiConversion","#eta vs #phi for converted clusters",netabins,etamin,etamax,nphibins,phimin,phimax);
1112 fhEtaPhiConversion->SetYTitle("#phi (rad)");
1113 fhEtaPhiConversion->SetXTitle("#eta");
1114 outputContainer->Add(fhEtaPhiConversion) ;
1115 if(GetMinPt() < 0.5){
1116 fhEtaPhi05Conversion = new TH2F
1117 ("hEtaPhi05Conversion","#eta vs #phi, E > 0.5, for converted clusters",netabins,etamin,etamax,nphibins,phimin,phimax);
1118 fhEtaPhi05Conversion->SetYTitle("#phi (rad)");
1119 fhEtaPhi05Conversion->SetXTitle("#eta");
1120 outputContainer->Add(fhEtaPhi05Conversion) ;
1121 }
1122
20218aea 1123 fhPtPhotonConv = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax);
1124 fhPtPhotonConv->SetYTitle("N");
1125 fhPtPhotonConv->SetXTitle("p_{T #gamma}(GeV/c)");
1126 outputContainer->Add(fhPtPhotonConv) ;
1127
1128 fhEtaPhiPhotonConv = new TH2F
1129 ("hEtaPhiPhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
1130 fhEtaPhiPhotonConv->SetYTitle("#phi (rad)");
1131 fhEtaPhiPhotonConv->SetXTitle("#eta");
1132 outputContainer->Add(fhEtaPhiPhotonConv) ;
1133 if(GetMinPt() < 0.5){
1134 fhEtaPhi05PhotonConv = new TH2F
1135 ("hEtaPhi05PhotonConv","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
1136 fhEtaPhi05PhotonConv->SetYTitle("#phi (rad)");
1137 fhEtaPhi05PhotonConv->SetXTitle("#eta");
1138 outputContainer->Add(fhEtaPhi05PhotonConv) ;
1139 }
1140
1141 fhConvDeltaEta = new TH2F
1142 ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins*2,-0.5,0.5);
1143 fhConvDeltaEta->SetYTitle("#Delta #eta");
1144 fhConvDeltaEta->SetXTitle("Pair Mass (GeV/c^2)");
1145 outputContainer->Add(fhConvDeltaEta) ;
1146
1147 fhConvDeltaPhi = new TH2F
1148 ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins*2,-0.5,0.5);
1149 fhConvDeltaPhi->SetYTitle("#Delta #phi");
1150 fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/c^2)");
1151 outputContainer->Add(fhConvDeltaPhi) ;
1152
1153 fhConvDeltaEtaPhi = new TH2F
1154 ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5);
1155 fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
1156 fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
1157 outputContainer->Add(fhConvDeltaEtaPhi) ;
1158
1159 fhConvAsym = new TH2F
1160 ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1);
1161 fhConvAsym->SetYTitle("Asymmetry");
1162 fhConvAsym->SetXTitle("Pair Mass (GeV/c^2)");
1163 outputContainer->Add(fhConvAsym) ;
1164
1165 fhConvPt = new TH2F
1166 ("hConvPt","p_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.);
1167 fhConvPt->SetYTitle("Pair p_{T} (GeV/c)");
1168 fhConvPt->SetXTitle("Pair Mass (GeV/c^2)");
1169 outputContainer->Add(fhConvPt) ;
c4a7d28a 1170
1171 fhConvDistEta = new TH2F
1172 ("hConvDistEta","distance to conversion vertex",100,-0.7,0.7,100,0.,5.);
1173 fhConvDistEta->SetXTitle("#eta");
1174 fhConvDistEta->SetYTitle(" distance (m)");
1175 outputContainer->Add(fhConvDistEta) ;
1176
1177 fhConvDistEn = new TH2F
1178 ("hConvDistEn","distance to conversion vertex",nptbins,ptmin,ptmax,100,0.,5.);
1179 fhConvDistEn->SetXTitle("E (GeV)");
1180 fhConvDistEn->SetYTitle(" distance (m)");
1181 outputContainer->Add(fhConvDistEn) ;
1182
1183 fhConvDistMass = new TH2F
1184 ("hConvDistMass","distance to conversion vertex",100,0,fMassCut,100,0.,5.);
1185 fhConvDistMass->SetXTitle("m (GeV/c^2)");
1186 fhConvDistMass->SetYTitle(" distance (m)");
1187 outputContainer->Add(fhConvDistMass) ;
1188
1189 fhConvDistEtaCutEta = new TH2F
1190 ("hConvDistEtaCutEta","distance to conversion vertex, dEta < 0.05",100,-0.7,0.7,100,0.,5.);
1191 fhConvDistEtaCutEta->SetXTitle("#eta");
1192 fhConvDistEtaCutEta->SetYTitle(" distance (m)");
1193 outputContainer->Add(fhConvDistEtaCutEta) ;
1194
1195 fhConvDistEnCutEta = new TH2F
1196 ("hConvDistEnCutEta","distance to conversion vertex, dEta < 0.05",nptbins,ptmin,ptmax,100,0.,5.);
1197 fhConvDistEnCutEta->SetXTitle("E (GeV)");
1198 fhConvDistEnCutEta->SetYTitle(" distance (m)");
1199 outputContainer->Add(fhConvDistEnCutEta) ;
1200
1201 fhConvDistMassCutEta = new TH2F
1202 ("hConvDistMassCutEta","distance to conversion vertex, dEta < 0.05",100,0,fMassCut,100,0.,5.);
1203 fhConvDistMassCutEta->SetXTitle("m (GeV/c^2)");
1204 fhConvDistMassCutEta->SetYTitle(" distance (m)");
1205 outputContainer->Add(fhConvDistMassCutEta) ;
1206
1207 fhConvDistEtaCutMass = new TH2F
1208 ("hConvDistEtaCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",100,-0.7,0.7,100,0.,5.);
1209 fhConvDistEtaCutMass->SetXTitle("#eta");
1210 fhConvDistEtaCutMass->SetYTitle(" distance (m)");
1211 outputContainer->Add(fhConvDistEtaCutMass) ;
1212
1213 fhConvDistEnCutMass = new TH2F
1214 ("hConvDistEnCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",nptbins,ptmin,ptmax,100,0.,5.);
1215 fhConvDistEnCutMass->SetXTitle("E (GeV)");
1216 fhConvDistEnCutMass->SetYTitle(" distance (m)");
1217 outputContainer->Add(fhConvDistEnCutMass) ;
1218
1219 fhConvDistEtaCutAsy = new TH2F
1220 ("hConvDistEtaCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",100,-0.7,0.7,100,0.,5.);
1221 fhConvDistEtaCutAsy->SetXTitle("#eta");
1222 fhConvDistEtaCutAsy->SetYTitle(" distance (m)");
1223 outputContainer->Add(fhConvDistEtaCutAsy) ;
1224
1225 fhConvDistEnCutAsy = new TH2F
1226 ("hConvDistEnCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",nptbins,ptmin,ptmax,100,0.,5.);
1227 fhConvDistEnCutAsy->SetXTitle("E (GeV)");
1228 fhConvDistEnCutAsy->SetYTitle(" distance (m)");
1229 outputContainer->Add(fhConvDistEnCutAsy) ;
1230
521636d2 1231 } // check conversion
1232
1233
1234 //Shower shape
1235 if(fFillSSHistograms){
1236
1237 fhLam0E = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1238 fhLam0E->SetYTitle("#lambda_{0}^{2}");
1239 fhLam0E->SetXTitle("E (GeV)");
1240 outputContainer->Add(fhLam0E);
1241
1242 fhLam1E = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1243 fhLam1E->SetYTitle("#lambda_{1}^{2}");
1244 fhLam1E->SetXTitle("E (GeV)");
1245 outputContainer->Add(fhLam1E);
1246
1247 fhDispE = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1248 fhDispE->SetYTitle("D^{2}");
1249 fhDispE->SetXTitle("E (GeV) ");
1250 outputContainer->Add(fhDispE);
1251
1252 fhdLam0E = new TH2F ("hdLam0E","#lambda_{0}^{2}/N_{cells}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1253 fhdLam0E->SetYTitle("d#lambda_{0}^{2}");
1254 fhdLam0E->SetXTitle("E (GeV)");
1255 outputContainer->Add(fhdLam0E);
1256
1257 fhdLam1E = new TH2F ("hdLam1E","#lambda_{1}^{2}/N_{cells}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1258 fhdLam1E->SetYTitle("d#lambda_{1}^{2}");
1259 fhdLam1E->SetXTitle("E (GeV)");
1260 outputContainer->Add(fhdLam1E);
1261
1262 fhdDispE = new TH2F ("hdDispE"," dispersion^{2}/N_{cells}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1263 fhdDispE->SetYTitle("dD^{2}");
1264 fhdDispE->SetXTitle("E (GeV) ");
1265 outputContainer->Add(fhdDispE);
1266
1267
1268 if(fCalorimeter == "EMCAL"){
1269 fhLam0ETRD = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1270 fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
1271 fhLam0ETRD->SetXTitle("E (GeV)");
1272 outputContainer->Add(fhLam0ETRD);
1273
1274 fhLam1ETRD = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1275 fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
1276 fhLam1ETRD->SetXTitle("E (GeV)");
1277 outputContainer->Add(fhLam1ETRD);
1278
1279 fhDispETRD = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1280 fhDispETRD->SetYTitle("Dispersion^{2}");
1281 fhDispETRD->SetXTitle("E (GeV) ");
1282 outputContainer->Add(fhDispETRD);
1283
1284 fhdLam0ETRD = new TH2F ("hdLam0ETRD","#lambda_{0}^{2}/N_{cells}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1285 fhdLam0ETRD->SetYTitle("d#lambda_{0}^{2}");
1286 fhdLam0ETRD->SetXTitle("E (GeV)");
1287 outputContainer->Add(fhdLam0ETRD);
1288
1289 fhdLam1ETRD = new TH2F ("hdLam1ETRD","#lambda_{1}^{2}/N_{cells}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1290 fhdLam1ETRD->SetYTitle("d#lambda_{1}^{2}");
1291 fhdLam1ETRD->SetXTitle("E (GeV)");
1292 outputContainer->Add(fhdLam1ETRD);
1293
1294 fhdDispETRD = new TH2F ("hdDispETRD"," dispersion^{2}/N_{cells}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1295 fhdDispETRD->SetYTitle("dD^{2}");
1296 fhdDispETRD->SetXTitle("E (GeV) ");
1297 outputContainer->Add(fhdDispETRD);
1298
1299 }
1300
1301 fhNCellsLam0LowE = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1302 fhNCellsLam0LowE->SetXTitle("N_{cells}");
1303 fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
1304 outputContainer->Add(fhNCellsLam0LowE);
1305
1306 fhNCellsLam0HighE = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1307 fhNCellsLam0HighE->SetXTitle("N_{cells}");
1308 fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
1309 outputContainer->Add(fhNCellsLam0HighE);
1310
1311 fhNCellsLam1LowE = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1312 fhNCellsLam1LowE->SetXTitle("N_{cells}");
1313 fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
1314 outputContainer->Add(fhNCellsLam1LowE);
1315
1316 fhNCellsLam1HighE = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1317 fhNCellsLam1HighE->SetXTitle("N_{cells}");
1318 fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
1319 outputContainer->Add(fhNCellsLam1HighE);
1320
1321 fhNCellsDispLowE = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1322 fhNCellsDispLowE->SetXTitle("N_{cells}");
1323 fhNCellsDispLowE->SetYTitle("D^{2}");
1324 outputContainer->Add(fhNCellsDispLowE);
1325
1326 fhNCellsDispHighE = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1327 fhNCellsDispHighE->SetXTitle("N_{cells}");
1328 fhNCellsDispHighE->SetYTitle("D^{2}");
1329 outputContainer->Add(fhNCellsDispHighE);
1330
1331
1332 fhNCellsdLam0LowE = new TH2F ("hNCellsdLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1333 fhNCellsdLam0LowE->SetXTitle("N_{cells}");
1334 fhNCellsdLam0LowE->SetYTitle("#lambda_{0}^{2}");
1335 outputContainer->Add(fhNCellsdLam0LowE);
1336
1337 fhNCellsdLam0HighE = new TH2F ("hNCellsdLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}/N_{cells}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1338 fhNCellsdLam0HighE->SetXTitle("N_{cells}");
1339 fhNCellsdLam0HighE->SetYTitle("#lambda_{0}^{2}");
1340 outputContainer->Add(fhNCellsdLam0HighE);
1341
1342 fhNCellsdLam1LowE = new TH2F ("hNCellsdLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1343 fhNCellsdLam1LowE->SetXTitle("N_{cells}");
1344 fhNCellsdLam1LowE->SetYTitle("#lambda_{0}^{2}");
1345 outputContainer->Add(fhNCellsdLam1LowE);
1346
1347 fhNCellsdLam1HighE = new TH2F ("hNCellsdLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}/N_{cells}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1348 fhNCellsdLam1HighE->SetXTitle("N_{cells}");
1349 fhNCellsdLam1HighE->SetYTitle("#lambda_{0}^{2}");
1350 outputContainer->Add(fhNCellsdLam1HighE);
1351
1352 fhNCellsdDispLowE = new TH2F ("hNCellsdDispLowE","N_{cells} in cluster vs dispersion^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1353 fhNCellsdDispLowE->SetXTitle("N_{cells}");
1354 fhNCellsdDispLowE->SetYTitle("D^{2}");
1355 outputContainer->Add(fhNCellsdDispLowE);
1356
1357 fhNCellsdDispHighE = new TH2F ("hNCellsdDispHighE","N_{cells} in cluster vs dispersion^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1358 fhNCellsdDispHighE->SetXTitle("N_{cells}");
1359 fhNCellsdDispHighE->SetYTitle("D^{2}");
1360 outputContainer->Add(fhNCellsdDispHighE);
1361
1362
1363 fhEtaLam0LowE = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1364 fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
1365 fhEtaLam0LowE->SetXTitle("#eta");
1366 outputContainer->Add(fhEtaLam0LowE);
1367
1368 fhPhiLam0LowE = new TH2F ("hPhiLam0LowE","#phi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1369 fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
1370 fhPhiLam0LowE->SetXTitle("#phi");
1371 outputContainer->Add(fhPhiLam0LowE);
1372
1373 fhEtaLam0HighE = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, E > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1374 fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
1375 fhEtaLam0HighE->SetXTitle("#eta");
1376 outputContainer->Add(fhEtaLam0HighE);
1377
1378 fhPhiLam0HighE = new TH2F ("hPhiLam0HighE","#phi vs #lambda_{0}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1379 fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
1380 fhPhiLam0HighE->SetXTitle("#phi");
1381 outputContainer->Add(fhPhiLam0HighE);
1382
1383 fhLam1Lam0LowE = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1384 fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
1385 fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
1386 outputContainer->Add(fhLam1Lam0LowE);
1387
1388 fhLam1Lam0HighE = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1389 fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
1390 fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
1391 outputContainer->Add(fhLam1Lam0HighE);
1392
1393 fhLam0DispLowE = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1394 fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
1395 fhLam0DispLowE->SetYTitle("D^{2}");
1396 outputContainer->Add(fhLam0DispLowE);
1397
1398 fhLam0DispHighE = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1399 fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
1400 fhLam0DispHighE->SetYTitle("D^{2}");
1401 outputContainer->Add(fhLam0DispHighE);
1402
1403 fhDispLam1LowE = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1404 fhDispLam1LowE->SetXTitle("D^{2}");
1405 fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
1406 outputContainer->Add(fhDispLam1LowE);
1407
1408 fhDispLam1HighE = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of E > 2 GeV", ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1409 fhDispLam1HighE->SetXTitle("D^{2}");
1410 fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
1411 outputContainer->Add(fhDispLam1HighE);
1412
1413
1414
1415 fhEtadLam0LowE = new TH2F ("hEtadLam0LowE","#eta vs #lambda_{0}^{2}/N_{cells}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax/50);
1416 fhEtadLam0LowE->SetYTitle("d#lambda_{0}^{2}");
1417 fhEtadLam0LowE->SetXTitle("#eta");
1418 outputContainer->Add(fhEtadLam0LowE);
1419
1420 fhPhidLam0LowE = new TH2F ("hPhidLam0LowE","#phi vs #lambda_{0}^{2}/N_{cells}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax/50);
1421 fhPhidLam0LowE->SetYTitle("d#lambda_{0}^{2}");
1422 fhPhidLam0LowE->SetXTitle("#phi");
1423 outputContainer->Add(fhPhidLam0LowE);
1424
1425 fhEtadLam0HighE = new TH2F ("hEtadLam0HighE","#eta vs #lambda_{0}^{2}/N_{cells}^{2}", netabins,etamin,etamax, ssbins,ssmin,ssmax/50);
1426 fhEtadLam0HighE->SetYTitle("d#lambda_{0}^{2}");
1427 fhEtadLam0HighE->SetXTitle("#eta");
1428 outputContainer->Add(fhEtadLam0HighE);
1429
1430 fhPhidLam0HighE = new TH2F ("hPhidLam0HighE","#phi vs #lambda_{0}^{2}/N_{cells}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax/50);
1431 fhPhidLam0HighE->SetYTitle("d#lambda_{0}^{2}");
1432 fhPhidLam0HighE->SetXTitle("#phi");
1433 outputContainer->Add(fhPhidLam0HighE);
1434
1435 fhdLam1dLam0LowE = new TH2F ("hdLam1dLam0LowE","#lambda_{0}^{2}/N_{cells}^{2} vs #lambda_{1}^{2}/N_{cells}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1436 fhdLam1dLam0LowE->SetYTitle("d#lambda_{0}^{2}");
1437 fhdLam1dLam0LowE->SetXTitle("d#lambda_{1}^{2}");
1438 outputContainer->Add(fhdLam1dLam0LowE);
1439
1440 fhdLam1dLam0HighE = new TH2F ("hdLam1dLam0HighE","#lambda_{0}^{2}/N_{cells}^{2}vs #lambda_{1}^{2}/N_{cells}^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1441 fhdLam1dLam0HighE->SetYTitle("d#lambda_{0}^{2}");
1442 fhdLam1dLam0HighE->SetXTitle("d#lambda_{1}^{2}");
1443 outputContainer->Add(fhdLam1dLam0HighE);
1444
1445 fhdLam0dDispLowE = new TH2F ("hdLam0dDispLowE","#lambda_{0}^{2}/N_{cells}^{2} vs dispersion^{2}/N_{cells}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1446 fhdLam0dDispLowE->SetXTitle("d#lambda_{0}^{2}");
1447 fhdLam0dDispLowE->SetYTitle("dD^{2}");
1448 outputContainer->Add(fhdLam0dDispLowE);
1449
1450 fhdLam0dDispHighE = new TH2F ("hdLam0dDispHighE","#lambda_{0}^{2}/N_{cells}^{2} vs dispersion^{2}/N_{cells}^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1451 fhdLam0dDispHighE->SetXTitle("d#lambda_{0}^{2}");
1452 fhdLam0dDispHighE->SetYTitle("dD^{2}");
1453 outputContainer->Add(fhdLam0dDispHighE);
1454
1455 fhdDispdLam1LowE = new TH2F ("hdDispddLam1LowE","Dispersion^{2}/N_{cells}^{2} vs #lambda_{1}^{2}/N_{cells}^{2} in cluster of E < 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1456 fhdDispdLam1LowE->SetXTitle("dD^{2}");
1457 fhdDispdLam1LowE->SetYTitle("d#lambda_{1}^{2}");
1458 outputContainer->Add(fhdDispdLam1LowE);
1459
1460 fhdDispdLam1HighE = new TH2F ("hdDispdLam1HighE","Dispersion^{2}/N_{cells}^{2} vs #lambda_{1^{2}}/N_{cells}^{2} in cluster of E > 2 GeV", ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1461 fhdDispdLam1HighE->SetXTitle("dD^{2}");
1462 fhdDispdLam1HighE->SetYTitle("d#lambda_{1}^{2}");
1463 outputContainer->Add(fhdDispdLam1HighE);
3d5d5078 1464
521636d2 1465
1466 } // Shower shape
1467
6175da48 1468
477d6cee 1469 if(IsDataMC()){
123fc3bd 1470 fhDeltaE = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
1471 fhDeltaE->SetXTitle("#Delta E (GeV)");
1472 outputContainer->Add(fhDeltaE);
1473
1474 fhDeltaPt = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
1475 fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
1476 outputContainer->Add(fhDeltaPt);
1477
1478 fhRatioE = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
1479 fhRatioE->SetXTitle("E_{reco}/E_{gen}");
1480 outputContainer->Add(fhRatioE);
477d6cee 1481
123fc3bd 1482 fhRatioPt = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
1483 fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
1484 outputContainer->Add(fhRatioPt);
1485
1486 fh2E = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 1487 fh2E->SetXTitle("E_{rec} (GeV)");
1488 fh2E->SetYTitle("E_{gen} (GeV)");
123fc3bd 1489 outputContainer->Add(fh2E);
1490
1491 fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 1492 fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
1493 fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
123fc3bd 1494 outputContainer->Add(fh2Pt);
1495
f66d95af 1496 TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
1497 "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P",
1498 "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}","String" } ;
3d5d5078 1499
f66d95af 1500 TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
1501 "Conversion", "Hadron", "AntiNeutron","AntiProton",
1502 "PhotonPrompt","PhotonFragmentation","PhotonISR","String" } ;
521636d2 1503
f66d95af 1504 for(Int_t i = 0; i < fNOriginHistograms; i++){
521636d2 1505
3d5d5078 1506 fhMCE[i] = new TH1F(Form("hE_MC%s",pname[i].Data()),
521636d2 1507 Form("cluster from %s : E ",ptype[i].Data()),
1508 nptbins,ptmin,ptmax);
3d5d5078 1509 fhMCE[i]->SetXTitle("E (GeV)");
1510 outputContainer->Add(fhMCE[i]) ;
521636d2 1511
1512 fhPtMC[i] = new TH1F(Form("hPt_MC%s",pname[i].Data()),
1513 Form("cluster from %s : p_{T} ",ptype[i].Data()),
1514 nptbins,ptmin,ptmax);
1515 fhPtMC[i]->SetXTitle("p_{T} (GeV/c)");
1516 outputContainer->Add(fhPtMC[i]) ;
1517
1518 fhEtaMC[i] = new TH2F(Form("hEta_MC%s",pname[i].Data()),
1519 Form("cluster from %s : #eta ",ptype[i].Data()),
1520 nptbins,ptmin,ptmax,netabins,etamin,etamax);
1521 fhEtaMC[i]->SetYTitle("#eta");
1522 fhEtaMC[i]->SetXTitle("E (GeV)");
1523 outputContainer->Add(fhEtaMC[i]) ;
1524
1525 fhPhiMC[i] = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
1526 Form("cluster from %s : #phi ",ptype[i].Data()),
1527 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1528 fhPhiMC[i]->SetYTitle("#phi (rad)");
1529 fhPhiMC[i]->SetXTitle("E (GeV)");
1530 outputContainer->Add(fhPhiMC[i]) ;
1531
1532 }
3d5d5078 1533
f66d95af 1534 TString pptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}","hadron?",
1535 "#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}"} ;
1536
1537 TString ppname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Hadron",
1538 "PhotonPrompt","PhotonFragmentation","PhotonISR"} ;
1539
1540 for(Int_t i = 0; i < fNPrimaryHistograms; i++){
1541 fhEPrimMC[i] = new TH1F(Form("hEPrim_MC%s",ppname[i].Data()),
1542 Form("primary photon %s : E ",pptype[i].Data()),
3d5d5078 1543 nptbins,ptmin,ptmax);
1544 fhEPrimMC[i]->SetXTitle("E (GeV)");
1545 outputContainer->Add(fhEPrimMC[i]) ;
1546
f66d95af 1547 fhPtPrimMC[i] = new TH1F(Form("hPtPrim_MC%s",ppname[i].Data()),
1548 Form("primary photon %s : p_{T} ",pptype[i].Data()),
3d5d5078 1549 nptbins,ptmin,ptmax);
1550 fhPtPrimMC[i]->SetXTitle("p_{T} (GeV/c)");
1551 outputContainer->Add(fhPtPrimMC[i]) ;
1552
f66d95af 1553 fhYPrimMC[i] = new TH2F(Form("hYPrim_MC%s",ppname[i].Data()),
1554 Form("primary photon %s : Rapidity ",pptype[i].Data()),
3d5d5078 1555 nptbins,ptmin,ptmax,800,-8,8);
1556 fhYPrimMC[i]->SetYTitle("Rapidity");
1557 fhYPrimMC[i]->SetXTitle("E (GeV)");
1558 outputContainer->Add(fhYPrimMC[i]) ;
1559
f66d95af 1560 fhPhiPrimMC[i] = new TH2F(Form("hPhiPrim_MC%s",ppname[i].Data()),
1561 Form("primary photon %s : #phi ",pptype[i].Data()),
3d5d5078 1562 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1563 fhPhiPrimMC[i]->SetYTitle("#phi (rad)");
1564 fhPhiPrimMC[i]->SetXTitle("E (GeV)");
1565 outputContainer->Add(fhPhiPrimMC[i]) ;
1566
1567
f66d95af 1568 fhEPrimMCAcc[i] = new TH1F(Form("hEPrimAcc_MC%s",ppname[i].Data()),
1569 Form("primary photon %s in acceptance: E ",pptype[i].Data()),
3d5d5078 1570 nptbins,ptmin,ptmax);
1571 fhEPrimMCAcc[i]->SetXTitle("E (GeV)");
1572 outputContainer->Add(fhEPrimMCAcc[i]) ;
1573
f66d95af 1574 fhPtPrimMCAcc[i] = new TH1F(Form("hPtPrimAcc_MC%s",ppname[i].Data()),
1575 Form("primary photon %s in acceptance: p_{T} ",pptype[i].Data()),
3d5d5078 1576 nptbins,ptmin,ptmax);
1577 fhPtPrimMCAcc[i]->SetXTitle("p_{T} (GeV/c)");
1578 outputContainer->Add(fhPtPrimMCAcc[i]) ;
1579
f66d95af 1580 fhYPrimMCAcc[i] = new TH2F(Form("hYPrimAcc_MC%s",ppname[i].Data()),
1581 Form("primary photon %s in acceptance: Rapidity ",pptype[i].Data()),
3d5d5078 1582 nptbins,ptmin,ptmax,100,-1,1);
1583 fhYPrimMCAcc[i]->SetYTitle("Rapidity");
1584 fhYPrimMCAcc[i]->SetXTitle("E (GeV)");
1585 outputContainer->Add(fhYPrimMCAcc[i]) ;
1586
f66d95af 1587 fhPhiPrimMCAcc[i] = new TH2F(Form("hPhiPrimAcc_MC%s",ppname[i].Data()),
1588 Form("primary photon %s in acceptance: #phi ",pptype[i].Data()),
3d5d5078 1589 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1590 fhPhiPrimMCAcc[i]->SetYTitle("#phi (rad)");
1591 fhPhiPrimMCAcc[i]->SetXTitle("E (GeV)");
1592 outputContainer->Add(fhPhiPrimMCAcc[i]) ;
1593
1594 }
1595
20218aea 1596 if(fCheckConversion){
1597 fhPtConversionTagged = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
1598 fhPtConversionTagged->SetYTitle("N");
1599 fhPtConversionTagged->SetXTitle("p_{T #gamma}(GeV/c)");
1600 outputContainer->Add(fhPtConversionTagged) ;
1601
1602
1603 fhPtAntiNeutronTagged = new TH1F("hPtMCAntiNeutronTagged","Number of AntiNeutron id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
1604 fhPtAntiNeutronTagged->SetYTitle("N");
1605 fhPtAntiNeutronTagged->SetXTitle("p_{T #gamma}(GeV/c)");
1606 outputContainer->Add(fhPtAntiNeutronTagged) ;
1607
1608 fhPtAntiProtonTagged = new TH1F("hPtMCAntiProtonTagged","Number of AntiProton id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
1609 fhPtAntiProtonTagged->SetYTitle("N");
1610 fhPtAntiProtonTagged->SetXTitle("p_{T #gamma}(GeV/c)");
1611 outputContainer->Add(fhPtAntiProtonTagged) ;
1612
1613 fhPtUnknownTagged = new TH1F("hPtMCUnknownTagged","Number of Unknown id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
1614 fhPtUnknownTagged->SetYTitle("N");
1615 fhPtUnknownTagged->SetXTitle("p_{T #gamma}(GeV/c)");
1616 outputContainer->Add(fhPtUnknownTagged) ;
1617
1618 fhConvDeltaEtaMCConversion = new TH2F
1619 ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5);
1620 fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
1621 fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
1622 outputContainer->Add(fhConvDeltaEtaMCConversion) ;
1623
1624 fhConvDeltaPhiMCConversion = new TH2F
1625 ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5);
1626 fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
1627 fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
1628 outputContainer->Add(fhConvDeltaPhiMCConversion) ;
1629
1630 fhConvDeltaEtaPhiMCConversion = new TH2F
1631 ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5);
1632 fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
1633 fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
1634 outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
1635
1636 fhConvAsymMCConversion = new TH2F
1637 ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1);
1638 fhConvAsymMCConversion->SetYTitle("Asymmetry");
1639 fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
1640 outputContainer->Add(fhConvAsymMCConversion) ;
1641
1642 fhConvPtMCConversion = new TH2F
1643 ("hConvPtMCConversion","p_{T} of selected conversion pairs from real conversions",100,0,fMassCut,100,0.,10.);
1644 fhConvPtMCConversion->SetYTitle("Pair p_{T} (GeV/c)");
1645 fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
1646 outputContainer->Add(fhConvPtMCConversion) ;
1647
1648 fhConvDispersionMCConversion = new TH2F
1649 ("hConvDispersionMCConversion","p_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.);
1650 fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
1651 fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
1652 outputContainer->Add(fhConvDispersionMCConversion) ;
1653
1654 fhConvM02MCConversion = new TH2F
1655 ("hConvM02MCConversion","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
1656 fhConvM02MCConversion->SetYTitle("M02 cluster 1");
1657 fhConvM02MCConversion->SetXTitle("M02 cluster 2");
1658 outputContainer->Add(fhConvM02MCConversion) ;
1659
1660 fhConvDeltaEtaMCAntiNeutron = new TH2F
1661 ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5);
1662 fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
1663 fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
1664 outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
1665
1666 fhConvDeltaPhiMCAntiNeutron = new TH2F
1667 ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5);
1668 fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
1669 fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
1670 outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
1671
1672 fhConvDeltaEtaPhiMCAntiNeutron = new TH2F
1673 ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
1674 fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
1675 fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
1676 outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;
1677
1678 fhConvAsymMCAntiNeutron = new TH2F
1679 ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1);
1680 fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
1681 fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
1682 outputContainer->Add(fhConvAsymMCAntiNeutron) ;
1683
1684 fhConvPtMCAntiNeutron = new TH2F
1685 ("hConvPtMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0.,10.);
1686 fhConvPtMCAntiNeutron->SetYTitle("Pair p_{T} (GeV/c)");
1687 fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
1688 outputContainer->Add(fhConvPtMCAntiNeutron) ;
1689
1690 fhConvDispersionMCAntiNeutron = new TH2F
1691 ("hConvDispersionMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.);
1692 fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
1693 fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
1694 outputContainer->Add(fhConvDispersionMCAntiNeutron) ;
1695
1696 fhConvM02MCAntiNeutron = new TH2F
1697 ("hConvM02MCAntiNeutron","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
1698 fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
1699 fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
1700 outputContainer->Add(fhConvM02MCAntiNeutron) ;
1701
1702 fhConvDeltaEtaMCAntiProton = new TH2F
1703 ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5);
1704 fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
1705 fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
1706 outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
1707
1708 fhConvDeltaPhiMCAntiProton = new TH2F
1709 ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5);
1710 fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
1711 fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
1712 outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
1713
1714 fhConvDeltaEtaPhiMCAntiProton = new TH2F
1715 ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
1716 fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
1717 fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
1718 outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;
1719
1720 fhConvAsymMCAntiProton = new TH2F
1721 ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1);
1722 fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
1723 fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
1724 outputContainer->Add(fhConvAsymMCAntiProton) ;
1725
1726 fhConvPtMCAntiProton = new TH2F
1727 ("hConvPtMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,100,0.,10.);
1728 fhConvPtMCAntiProton->SetYTitle("Pair p_{T} (GeV/c)");
1729 fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
1730 outputContainer->Add(fhConvPtMCAntiProton) ;
1731
1732 fhConvDispersionMCAntiProton = new TH2F
1733 ("hConvDispersionMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.);
1734 fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
1735 fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
1736 outputContainer->Add(fhConvDispersionMCAntiProton) ;
1737
1738 fhConvM02MCAntiProton = new TH2F
1739 ("hConvM02MCAntiProton","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
1740 fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
1741 fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
1742 outputContainer->Add(fhConvM02MCAntiProton) ;
1743
1744 fhConvDeltaEtaMCString = new TH2F
1745 ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5);
1746 fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
1747 fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/c^2)");
1748 outputContainer->Add(fhConvDeltaEtaMCString) ;
1749
1750 fhConvDeltaPhiMCString = new TH2F
1751 ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5);
1752 fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
1753 fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/c^2)");
1754 outputContainer->Add(fhConvDeltaPhiMCString) ;
1755
1756 fhConvDeltaEtaPhiMCString = new TH2F
1757 ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5);
1758 fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
1759 fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
1760 outputContainer->Add(fhConvDeltaEtaPhiMCString) ;
1761
1762 fhConvAsymMCString = new TH2F
1763 ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1);
1764 fhConvAsymMCString->SetYTitle("Asymmetry");
1765 fhConvAsymMCString->SetXTitle("Pair Mass (GeV/c^2)");
1766 outputContainer->Add(fhConvAsymMCString) ;
1767
1768 fhConvPtMCString = new TH2F
1769 ("hConvPtMCString","p_{T} of selected conversion pairs from string",100,0,fMassCut,100,0.,10.);
1770 fhConvPtMCString->SetYTitle("Pair p_{T} (GeV/c)");
1771 fhConvPtMCString->SetXTitle("Pair Mass (GeV/c^2)");
1772 outputContainer->Add(fhConvPtMCString) ;
1773
1774 fhConvDispersionMCString = new TH2F
1775 ("hConvDispersionMCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
1776 fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
1777 fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
1778 outputContainer->Add(fhConvDispersionMCString) ;
1779
1780 fhConvM02MCString = new TH2F
1781 ("hConvM02MCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
1782 fhConvM02MCString->SetYTitle("M02 cluster 1");
1783 fhConvM02MCString->SetXTitle("M02 cluster 2");
c4a7d28a 1784 outputContainer->Add(fhConvM02MCString) ;
1785
1786 fhConvDistMCConversion = new TH2F
1787 ("hConvDistMCConversion","calculated conversion distance vs real vertes for MC conversion",100,0.,5.,100,0.,5.);
1788 fhConvDistMCConversion->SetYTitle("distance");
1789 fhConvDistMCConversion->SetXTitle("vertex R");
1790 outputContainer->Add(fhConvDistMCConversion) ;
1791
1792 fhConvDistMCConversionCuts = new TH2F
1793 ("hConvDistMCConversionCuts","calculated conversion distance vs real vertes for MC conversion, deta < 0.05, m < 10 MeV, asym < 0.1",100,0.,5.,100,0.,5.);
1794 fhConvDistMCConversionCuts->SetYTitle("distance");
1795 fhConvDistMCConversionCuts->SetXTitle("vertex R");
1796 outputContainer->Add(fhConvDistMCConversionCuts) ;
20218aea 1797 }
6175da48 1798
521636d2 1799 if(fFillSSHistograms){
1800
3d5d5078 1801 TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
1802
1803 TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
1804
1805 for(Int_t i = 0; i < 6; i++){
521636d2 1806
3d5d5078 1807 fhMCELambda0[i] = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
1808 Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
521636d2 1809 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1810 fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
1811 fhMCELambda0[i]->SetXTitle("E (GeV)");
1812 outputContainer->Add(fhMCELambda0[i]) ;
521636d2 1813
3d5d5078 1814 fhMCEdLambda0[i] = new TH2F(Form("hEdLambda0_MC%s",pnamess[i].Data()),
1815 Form("cluster from %s : E vs #lambda_{0}^{2}/N_{cells}^{2}",ptypess[i].Data()),
521636d2 1816 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
3d5d5078 1817 fhMCEdLambda0[i]->SetYTitle("d#lambda_{0}^{2}");
1818 fhMCEdLambda0[i]->SetXTitle("E (GeV)");
1819 outputContainer->Add(fhMCEdLambda0[i]) ;
521636d2 1820
3d5d5078 1821 fhMCELambda1[i] = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
1822 Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
521636d2 1823 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1824 fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
1825 fhMCELambda1[i]->SetXTitle("E (GeV)");
1826 outputContainer->Add(fhMCELambda1[i]) ;
521636d2 1827
3d5d5078 1828 fhMCEdLambda1[i] = new TH2F(Form("hEdLambda1_MC%s",pnamess[i].Data()),
1829 Form("cluster from %s : E vs d#lambda_{1}^{2}/N_{cells}^{2}",ptypess[i].Data()),
521636d2 1830 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
3d5d5078 1831 fhMCEdLambda1[i]->SetYTitle("d#lambda_{1}^{2}");
1832 fhMCEdLambda1[i]->SetXTitle("E (GeV)");
1833 outputContainer->Add(fhMCEdLambda1[i]) ;
521636d2 1834
3d5d5078 1835 fhMCEDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
1836 Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
521636d2 1837 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3d5d5078 1838 fhMCEDispersion[i]->SetYTitle("D^{2}");
1839 fhMCEDispersion[i]->SetXTitle("E (GeV)");
1840 outputContainer->Add(fhMCEDispersion[i]) ;
521636d2 1841
3d5d5078 1842 fhMCEdDispersion[i] = new TH2F(Form("hEdDispersion_MC%s",pnamess[i].Data()),
1843 Form("cluster from %s : E vs dispersion^{2}/N_{cells}^{2}",ptypess[i].Data()),
521636d2 1844 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
3d5d5078 1845 fhMCEdDispersion[i]->SetYTitle("dD^{2}");
1846 fhMCEdDispersion[i]->SetXTitle("E (GeV)");
1847 outputContainer->Add(fhMCEdDispersion[i]) ;
1848
f66d95af 1849 fhMCNCellsE[i] = new TH2F (Form("hNCellsE_MC%s",pnamess[i].Data()),
1850 Form("# of cells in cluster from %s vs E of clusters",ptypess[i].Data()),
1851 nptbins,ptmin,ptmax, nbins,nmin,nmax);
1852 fhMCNCellsE[i]->SetXTitle("E (GeV)");
1853 fhMCNCellsE[i]->SetYTitle("# of cells in cluster");
1854 outputContainer->Add(fhMCNCellsE[i]);
1855
1856 fhMCMaxCellDiffClusterE[i] = new TH2F (Form("hMaxCellDiffClusterE_MC%s",pnamess[i].Data()),
1857 Form("energy vs difference of cluster energy from %s - max cell energy / cluster energy, good clusters",ptypess[i].Data()),
1858 nptbins,ptmin,ptmax, 500,0,1.);
1859 fhMCMaxCellDiffClusterE[i]->SetXTitle("E_{cluster} (GeV) ");
1860 fhMCMaxCellDiffClusterE[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1861 outputContainer->Add(fhMCMaxCellDiffClusterE[i]);
1862
1863 fhMCLambda0vsClusterMaxCellDiffE0[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
1864 Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
1865 ssbins,ssmin,ssmax,500,0,1.);
1866 fhMCLambda0vsClusterMaxCellDiffE0[i]->SetXTitle("#lambda_{0}^{2}");
1867 fhMCLambda0vsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1868 outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE0[i]) ;
1869
1870 fhMCLambda0vsClusterMaxCellDiffE2[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
1871 Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
1872 ssbins,ssmin,ssmax,500,0,1.);
1873 fhMCLambda0vsClusterMaxCellDiffE2[i]->SetXTitle("#lambda_{0}^{2}");
1874 fhMCLambda0vsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1875 outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE2[i]) ;
1876
1877 fhMCLambda0vsClusterMaxCellDiffE6[i] = new TH2F(Form("hLambda0vsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
1878 Form("cluster from %s : #lambda^{2}_{0} vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
1879 ssbins,ssmin,ssmax,500,0,1.);
1880 fhMCLambda0vsClusterMaxCellDiffE6[i]->SetXTitle("#lambda_{0}^{2}");
1881 fhMCLambda0vsClusterMaxCellDiffE6[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1882 outputContainer->Add(fhMCLambda0vsClusterMaxCellDiffE6[i]) ;
1883
1884 fhMCNCellsvsClusterMaxCellDiffE0[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE0_MC%s",pnamess[i].Data()),
1885 Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E < 2 GeV",ptypess[i].Data()),
1886 nbins/5,nmin,nmax/5,500,0,1.);
1887 fhMCNCellsvsClusterMaxCellDiffE0[i]->SetXTitle("N cells in cluster");
1888 fhMCNCellsvsClusterMaxCellDiffE0[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1889 outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE0[i]) ;
1890
1891 fhMCNCellsvsClusterMaxCellDiffE2[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE2_MC%s",pnamess[i].Data()),
1892 Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, 2< E < 6 GeV",ptypess[i].Data()),
1893 nbins/5,nmin,nmax/5,500,0,1.);
1894 fhMCNCellsvsClusterMaxCellDiffE2[i]->SetXTitle("N cells in cluster");
1895 fhMCNCellsvsClusterMaxCellDiffE2[i]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1896 outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE2[i]) ;
1897
1898 fhMCNCellsvsClusterMaxCellDiffE6[i] = new TH2F(Form("hNCellsvsClusterMaxCellDiffE6_MC%s",pnamess[i].Data()),
1899 Form("cluster from %s : N cells in cluster vs fraction of energy carried by max cell, E > 6 GeV",ptypess[i].Data()),
1900 nbins/5,nmin,nmax/5,500,0,1.);
1901 fhMCNCellsvsClusterMaxCellDiffE6[i]->SetXTitle("N cells in cluster");
1902 fhMCNCellsvsClusterMaxCellDiffE6[i]->SetYTitle("E (GeV)");
1903 outputContainer->Add(fhMCNCellsvsClusterMaxCellDiffE6[i]) ;
1904
3d5d5078 1905 }// loop
1906
1907 if(!GetReader()->IsEmbeddedClusterSelectionOn())
1908 {
1909 fhMCPhotonELambda0NoOverlap = new TH2F("hELambda0_MCPhoton_NoOverlap",
1910 "cluster from Photon : E vs #lambda_{0}^{2}",
1911 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1912 fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
1913 fhMCPhotonELambda0NoOverlap->SetXTitle("E (GeV)");
1914 outputContainer->Add(fhMCPhotonELambda0NoOverlap) ;
1915
1916 fhMCPhotonEdLambda0NoOverlap = new TH2F("hEdLambda0_MCPhoton_NoOverlap",
1917 "cluster from Photon : E vs #lambda_{0}^{2}/N_{cells}^{2}",
1918 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1919 fhMCPhotonEdLambda0NoOverlap->SetYTitle("d#lambda_{0}^{2}");
1920 fhMCPhotonEdLambda0NoOverlap->SetXTitle("E (GeV)");
1921 outputContainer->Add(fhMCPhotonEdLambda0NoOverlap) ;
1922
1923 fhMCPhotonELambda0TwoOverlap = new TH2F("hELambda0_MCPhoton_TwoOverlap",
1924 "cluster from Photon : E vs #lambda_{0}^{2}",
1925 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1926 fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
1927 fhMCPhotonELambda0TwoOverlap->SetXTitle("E (GeV)");
1928 outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ;
1929
1930 fhMCPhotonEdLambda0TwoOverlap = new TH2F("hEdLambda0_MCPhoton_TwoOverlap",
1931 "cluster from Photon : E vs #lambda_{0}^{2}/N_{cells}^{2}",
1932 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1933 fhMCPhotonEdLambda0TwoOverlap->SetYTitle("d#lambda_{0}^{2}");
1934 fhMCPhotonEdLambda0TwoOverlap->SetXTitle("E (GeV)");
1935 outputContainer->Add(fhMCPhotonEdLambda0TwoOverlap) ;
1936
1937 fhMCPhotonELambda0NOverlap = new TH2F("hELambda0_MCPhoton_NOverlap",
1938 "cluster from Photon : E vs #lambda_{0}^{2}",
1939 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1940 fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
1941 fhMCPhotonELambda0NOverlap->SetXTitle("E (GeV)");
1942 outputContainer->Add(fhMCPhotonELambda0NOverlap) ;
521636d2 1943
3d5d5078 1944 fhMCPhotonEdLambda0NOverlap = new TH2F("hEdLambda0_MCPhoton_NOverlap",
1945 "cluster from Photon : E vs #lambda_{0}^{2}/N_{cells}^{2}",
1946 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1947 fhMCPhotonEdLambda0NOverlap->SetYTitle("d#lambda_{0}^{2}");
1948 fhMCPhotonEdLambda0NOverlap->SetXTitle("E (GeV)");
1949 outputContainer->Add(fhMCPhotonEdLambda0NOverlap) ;
1950 } //No embedding
1951
1952 //Fill histograms to check shape of embedded clusters
1953 if(GetReader()->IsEmbeddedClusterSelectionOn())
1954 {
1955
1956 fhEmbeddedSignalFractionEnergy = new TH2F("hEmbeddedSignal_FractionEnergy",
1957 "Energy Fraction of embedded signal versus cluster energy",
1958 nptbins,ptmin,ptmax,100,0.,1.);
1959 fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1960 fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1961 outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
1962
1963 fhEmbedPhotonELambda0FullSignal = new TH2F("hELambda0_EmbedPhoton_FullSignal",
1964 "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1965 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1966 fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1967 fhEmbedPhotonELambda0FullSignal->SetXTitle("E (GeV)");
1968 outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ;
1969
1970 fhEmbedPhotonEdLambda0FullSignal = new TH2F("hEdLambda0_EmbedPhoton_FullSignal",
1971 "cluster from Photon with more than 90% energy in cluster: E vs d#lambda_{0}^{2}",
1972 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1973 fhEmbedPhotonEdLambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1974 fhEmbedPhotonEdLambda0FullSignal->SetXTitle("E (GeV)");
1975 outputContainer->Add(fhEmbedPhotonEdLambda0FullSignal) ;
1976
1977
1978 fhEmbedPhotonELambda0MostlySignal = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
1979 "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1980 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1981 fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1982 fhEmbedPhotonELambda0MostlySignal->SetXTitle("E (GeV)");
1983 outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ;
1984
1985 fhEmbedPhotonEdLambda0MostlySignal = new TH2F("hEdLambda0_EmbedPhoton_MostlySignal",
1986 "cluster from Photon with 50% to 90% energy in cluster: E vs d#lambda_{0}^{2}",
1987 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1988 fhEmbedPhotonEdLambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1989 fhEmbedPhotonEdLambda0MostlySignal->SetXTitle("E (GeV)");
1990 outputContainer->Add(fhEmbedPhotonEdLambda0MostlySignal) ;
1991
1992
1993 fhEmbedPhotonELambda0MostlyBkg = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
1994 "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1995 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1996 fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1997 fhEmbedPhotonELambda0MostlyBkg->SetXTitle("E (GeV)");
1998 outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ;
1999
2000 fhEmbedPhotonEdLambda0MostlyBkg = new TH2F("hEdLambda0_EmbedPhoton_MostlyBkg",
2001 "cluster from Photon with 10% to 50% energy in cluster: E vs d#lambda_{0}^{2}",
2002 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2003 fhEmbedPhotonEdLambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
2004 fhEmbedPhotonEdLambda0MostlyBkg->SetXTitle("E (GeV)");
2005 outputContainer->Add(fhEmbedPhotonEdLambda0MostlyBkg) ;
2006
2007
2008 fhEmbedPhotonELambda0FullBkg = new TH2F("hELambda0_EmbedPhoton_FullBkg",
2009 "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
2010 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2011 fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
2012 fhEmbedPhotonELambda0FullBkg->SetXTitle("E (GeV)");
2013 outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ;
2014
2015 fhEmbedPhotonEdLambda0FullBkg = new TH2F("hEdLambda0_EmbedPhoton_FullBkg",
2016 "cluster from Photon with 0% to 10% energy in cluster: E vs d#lambda_{0}^{2}",
2017 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2018 fhEmbedPhotonEdLambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
2019 fhEmbedPhotonEdLambda0FullBkg->SetXTitle("E (GeV)");
2020 outputContainer->Add(fhEmbedPhotonEdLambda0FullBkg) ;
2021
2022
2023 fhEmbedPi0ELambda0FullSignal = new TH2F("hELambda0_EmbedPi0_FullSignal",
2024 "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
2025 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2026 fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
2027 fhEmbedPi0ELambda0FullSignal->SetXTitle("E (GeV)");
2028 outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ;
2029
2030 fhEmbedPi0EdLambda0FullSignal = new TH2F("hEdLambda0_EmbedPi0_FullSignal",
2031 "cluster from Pi0 with more than 90% energy in cluster: E vs d#lambda_{0}^{2}",
2032 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2033 fhEmbedPi0EdLambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
2034 fhEmbedPi0EdLambda0FullSignal->SetXTitle("E (GeV)");
2035 outputContainer->Add(fhEmbedPi0EdLambda0FullSignal) ;
2036
2037
2038 fhEmbedPi0ELambda0MostlySignal = new TH2F("hELambda0_EmbedPi0_MostlySignal",
2039 "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
2040 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2041 fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
2042 fhEmbedPi0ELambda0MostlySignal->SetXTitle("E (GeV)");
2043 outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ;
2044
2045 fhEmbedPi0EdLambda0MostlySignal = new TH2F("hEdLambda0_EmbedPi0_MostlySignal",
2046 "cluster from Pi0 with 50% to 90% energy in cluster: E vs d#lambda_{0}^{2}",
2047 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2048 fhEmbedPi0EdLambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
2049 fhEmbedPi0EdLambda0MostlySignal->SetXTitle("E (GeV)");
2050 outputContainer->Add(fhEmbedPi0EdLambda0MostlySignal) ;
2051
2052
2053 fhEmbedPi0ELambda0MostlyBkg = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
2054 "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
2055 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2056 fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
2057 fhEmbedPi0ELambda0MostlyBkg->SetXTitle("E (GeV)");
2058 outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ;
2059
2060 fhEmbedPi0EdLambda0MostlyBkg = new TH2F("hEdLambda0_EmbedPi0_MostlyBkg",
2061 "cluster from Pi0 with 10% to 50% energy in cluster: E vs d#lambda_{0}^{2}",
2062 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2063 fhEmbedPi0EdLambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
2064 fhEmbedPi0EdLambda0MostlyBkg->SetXTitle("E (GeV)");
2065 outputContainer->Add(fhEmbedPi0EdLambda0MostlyBkg) ;
2066
2067
2068 fhEmbedPi0ELambda0FullBkg = new TH2F("hELambda0_EmbedPi0_FullBkg",
2069 "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
2070 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2071 fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
2072 fhEmbedPi0ELambda0FullBkg->SetXTitle("E (GeV)");
2073 outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ;
2074
2075 fhEmbedPi0EdLambda0FullBkg = new TH2F("hEdLambda0_EmbedPi0_FullBkg",
2076 "cluster from Pi0 with 0% to 10% energy in cluster: E vs d#lambda_{0}^{2}",
2077 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2078 fhEmbedPi0EdLambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
2079 fhEmbedPi0EdLambda0FullBkg->SetXTitle("E (GeV)");
2080 outputContainer->Add(fhEmbedPi0EdLambda0FullBkg) ;
2081
2082 }// embedded histograms
2083
521636d2 2084
2085 }// Fill SS MC histograms
2086
477d6cee 2087 }//Histos with MC
0c1383b5 2088
d39cba7e 2089 //Store calo PID histograms
2090 if(fRejectTrackMatch){
2091 TList * caloPIDHistos = GetCaloPID()->GetCreateOutputObjects() ;
2092 for(Int_t i = 0; i < caloPIDHistos->GetEntries(); i++) {
2093 outputContainer->Add(caloPIDHistos->At(i)) ;
2094 }
2095 delete caloPIDHistos;
2096 }
2097
477d6cee 2098 return outputContainer ;
2099
1c5acb87 2100}
2101
6639984f 2102//____________________________________________________________________________
2103void AliAnaPhoton::Init()
2104{
2105
2106 //Init
2107 //Do some checks
1e86c71e 2108 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
591cc579 2109 printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 2110 abort();
2111 }
1e86c71e 2112 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
591cc579 2113 printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 2114 abort();
2115 }
2116
2117}
2118
1c5acb87 2119//____________________________________________________________________________
2120void AliAnaPhoton::InitParameters()
2121{
2122
2123 //Initialize the parameters of the analysis.
a3aebfff 2124 AddToHistogramsName("AnaPhoton_");
521636d2 2125
6175da48 2126 fCalorimeter = "EMCAL" ;
2127 fMinDist = 2.;
2128 fMinDist2 = 4.;
2129 fMinDist3 = 5.;
2130 fMassCut = 0.03; //30 MeV
1e86c71e 2131
4cf55759 2132 fTimeCutMin = -1;
2133 fTimeCutMax = 9999999;
6175da48 2134 fNCellsCut = 0;
2ac125bf 2135
1e86c71e 2136 fRejectTrackMatch = kTRUE ;
2137 fCheckConversion = kFALSE;
20218aea 2138 fRemoveConvertedPair = kFALSE;
1e86c71e 2139 fAddConvertedPairsToAOD = kFALSE;
2140
1c5acb87 2141}
2142
2143//__________________________________________________________________
2144void AliAnaPhoton::MakeAnalysisFillAOD()
2145{
f8006433 2146 //Do photon analysis and fill aods
f37fa8d2 2147
6175da48 2148 //Get the vertex
5025c139 2149 Double_t v[3] = {0,0,0}; //vertex ;
2150 GetReader()->GetVertex(v);
f8006433 2151
f37fa8d2 2152 //Select the Calorimeter of the photon
c8fe2783 2153 TObjArray * pl = 0x0;
477d6cee 2154 if(fCalorimeter == "PHOS")
be518ab0 2155 pl = GetPHOSClusters();
477d6cee 2156 else if (fCalorimeter == "EMCAL")
be518ab0 2157 pl = GetEMCALClusters();
5ae09196 2158
2159 if(!pl) {
2160 Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2161 return;
2162 }
521636d2 2163
6175da48 2164 //Init arrays, variables, get number of clusters
1e86c71e 2165 TLorentzVector mom, mom2 ;
2166 Int_t nCaloClusters = pl->GetEntriesFast();
6175da48 2167 //List to be used in conversion analysis, to tag the cluster as candidate for conversion
20218aea 2168 Bool_t * indexConverted = 0x0;
2169 if(fCheckConversion){
2170 indexConverted = new Bool_t[nCaloClusters];
2171 for (Int_t i = 0; i < nCaloClusters; i++)
2172 indexConverted[i] = kFALSE;
2173 }
2174
6175da48 2175 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
521636d2 2176
6175da48 2177 //----------------------------------------------------
2178 // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
2179 //----------------------------------------------------
2180 // Loop on clusters
1e86c71e 2181 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
2182
0ae57829 2183 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
2184 //printf("calo %d, %f\n",icalo,calo->E());
521636d2 2185
f8006433 2186 //Get the index where the cluster comes, to retrieve the corresponding vertex
c8fe2783 2187 Int_t evtIndex = 0 ;
2188 if (GetMixedEvent()) {
2189 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
5025c139 2190 //Get the vertex and check it is not too large in z
96539743 2191 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
c8fe2783 2192 }
521636d2 2193
2194 //Cluster selection, not charged, with photon id and in fiducial cut
f8006433 2195 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
2196 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
2197 else{
2198 Double_t vertex[]={0,0,0};
2199 calo->GetMomentum(mom,vertex) ;
2200 }
c8fe2783 2201
6175da48 2202 //--------------------------------------
2203 // Cluster selection
2204 //--------------------------------------
c4a7d28a 2205 if(!ClusterSelected(calo,mom)) continue;
6175da48 2206
2207 //----------------------------
2208 //Create AOD for analysis
2209 //----------------------------
2210 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
2211
2212 //...............................................
2213 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
2214 Int_t label = calo->GetLabel();
2215 aodph.SetLabel(label);
6175da48 2216 aodph.SetCaloLabel(calo->GetID(),-1);
2217 aodph.SetDetector(fCalorimeter);
c4a7d28a 2218 //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
521636d2 2219
6175da48 2220 //...............................................
2221 //Set bad channel distance bit
c4a7d28a 2222 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
f37fa8d2 2223 if (distBad > fMinDist3) aodph.SetDistToBad(2) ;
477d6cee 2224 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
f37fa8d2 2225 else aodph.SetDistToBad(0) ;
af7b3903 2226 //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
c8fe2783 2227
521636d2 2228 //--------------------------------------------------------------------------------------
2229 //Play with the MC stack if available
2230 //--------------------------------------------------------------------------------------
2231
2232 //Check origin of the candidates
2233 if(IsDataMC()){
2234 aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
3d5d5078 2235
521636d2 2236 if(GetDebug() > 0)
2237 printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
2238 }//Work with stack also
2239
2240 //--------------------------------------------------------------------------------------
2241 //Fill some shower shape histograms before PID is applied
2242 //--------------------------------------------------------------------------------------
2243
2244 FillShowerShapeHistograms(calo,aodph.GetTag());
6175da48 2245
2246 //-------------------------------------
f37fa8d2 2247 //PID selection or bit setting
6175da48 2248 //-------------------------------------
2249 // MC
477d6cee 2250 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
f37fa8d2 2251 //Get most probable PID, check PID weights (in MC this option is mandatory)
21a4b1c0 2252 aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
2253 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
f37fa8d2 2254 //If primary is not photon, skip it.
21a4b1c0 2255 if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
6175da48 2256 }
2257 //...............................................
2258 // Data, PID check on
477d6cee 2259 else if(IsCaloPIDOn()){
f37fa8d2 2260 //Get most probable PID, 2 options check PID weights
3d5d5078 2261 //or redo PID, recommended option for MCEal.
477d6cee 2262 if(!IsCaloPIDRecalculationOn())
21a4b1c0 2263 aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
477d6cee 2264 else
21a4b1c0 2265 aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
477d6cee 2266
21a4b1c0 2267 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
477d6cee 2268
f37fa8d2 2269 //If cluster does not pass pid, not photon, skip it.
21a4b1c0 2270 if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
477d6cee 2271
2272 }
6175da48 2273 //...............................................
2274 // Data, PID check off
477d6cee 2275 else{
f37fa8d2 2276 //Set PID bits for later selection (AliAnaPi0 for example)
2277 //GetPDG already called in SetPIDBits.
f2ccb5b8 2278 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
a3aebfff 2279 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
477d6cee 2280 }
2281
21a4b1c0 2282 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetIdentifiedParticleType());
477d6cee 2283
477d6cee 2284
6175da48 2285 //--------------------------------------------------------------------------------------
2286 // Conversions pairs analysis
f37fa8d2 2287 // Check if cluster comes from a conversion in the material in front of the calorimeter
2288 // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
6175da48 2289 //--------------------------------------------------------------------------------------
521636d2 2290
6175da48 2291 // Do analysis only if there are more than one cluster
20218aea 2292 if( nCaloClusters > 1 && fCheckConversion){
c8fe2783 2293 Bool_t bConverted = kFALSE;
2294 Int_t id2 = -1;
1e86c71e 2295
f37fa8d2 2296 //Check if set previously as converted couple, if so skip its use.
20218aea 2297 if (indexConverted[icalo]) continue;
1e86c71e 2298
6175da48 2299 // Second cluster loop
c8fe2783 2300 for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
f37fa8d2 2301 //Check if set previously as converted couple, if so skip its use.
c8fe2783 2302 if (indexConverted[jcalo]) continue;
f37fa8d2 2303 //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
521636d2 2304 AliVCluster * calo2 = (AliVCluster*) (pl->At(jcalo)); //Get cluster kinematics
6175da48 2305
c4a7d28a 2306
6175da48 2307 //Mixed event, get index of event
c8fe2783 2308 Int_t evtIndex2 = 0 ;
2309 if (GetMixedEvent()) {
2310 evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
f8006433 2311
6175da48 2312 }
2313
2314 //Get kinematics of second cluster
f8006433 2315 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
8cdc266d 2316 calo2->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
f8006433 2317 else{
2318 Double_t vertex[]={0,0,0};
8cdc266d 2319 calo2->GetMomentum(mom2,vertex) ;
f8006433 2320 }
2321
c4a7d28a 2322 //--------------------------------------
2323 // Cluster selection
2324 //--------------------------------------
2325
2326 if(!ClusterSelected(calo2,mom2)) continue;
2327
2328 //................................................
2329 // Get TOF of each cluster in pair, calculate difference if small,
2330 // take this pair. Way to reject clusters from hadrons (or pileup?)
2331 Double_t t12diff = calo2->GetTOF()-calo->GetTOF()*1e9;
2332 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
c8fe2783 2333
6175da48 2334 //................................................
f37fa8d2 2335 //Get mass of pair, if small, take this pair.
41121cfe 2336 Float_t pairM = (mom+mom2).M();
2337 //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
2338 if(pairM < fMassCut){
2339 aodph.SetTagged(kFALSE);
c8fe2783 2340 id2 = calo2->GetID();
6175da48 2341 indexConverted[icalo]=kTRUE;
c4a7d28a 2342 indexConverted[jcalo]=kTRUE;
41121cfe 2343 Float_t asymmetry = TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E());
2344 Float_t dPhi = mom.Phi()-mom2.Phi();
d39cba7e 2345 Float_t dEta = mom.Eta()-mom2.Eta();
41121cfe 2346
c4a7d28a 2347 //...............................................
2348 //Fill few histograms with kinematics of the pair
2349 //FIXME, move all this to MakeAnalysisFillHistograms ...
2350
2351 fhConvDeltaEta ->Fill( pairM, dPhi );
2352 fhConvDeltaPhi ->Fill( pairM, dEta );
2353 fhConvAsym ->Fill( pairM, asymmetry );
2354 fhConvDeltaEtaPhi->Fill( dEta , dPhi );
2355 fhConvPt ->Fill( pairM, (mom+mom2).Pt());
2356
d39cba7e 2357 //Estimate conversion distance, T. Awes, M. Ivanov
2358 //Under the assumption that the pair has zero mass, and that each electron
2359 //of the pair has the same momentum, they will each have the same bend radius
2360 //given by R=p/(qB) = p / (300 B) with p in [MeV/c], B in [Tesla] and R in [m].
2361 //With nominal ALICE magnet current of 30kA B=0.5T, and so with E_cluster=p,
2362 //R = E/1.5 [cm]. Under these assumptions, the distance from the conversion
3d5d5078 2363 //point to the MCEal can be related to the separation distance, L=2y, on the MCEal
d39cba7e 2364 //as d = sqrt(R^2 -(R-y)^2) = sqrt(2Ry - y^2). And since R>>y we can write as
2365 //d = sqrt(E*L/1.5) where E is the cluster energy and L is the distance in cm between
2366 //the clusters.
c4a7d28a 2367 Float_t pos1[3];
2368 calo->GetPosition(pos1);
2369 Float_t pos2[3];
2370 calo2->GetPosition(pos2);
2371 Float_t clustDist = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0])+
521636d2 2372 (pos1[1]-pos2[1])*(pos1[1]-pos2[1])+
2373 (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
c4a7d28a 2374
2375 Float_t convDist = TMath::Sqrt(mom.E() *clustDist*0.01/0.15);
2376 Float_t convDist2 = TMath::Sqrt(mom2.E()*clustDist*0.01/0.15);
2377 //printf("l = %f, e1 = %f, d1=%f, e2 = %f, d2=%f\n",clustDist,mom.E(),convDist,mom2.E(),convDist2);
8cdc266d 2378 if(GetDebug() > 2)
41121cfe 2379 printf("AliAnaPhoton::MakeAnalysisFillAOD(): Pair with mass %2.3f < %2.3f, %1.2f < dPhi %2.2f < %2.2f, dEta %f < %2.2f, asymmetry %2.2f< %2.2f; \n cluster1 id %d, e %2.3f SM %d, eta %2.3f, phi %2.3f ; \n cluster2 id %d, e %2.3f, SM %d,eta %2.3f, phi %2.3f\n",
2380 pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
8cdc266d 2381 calo->GetID(),calo->E(),GetCaloUtils()->GetModuleNumber(calo), mom.Eta(), mom.Phi(),
2382 id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2),mom2.Eta(), mom2.Phi());
6175da48 2383
c4a7d28a 2384 fhConvDistEta ->Fill(mom .Eta(), convDist );
2385 fhConvDistEta ->Fill(mom2.Eta(), convDist2);
521636d2 2386 fhConvDistEn ->Fill(mom .E(), convDist );
2387 fhConvDistEn ->Fill(mom2.E(), convDist2);
c4a7d28a 2388 fhConvDistMass->Fill((mom+mom2).M(), convDist );
2389 //dEta cut
2390 if(dEta<0.05){
2391 fhConvDistEtaCutEta ->Fill(mom .Eta(), convDist );
2392 fhConvDistEtaCutEta ->Fill(mom2.Eta(), convDist2);
2393 fhConvDistEnCutEta ->Fill(mom .E(), convDist );
2394 fhConvDistEnCutEta ->Fill(mom2.E(), convDist2);
2395 fhConvDistMassCutEta->Fill((mom+mom2).M(), convDist );
2396 //mass cut
2397 if(pairM<0.01){//10 MeV
2398 fhConvDistEtaCutMass ->Fill(mom .Eta(), convDist );
2399 fhConvDistEtaCutMass ->Fill(mom2.Eta(), convDist2);
2400 fhConvDistEnCutMass ->Fill(mom .E(), convDist );
2401 fhConvDistEnCutMass ->Fill(mom2.E(), convDist2);
2402 // asymmetry cut
2403 if(asymmetry<0.1){
2404 fhConvDistEtaCutAsy ->Fill(mom .Eta(), convDist );
2405 fhConvDistEtaCutAsy ->Fill(mom2.Eta(), convDist2);
2406 fhConvDistEnCutAsy ->Fill(mom .E(), convDist );
2407 fhConvDistEnCutAsy ->Fill(mom2.E(), convDist2);
2408 }//asymmetry cut
2409 }//mass cut
2410 }//dEta cut
6175da48 2411
2412 //...............................................
2413 //Select pairs in a eta-phi window
41121cfe 2414 if(TMath::Abs(dEta) < fConvDEtaCut &&
2415 TMath::Abs(dPhi) < fConvDPhiMaxCut &&
2416 TMath::Abs(dPhi) > fConvDPhiMinCut &&
2417 asymmetry < fConvAsymCut ){
521636d2 2418 bConverted = kTRUE;
41121cfe 2419 }
2420 //printf("Accepted? %d\n",bConverted);
6175da48 2421 //...........................................
2422 //Fill more histograms, simulated data
2423 //FIXME, move all this to MakeAnalysisFillHistograms ...
2424 if(IsDataMC()){
2425
2426 //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
2427 Int_t ancPDG = 0;
2428 Int_t ancStatus = 0;
2429 TLorentzVector momentum;
c4a7d28a 2430 TVector3 prodVertex;
6175da48 2431 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(calo->GetLabel(), calo2->GetLabel(),
c4a7d28a 2432 GetReader(), ancPDG, ancStatus, momentum, prodVertex);
6175da48 2433
2434 // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
2435 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2436
2437 Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(calo2->GetLabels(),calo2->GetNLabels(),GetReader(), 0);
2438 if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCConversion)){
2439 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
41121cfe 2440 fhConvDeltaEtaMCConversion ->Fill( pairM, dEta );
2441 fhConvDeltaPhiMCConversion ->Fill( pairM, dPhi );
2442 fhConvAsymMCConversion ->Fill( pairM, asymmetry );
2443 fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi );
2444 fhConvPtMCConversion ->Fill( pairM, (mom+mom2).Pt());
6175da48 2445 fhConvDispersionMCConversion ->Fill( calo->GetDispersion(), calo2->GetDispersion());
2446 fhConvM02MCConversion ->Fill( calo->GetM02(), calo2->GetM02());
c4a7d28a 2447 fhConvDistMCConversion ->Fill( convDist , prodVertex.Mag() );
2448 fhConvDistMCConversion ->Fill( convDist2, prodVertex.Mag() );
2449
2450 if(dEta<0.05 && pairM<0.01 && asymmetry<0.1){
2451 fhConvDistMCConversionCuts->Fill( convDist , prodVertex.Mag() );
2452 fhConvDistMCConversionCuts->Fill( convDist2, prodVertex.Mag() );
2453 }
521636d2 2454
6175da48 2455 }
2456 }
2457 else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiNeutron)){
2458 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
41121cfe 2459 fhConvDeltaEtaMCAntiNeutron ->Fill( pairM, dEta );
2460 fhConvDeltaPhiMCAntiNeutron ->Fill( pairM, dPhi );
2461 fhConvAsymMCAntiNeutron ->Fill( pairM, asymmetry );
2462 fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi );
2463 fhConvPtMCAntiNeutron ->Fill( pairM, (mom+mom2).Pt());
6175da48 2464 fhConvDispersionMCAntiNeutron ->Fill( calo->GetDispersion(), calo2->GetDispersion());
2465 fhConvM02MCAntiNeutron ->Fill( calo->GetM02(), calo2->GetM02());
2466 }
2467 }
2468 else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiProton)){
2469 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
41121cfe 2470 fhConvDeltaEtaMCAntiProton ->Fill( pairM, dEta );
2471 fhConvDeltaPhiMCAntiProton ->Fill( pairM, dPhi );
2472 fhConvAsymMCAntiProton ->Fill( pairM, asymmetry );
2473 fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi );
2474 fhConvPtMCAntiProton ->Fill( pairM, (mom+mom2).Pt());
6175da48 2475 fhConvDispersionMCAntiProton ->Fill( calo->GetDispersion(), calo2->GetDispersion());
2476 fhConvM02MCAntiProton ->Fill( calo->GetM02(), calo2->GetM02());
2477 }
2478 }
2479
2480 //Pairs coming from fragmenting pairs.
2481 if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
41121cfe 2482 fhConvDeltaEtaMCString ->Fill( pairM, dPhi);
2483 fhConvDeltaPhiMCString ->Fill( pairM, dPhi);
2484 fhConvAsymMCString ->Fill( pairM, TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
2485 fhConvDeltaEtaPhiMCString ->Fill( dPhi, dPhi );
2486 fhConvPtMCString ->Fill( pairM, (mom+mom2).Pt());
6175da48 2487 fhConvDispersionMCString ->Fill( calo->GetDispersion(), calo2->GetDispersion());
2488 fhConvM02MCString ->Fill( calo->GetM02(), calo2->GetM02());
2489 }
2490
2491 }// Data MC
521636d2 2492
c8fe2783 2493 break;
2494 }
1e86c71e 2495
c8fe2783 2496 }//Mass loop
1e86c71e 2497
6175da48 2498 //..........................................................................................................
2499 //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
c8fe2783 2500 if(bConverted){
6175da48 2501 //Add to AOD
c8fe2783 2502 if(fAddConvertedPairsToAOD){
f37fa8d2 2503 //Create AOD of pair analysis
c8fe2783 2504 TLorentzVector mpair = mom+mom2;
2505 AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
2506 aodpair.SetLabel(aodph.GetLabel());
f8006433 2507 //aodpair.SetInputFileIndex(input);
c8fe2783 2508
f37fa8d2 2509 //printf("Index %d, Id %d\n",icalo, calo->GetID());
2510 //Set the indeces of the original caloclusters
c8fe2783 2511 aodpair.SetCaloLabel(calo->GetID(),id2);
2512 aodpair.SetDetector(fCalorimeter);
21a4b1c0 2513 aodpair.SetIdentifiedParticleType(aodph.GetIdentifiedParticleType());
c8fe2783 2514 aodpair.SetTag(aodph.GetTag());
6175da48 2515 aodpair.SetTagged(kTRUE);
f37fa8d2 2516 //Add AOD with pair object to aod branch
c8fe2783 2517 AddAODParticle(aodpair);
f37fa8d2 2518 //printf("\t \t both added pair\n");
c8fe2783 2519 }
2520
f37fa8d2 2521 //Do not add the current calocluster
20218aea 2522 if(fRemoveConvertedPair) continue;
6175da48 2523 else {
2524 //printf("TAGGED\n");
2525 //Tag this cluster as likely conversion
2526 aodph.SetTagged(kTRUE);
2527 }
c8fe2783 2528 }//converted pair
2529 }//check conversion
f37fa8d2 2530 //printf("\t \t added single cluster %d\n",icalo);
1e86c71e 2531
6175da48 2532 //FIXME, this to MakeAnalysisFillHistograms ...
f66d95af 2533 Int_t absID = 0;
2534 Float_t maxCellFraction = 0;
2535 AliVCaloCells* cells = 0;
6175da48 2536
f66d95af 2537 if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
2538 else cells = GetPHOSCells();
2539
2540 absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
2541
2542 fhMaxCellDiffClusterE->Fill(aodph.E(),maxCellFraction);
2543 fhNCellsE ->Fill(aodph.E(),calo->GetNCells());
2544
f37fa8d2 2545 //Add AOD with photon object to aod branch
477d6cee 2546 AddAODParticle(aodph);
2547
2548 }//loop
2549
4a745797 2550 delete [] indexConverted;
2551
f37fa8d2 2552 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
477d6cee 2553
1c5acb87 2554}
2555
2556//__________________________________________________________________
2557void AliAnaPhoton::MakeAnalysisFillHistograms()
2558{
6175da48 2559 //Fill histograms
f8006433 2560
6175da48 2561 //-------------------------------------------------------------------
577d9801 2562 // Access MC information in stack if requested, check that it exists.
521636d2 2563 AliStack * stack = 0x0;
2564 TParticle * primary = 0x0;
2565 TClonesArray * mcparticles = 0x0;
2566 AliAODMCParticle * aodprimary = 0x0;
2567
577d9801 2568 if(IsDataMC()){
521636d2 2569
577d9801 2570 if(GetReader()->ReadStack()){
2571 stack = GetMCStack() ;
2572 if(!stack) {
3d5d5078 2573 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
2574 abort();
577d9801 2575 }
f8006433 2576
577d9801 2577 }
2578 else if(GetReader()->ReadAODMCParticles()){
f8006433 2579
577d9801 2580 //Get the list of MC particles
521636d2 2581 mcparticles = GetReader()->GetAODMCParticles(0);
2582 if(!mcparticles && GetDebug() > 0) {
3d5d5078 2583 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
577d9801 2584 }
577d9801 2585 }
2586 }// is data and MC
521636d2 2587
6175da48 2588
2589 // Get vertex
2244659d 2590 Double_t v[3] = {0,0,0}; //vertex ;
2591 GetReader()->GetVertex(v);
6175da48 2592 //fhVertex->Fill(v[0],v[1],v[2]);
2593 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
2594
2595 //----------------------------------
577d9801 2596 //Loop on stored AOD photons
2597 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
577d9801 2598 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
521636d2 2599
577d9801 2600 for(Int_t iaod = 0; iaod < naod ; iaod++){
2601 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2602 Int_t pdg = ph->GetIdentifiedParticleType();
521636d2 2603
577d9801 2604 if(GetDebug() > 3)
2605 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
521636d2 2606
577d9801 2607 //If PID used, fill histos with photons in Calorimeter fCalorimeter
2608 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
2609 if(ph->GetDetector() != fCalorimeter) continue;
521636d2 2610
577d9801 2611 if(GetDebug() > 2)
2612 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
521636d2 2613
6175da48 2614 //................................
577d9801 2615 //Fill photon histograms
2616 Float_t ptcluster = ph->Pt();
2617 Float_t phicluster = ph->Phi();
2618 Float_t etacluster = ph->Eta();
2619 Float_t ecluster = ph->E();
521636d2 2620
20218aea 2621 fhEPhoton ->Fill(ecluster);
577d9801 2622 fhPtPhoton ->Fill(ptcluster);
2623 fhPhiPhoton ->Fill(ptcluster,phicluster);
2624 fhEtaPhoton ->Fill(ptcluster,etacluster);
521636d2 2625 if (ecluster > 0.5) fhEtaPhiPhoton ->Fill(etacluster, phicluster);
20218aea 2626 else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
2627
2628 if(fCheckConversion &&ph->IsTagged()){
2629 fhPtPhotonConv->Fill(ptcluster);
2630 if(ecluster > 0.5) fhEtaPhiPhotonConv ->Fill(etacluster, phicluster);
2631 else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
6175da48 2632 }
20218aea 2633
6175da48 2634 //.......................................
577d9801 2635 //Play with the MC data if available
2636 if(IsDataMC()){
521636d2 2637
3d5d5078 2638 FillAcceptanceHistograms();
2639
577d9801 2640 Int_t tag =ph->GetTag();
521636d2 2641
f66d95af 2642 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[mcPhoton])
3d5d5078 2643 {
2644 fhMCE [mcPhoton] ->Fill(ecluster);
2645 fhPtMC [mcPhoton] ->Fill(ptcluster);
2646 fhPhiMC[mcPhoton] ->Fill(ecluster,phicluster);
2647 fhEtaMC[mcPhoton] ->Fill(ecluster,etacluster);
2648
f66d95af 2649 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[mcConversion])
3d5d5078 2650 {
2651 fhMCE [mcConversion] ->Fill(ecluster);
2652 fhPtMC [mcConversion] ->Fill(ptcluster);
2653 fhPhiMC[mcConversion] ->Fill(ecluster,phicluster);
2654 fhEtaMC[mcConversion] ->Fill(ecluster,etacluster);
2655
2656 if(fCheckConversion){
2657 if(ph->IsTagged()) fhPtConversionTagged ->Fill(ptcluster);
2658 if(ptcluster > 0.5)fhEtaPhiConversion ->Fill(etacluster,phicluster);
2659 else fhEtaPhi05Conversion ->Fill(etacluster,phicluster);
2660 }
2661 }
2662
f66d95af 2663 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt) && fhMCE[mcPrompt]){
3d5d5078 2664 fhMCE [mcPrompt] ->Fill(ecluster);
2665 fhPtMC [mcPrompt] ->Fill(ptcluster);
2666 fhPhiMC[mcPrompt] ->Fill(ecluster,phicluster);
2667 fhEtaMC[mcPrompt] ->Fill(ecluster,etacluster);
2668 }
f66d95af 2669 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation)&& fhMCE[mcFragmentation])
3d5d5078 2670 {
2671 fhMCE [mcFragmentation] ->Fill(ecluster);
2672 fhPtMC [mcFragmentation] ->Fill(ptcluster);
2673 fhPhiMC[mcFragmentation] ->Fill(ecluster,phicluster);
2674 fhEtaMC[mcFragmentation] ->Fill(ecluster,etacluster);
2675
2676 }
f66d95af 2677 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR)&& fhMCE[mcISR])
3d5d5078 2678 {
2679 fhMCE [mcISR] ->Fill(ecluster);
2680 fhPtMC [mcISR] ->Fill(ptcluster);
2681 fhPhiMC[mcISR] ->Fill(ecluster,phicluster);
2682 fhEtaMC[mcISR] ->Fill(ecluster,etacluster);
2683 }
2684 else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
f66d95af 2685 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[mcPi0Decay])
3d5d5078 2686 {
2687 fhMCE [mcPi0Decay] ->Fill(ecluster);
2688 fhPtMC [mcPi0Decay] ->Fill(ptcluster);
2689 fhPhiMC[mcPi0Decay] ->Fill(ecluster,phicluster);
2690 fhEtaMC[mcPi0Decay] ->Fill(ecluster,etacluster);
2691 }
2692 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
f66d95af 2693 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) && fhMCE[mcOtherDecay])
3d5d5078 2694 {
2695 fhMCE [mcOtherDecay] ->Fill(ecluster);
2696 fhPtMC [mcOtherDecay] ->Fill(ptcluster);
2697 fhPhiMC[mcOtherDecay] ->Fill(ecluster,phicluster);
2698 fhEtaMC[mcOtherDecay] ->Fill(ecluster,etacluster);
2699 }
f66d95af 2700 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE [mcPi0])
3d5d5078 2701 {
2702 fhMCE [mcPi0] ->Fill(ecluster);
2703 fhPtMC [mcPi0] ->Fill(ptcluster);
2704 fhPhiMC[mcPi0] ->Fill(ecluster,phicluster);
2705 fhEtaMC[mcPi0] ->Fill(ecluster,etacluster);
f66d95af 2706 }
2707 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[mcEta])
2708 {
2709 fhMCE [mcEta] ->Fill(ecluster);
2710 fhPtMC [mcEta] ->Fill(ptcluster);
2711 fhPhiMC[mcEta] ->Fill(ecluster,phicluster);
2712 fhEtaMC[mcEta] ->Fill(ecluster,etacluster);
2713 }
3d5d5078 2714 }
f66d95af 2715 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[mcAntiNeutron])
3d5d5078 2716 {
2717 fhMCE [mcAntiNeutron] ->Fill(ecluster);
2718 fhPtMC [mcAntiNeutron] ->Fill(ptcluster);
2719 fhPhiMC[mcAntiNeutron] ->Fill(ecluster,phicluster);
2720 fhEtaMC[mcAntiNeutron] ->Fill(ecluster,etacluster);
2721 if(ph->IsTagged() && fCheckConversion) fhPtAntiNeutronTagged ->Fill(ptcluster);
2722
2723 }
f66d95af 2724 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[mcAntiProton])
3d5d5078 2725 {
2726 fhMCE [mcAntiProton] ->Fill(ecluster);
2727 fhPtMC [mcAntiProton] ->Fill(ptcluster);
2728 fhPhiMC[mcAntiProton] ->Fill(ecluster,phicluster);
2729 fhEtaMC[mcAntiProton] ->Fill(ecluster,etacluster);
2730 if(ph->IsTagged() && fCheckConversion) fhPtAntiProtonTagged ->Fill(ptcluster);
2731 }
f66d95af 2732 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[mcElectron])
3d5d5078 2733 {
2734 fhMCE [mcElectron] ->Fill(ecluster);
2735 fhPtMC [mcElectron] ->Fill(ptcluster);
2736 fhPhiMC[mcElectron] ->Fill(ecluster,phicluster);
2737 fhEtaMC[mcElectron] ->Fill(ecluster,etacluster);
2738 }
f66d95af 2739 else if( fhMCE[mcOther]){
3d5d5078 2740 fhMCE [mcOther] ->Fill(ecluster);
521636d2 2741 fhPtMC [mcOther] ->Fill(ptcluster);
2742 fhPhiMC[mcOther] ->Fill(ecluster,phicluster);
2743 fhEtaMC[mcOther] ->Fill(ecluster,etacluster);
20218aea 2744 if(ph->IsTagged() && fCheckConversion) fhPtUnknownTagged ->Fill(ptcluster);
521636d2 2745
3d5d5078 2746
f8006433 2747 // printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
2748 // ph->GetLabel(),ph->Pt());
2749 // for(Int_t i = 0; i < 20; i++) {
2750 // if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
2751 // }
2752 // printf("\n");
2753
577d9801 2754 }
521636d2 2755
577d9801 2756 //....................................................................
2757 // Access MC information in stack if requested, check that it exists.
2758 Int_t label =ph->GetLabel();
2759 if(label < 0) {
3d5d5078 2760 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
2761 continue;
577d9801 2762 }
521636d2 2763
577d9801 2764 Float_t eprim = 0;
2765 Float_t ptprim = 0;
2766 if(GetReader()->ReadStack()){
3d5d5078 2767
2768 if(label >= stack->GetNtrack()) {
f8006433 2769 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
2770 continue ;
3d5d5078 2771 }
2772
2773 primary = stack->Particle(label);
2774 if(!primary){
2775 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
2776 continue;
2777 }
2778 eprim = primary->Energy();
2779 ptprim = primary->Pt();
2780
577d9801 2781 }
2782 else if(GetReader()->ReadAODMCParticles()){
3d5d5078 2783 //Check which is the input
2784 if(ph->GetInputFileIndex() == 0){
2785 if(!mcparticles) continue;
2786 if(label >= mcparticles->GetEntriesFast()) {
2787 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
521636d2 2788 label, mcparticles->GetEntriesFast());
3d5d5078 2789 continue ;
f8006433 2790 }
2791 //Get the particle
521636d2 2792 aodprimary = (AliAODMCParticle*) mcparticles->At(label);
f8006433 2793
3d5d5078 2794 }
2795
2796 if(!aodprimary){
2797 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
2798 continue;
2799 }
2800
2801 eprim = aodprimary->E();
2802 ptprim = aodprimary->Pt();
2803
577d9801 2804 }
521636d2 2805
577d9801 2806 fh2E ->Fill(ecluster, eprim);
2807 fh2Pt ->Fill(ptcluster, ptprim);
2808 fhDeltaE ->Fill(eprim-ecluster);
2809 fhDeltaPt->Fill(ptprim-ptcluster);
2810 if(eprim > 0) fhRatioE ->Fill(ecluster/eprim);
2811 if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
521636d2 2812
577d9801 2813 }//Histograms with MC
521636d2 2814
577d9801 2815 }// aod loop
521636d2 2816
1c5acb87 2817}
2818
2819
2820//__________________________________________________________________
2821void AliAnaPhoton::Print(const Option_t * opt) const
2822{
477d6cee 2823 //Print some relevant parameters set for the analysis
2824
2825 if(! opt)
2826 return;
2827
2828 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2829 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 2830
477d6cee 2831 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
2832 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
2833 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2834 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 2835 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1e86c71e 2836 printf("Check Pair Conversion = %d\n",fCheckConversion);
2837 printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD);
2838 printf("Conversion pair mass cut = %f\n",fMassCut);
41121cfe 2839 printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
2840 fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);
4cf55759 2841 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2ac125bf 2842 printf("Number of cells in cluster is > %d \n", fNCellsCut);
477d6cee 2843 printf(" \n") ;
1c5acb87 2844
2845}