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