]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaPhoton.cxx
Coverity fixed
[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//____________________________________________________________________________
53 AliAnaPhoton::AliAnaPhoton() :
54 AliAnaPartCorrBaseClass(), fCalorimeter(""),
a3aebfff 55 fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
c4a7d28a 56 fTimeCutMin(-1), fTimeCutMax(9999999), fNCellsCut(0), fFillSSHistograms(kFALSE),
20218aea 57 fCheckConversion(kFALSE), fRemoveConvertedPair(kFALSE), fAddConvertedPairsToAOD(kFALSE), fMassCut(0),
41121cfe 58 fConvAsymCut(1.), fConvDEtaCut(2.),fConvDPhiMinCut(-1.), fConvDPhiMaxCut(7.),
c4a7d28a 59 // Histograms
60 fhNtraNclu(0), fhNCellsE(0),
61 //Shower shape histograms
62 fhNCellsLam0LowE(0), fhNCellsLam1LowE(0), fhNCellsDispLowE(0),
63 fhNCellsLam0HighE(0), fhNCellsLam1HighE(0), fhNCellsDispHighE(0),
64 fhEtaLam0(0), fhPhiLam0(0),
65 fhLam1Lam0LowE(0), fhLam1Lam0HighE(0),
66 fhLam0E(0), fhLam1E(0),
67 // Spectra histograms
20218aea 68 fhEPhoton(0), fhPtPhoton(0), fhPhiPhoton(0), fhEtaPhoton(0), fhEtaPhiPhoton(0), fhEtaPhi05Photon(0),
c4a7d28a 69 // Conversion histograms
6175da48 70 fhPtPhotonConv(0), fhEtaPhiPhotonConv(0),fhEtaPhi05PhotonConv(0),
41121cfe 71 fhConvDeltaEta(0), fhConvDeltaPhi(0), fhConvDeltaEtaPhi(0), fhConvAsym(0), fhConvPt(0),
c4a7d28a 72 fhConvDistEta(0), fhConvDistEn(0), fhConvDistMass(0),
73 fhConvDistEtaCutEta(0), fhConvDistEnCutEta(0), fhConvDistMassCutEta(0),
74 fhConvDistEtaCutMass(0), fhConvDistEnCutMass(0), fhConvDistEtaCutAsy(0), fhConvDistEnCutAsy(0),
7175a03a 75 //MC
123fc3bd 76 fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
0ae57829 77 fhPtMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0),
1c5acb87 78 fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
79 fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
6639984f 80 fhPtISR(0),fhPhiISR(0),fhEtaISR(0),
1c5acb87 81 fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
6175da48 82 fhPtOtherDecay(0), fhPhiOtherDecay(0), fhEtaOtherDecay(0),
41121cfe 83 fhPtConversion(0), fhPhiConversion(0), fhEtaConversion(0),fhEtaPhiConversion(0),fhEtaPhi05Conversion(0),
6175da48 84 fhPtAntiNeutron(0), fhPhiAntiNeutron(0), fhEtaAntiNeutron(0),
85 fhPtAntiProton(0), fhPhiAntiProton(0), fhEtaAntiProton(0),
86 fhPtUnknown(0), fhPhiUnknown(0), fhEtaUnknown(0),
41121cfe 87 fhPtConversionTagged(0), fhPtAntiNeutronTagged(0), fhPtAntiProtonTagged(0), fhPtUnknownTagged(0),
6175da48 88 fhConvDeltaEtaMCConversion(0), fhConvDeltaPhiMCConversion(0), fhConvDeltaEtaPhiMCConversion(0), fhConvAsymMCConversion(0), fhConvPtMCConversion(0), fhConvDispersionMCConversion(0), fhConvM02MCConversion(0),
89 fhConvDeltaEtaMCAntiNeutron(0), fhConvDeltaPhiMCAntiNeutron(0), fhConvDeltaEtaPhiMCAntiNeutron(0), fhConvAsymMCAntiNeutron(0), fhConvPtMCAntiNeutron(0), fhConvDispersionMCAntiNeutron(0),fhConvM02MCAntiNeutron(0),
90 fhConvDeltaEtaMCAntiProton(0), fhConvDeltaPhiMCAntiProton(0), fhConvDeltaEtaPhiMCAntiProton(0), fhConvAsymMCAntiProton(0), fhConvPtMCAntiProton(0), fhConvDispersionMCAntiProton(0), fhConvM02MCAntiProton(0),
c4a7d28a 91 fhConvDeltaEtaMCString(0), fhConvDeltaPhiMCString(0), fhConvDeltaEtaPhiMCString(0), fhConvAsymMCString(0), fhConvPtMCString(0), fhConvDispersionMCString(0), fhConvM02MCString(0),
92 fhConvDistMCConversion(0), fhConvDistMCConversionCuts(0)
1c5acb87 93{
94 //default ctor
95
96 //Initialize parameters
97 InitParameters();
98
5ae09196 99}//____________________________________________________________________________
1c5acb87 100AliAnaPhoton::~AliAnaPhoton()
101{
102 //dtor
103
104}
105
c4a7d28a 106//__________________________________________________________________
107Bool_t AliAnaPhoton::ClusterSelected(AliVCluster* calo, TLorentzVector mom)
108{
109 //Select clusters if they pass different cuts
110 if(GetDebug() > 2)
111 printf("AliAnaPhoton::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
112 GetReader()->GetEventNumber(),
113 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
114
115 //.......................................
116 //If too small or big energy, skip it
117 if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ;
118 if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
119
120 //.......................................
121 // TOF cut, BE CAREFUL WITH THIS CUT
122 Double_t tof = calo->GetTOF()*1e9;
123 if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
124 if(GetDebug() > 2) printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
125
126 //.......................................
127 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
128 if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
129
130 //.......................................
131 //Check acceptance selection
132 if(IsFiducialCutOn()){
133 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
134 if(! in ) return kFALSE ;
135 }
136 if(GetDebug() > 2) printf("Fiducial cut passed \n");
137
138 //.......................................
139 //Skip matched clusters with tracks
140 if(fRejectTrackMatch){
141 if(IsTrackMatched(calo)) {
142 if(GetDebug() > 2) printf("\t Reject track-matched clusters\n");
143 return kFALSE ;
144 }
145 else
146 if(GetDebug() > 2) printf(" Track-matching cut passed \n");
147 }// reject matched clusters
148
149 //.......................................
150 //Check Distance to Bad channel, set bit.
151 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
152 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
153 if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
154 return kFALSE ;
155 }
156 else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
157 //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
158
159 if(GetDebug() > 0)
160 printf("AliAnaPhoton::ClusterSelected() Current Event %d; After selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
161 GetReader()->GetEventNumber(),
162 mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
163
164 //All checks passed, cluster selected
165 return kTRUE;
166
167}
168
0c1383b5 169//________________________________________________________________________
170TObjString * AliAnaPhoton::GetAnalysisCuts()
171{
172 //Save parameters used for analysis
173 TString parList ; //this will be list of parameters used for this analysis.
5ae09196 174 const Int_t buffersize = 255;
175 char onePar[buffersize] ;
0c1383b5 176
5ae09196 177 snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
0c1383b5 178 parList+=onePar ;
5ae09196 179 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
0c1383b5 180 parList+=onePar ;
5ae09196 181 snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
0c1383b5 182 parList+=onePar ;
5ae09196 183 snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
0c1383b5 184 parList+=onePar ;
5ae09196 185 snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
0c1383b5 186 parList+=onePar ;
5ae09196 187 snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
0c1383b5 188 parList+=onePar ;
41121cfe 189 snprintf(onePar,buffersize,"Conversion Selection: fConvAsymCut %1.2f, fConvDEtaCut %1.2f fConvDPhiCut (%1.2f,%1.2f)\n",
190 fConvAsymCut, fConvDEtaCut, fConvDPhiMinCut, fConvDPhiMaxCut) ;
6175da48 191 parList+=onePar ;
0c1383b5 192
193 //Get parameters set in base class.
194 parList += GetBaseParametersList() ;
195
196 //Get parameters set in PID class.
197 parList += GetCaloPID()->GetPIDParametersList() ;
198
199 //Get parameters set in FiducialCut class (not available yet)
200 //parlist += GetFidCut()->GetFidCutParametersList()
201
202 return new TObjString(parList) ;
203}
204
1c5acb87 205
206//________________________________________________________________________
207TList * AliAnaPhoton::GetCreateOutputObjects()
208{
477d6cee 209 // Create histograms to be saved in output file and
210 // store them in outputContainer
211 TList * outputContainer = new TList() ;
212 outputContainer->SetName("PhotonHistos") ;
4a745797 213
5a2dbc3c 214 Int_t nptbins = GetHistoPtBins();
215 Int_t nphibins = GetHistoPhiBins();
216 Int_t netabins = GetHistoEtaBins();
477d6cee 217 Float_t ptmax = GetHistoPtMax();
218 Float_t phimax = GetHistoPhiMax();
219 Float_t etamax = GetHistoEtaMax();
220 Float_t ptmin = GetHistoPtMin();
221 Float_t phimin = GetHistoPhiMin();
222 Float_t etamin = GetHistoEtaMin();
223
2244659d 224 fhNtraNclu = new TH2F ("hNtracksNcluster","# of tracks vs # of clusters", 500,0,500, 500,0,500);
225 fhNtraNclu->SetXTitle("# of tracks");
226 fhNtraNclu->SetYTitle("# of clusters");
227 outputContainer->Add(fhNtraNclu);
228
c4a7d28a 229
230 fhNCellsE = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", 100,0, 20, 20,0,20);
231 fhNCellsE->SetXTitle("E (GeV)");
232 fhNCellsE->SetYTitle("# of cells in cluster");
233 outputContainer->Add(fhNCellsE);
234
235 if(fFillSSHistograms){
236
237 fhNCellsLam0LowE = new TH2F ("hNCellsLam0LowE","# of cells in cluster vs #lambda_{0}, E < 2 GeV", 20,0, 20, 200,0,1);
238 fhNCellsLam0LowE->SetXTitle("N Cells");
239 fhNCellsLam0LowE->SetYTitle("#lambda_{0}");
240 outputContainer->Add(fhNCellsLam0LowE);
241
242
243 fhNCellsLam0HighE = new TH2F ("hNCellsLam0HighE","# of cells in cluster vs #lambda_{0}, E > 2 GeV", 20,0, 20, 200,0,1);
244 fhNCellsLam0HighE->SetXTitle("N Cells");
245 fhNCellsLam0HighE->SetYTitle("#lambda_{0}");
246 outputContainer->Add(fhNCellsLam0HighE);
247
248 fhNCellsLam1LowE = new TH2F ("hNCellsLam1LowE","# of cells in cluster vs #lambda_{1}, E < 2 GeV", 20,0, 20, 200,0,1);
249 fhNCellsLam1LowE->SetXTitle("N Cells");
250 fhNCellsLam1LowE->SetYTitle("#lambda_{0}");
251 outputContainer->Add(fhNCellsLam1LowE);
252
253 fhNCellsLam1HighE = new TH2F ("hNCellsLam1HighE","# of cells in cluster vs #lambda_{1}, E > 2 GeV", 20,0, 20, 200,0,1);
254 fhNCellsLam1HighE->SetXTitle("N Cells");
255 fhNCellsLam1HighE->SetYTitle("#lambda_{0}");
256 outputContainer->Add(fhNCellsLam1HighE);
257
258
259 fhNCellsDispLowE = new TH2F ("hNCellsDispLowE","# of cells in cluster vs dispersion, E < 2 GeV", 20,0, 20, 200,0,2);
260 fhNCellsDispLowE->SetXTitle("N Cells");
261 fhNCellsDispLowE->SetYTitle("dispersion");
262 outputContainer->Add(fhNCellsDispLowE);
263
264 fhNCellsDispHighE = new TH2F ("hNCellsDispHighE","# of cells in cluster vs dispersion, E > 2 GeV", 20,0, 20, 200,0,2);
265 fhNCellsDispHighE->SetXTitle("N Cells");
266 fhNCellsDispHighE->SetYTitle("dispersion");
267 outputContainer->Add(fhNCellsDispHighE);
268
269 fhEtaLam0 = new TH2F ("hEtaLam0","#eta vs #lambda_{0}",200,-0.8,0.8, 200,0,1);
270 fhEtaLam0->SetYTitle("#lambda_{0}");
271 fhEtaLam0->SetXTitle("#eta");
272 outputContainer->Add(fhEtaLam0);
273
274 fhPhiLam0 = new TH2F ("hPhiLam0","#phi vs #lambda_{0}", 200,80*TMath::DegToRad(),120*TMath::DegToRad(), 200,0,1);
275 fhPhiLam0->SetYTitle("#lambda_{0}");
276 fhPhiLam0->SetXTitle("#phi");
277 outputContainer->Add(fhPhiLam0);
278
279
280 fhLam1Lam0LowE = new TH2F ("hLam1Lam0LowE","#lambda_{0} vs #lambda_{1} in cluster of E < 2 GeV", 200,0,1, 200,0,1);
281 fhLam1Lam0LowE->SetYTitle("#lambda_{0}");
282 fhLam1Lam0LowE->SetXTitle("#lambda_{1}");
283 outputContainer->Add(fhLam1Lam0LowE);
284
285
286 fhLam1Lam0HighE = new TH2F ("hLam1Lam0HighE","#lambda_{0} vs #lambda_{1} in cluster of E > 2 GeV", 200,0,1, 200,0,1);
287 fhLam1Lam0HighE->SetYTitle("#lambda_{0}");
288 fhLam1Lam0HighE->SetXTitle("#lambda_{1}");
289 outputContainer->Add(fhLam1Lam0HighE);
290
291
292 fhLam0E = new TH2F ("hLam0E","#lambda_{0} vs E", 100,0, 20, 200,0,1);
293 fhLam0E->SetYTitle("#lambda_{0}");
294 fhLam0E->SetXTitle("E (GeV)");
295 outputContainer->Add(fhLam0E);
296
297 fhLam1E = new TH2F ("hLam1E","#lambda_{1} vs E", 100,0, 20, 200,0,1);
298 fhLam1E->SetYTitle("#lambda_{1}");
299 fhLam1E->SetXTitle("E (GeV)");
300 outputContainer->Add(fhLam1E);
301 }
6175da48 302
20218aea 303 fhEPhoton = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
304 fhEPhoton->SetYTitle("N");
305 fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
306 outputContainer->Add(fhEPhoton) ;
307
308 fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax);
477d6cee 309 fhPtPhoton->SetYTitle("N");
310 fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
311 outputContainer->Add(fhPtPhoton) ;
312
313 fhPhiPhoton = new TH2F
20218aea 314 ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
6175da48 315 fhPhiPhoton->SetYTitle("#phi (rad)");
477d6cee 316 fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
317 outputContainer->Add(fhPhiPhoton) ;
318
319 fhEtaPhoton = new TH2F
20218aea 320 ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 321 fhEtaPhoton->SetYTitle("#eta");
322 fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
323 outputContainer->Add(fhEtaPhoton) ;
324
6175da48 325 fhEtaPhiPhoton = new TH2F
326 ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
327 fhEtaPhiPhoton->SetYTitle("#phi (rad)");
328 fhEtaPhiPhoton->SetXTitle("#eta");
329 outputContainer->Add(fhEtaPhiPhoton) ;
20218aea 330 if(GetMinPt() < 0.5){
331 fhEtaPhi05Photon = new TH2F
332 ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
333 fhEtaPhi05Photon->SetYTitle("#phi (rad)");
334 fhEtaPhi05Photon->SetXTitle("#eta");
335 outputContainer->Add(fhEtaPhi05Photon) ;
336 }
6175da48 337
338 //Conversion
20218aea 339 if(fCheckConversion){
340 fhPtPhotonConv = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax);
341 fhPtPhotonConv->SetYTitle("N");
342 fhPtPhotonConv->SetXTitle("p_{T #gamma}(GeV/c)");
343 outputContainer->Add(fhPtPhotonConv) ;
344
345 fhEtaPhiPhotonConv = new TH2F
346 ("hEtaPhiPhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
347 fhEtaPhiPhotonConv->SetYTitle("#phi (rad)");
348 fhEtaPhiPhotonConv->SetXTitle("#eta");
349 outputContainer->Add(fhEtaPhiPhotonConv) ;
350 if(GetMinPt() < 0.5){
351 fhEtaPhi05PhotonConv = new TH2F
352 ("hEtaPhi05PhotonConv","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
353 fhEtaPhi05PhotonConv->SetYTitle("#phi (rad)");
354 fhEtaPhi05PhotonConv->SetXTitle("#eta");
355 outputContainer->Add(fhEtaPhi05PhotonConv) ;
356 }
357
358 fhConvDeltaEta = new TH2F
359 ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins*2,-0.5,0.5);
360 fhConvDeltaEta->SetYTitle("#Delta #eta");
361 fhConvDeltaEta->SetXTitle("Pair Mass (GeV/c^2)");
362 outputContainer->Add(fhConvDeltaEta) ;
363
364 fhConvDeltaPhi = new TH2F
365 ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins*2,-0.5,0.5);
366 fhConvDeltaPhi->SetYTitle("#Delta #phi");
367 fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/c^2)");
368 outputContainer->Add(fhConvDeltaPhi) ;
369
370 fhConvDeltaEtaPhi = new TH2F
371 ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5);
372 fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
373 fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
374 outputContainer->Add(fhConvDeltaEtaPhi) ;
375
376 fhConvAsym = new TH2F
377 ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1);
378 fhConvAsym->SetYTitle("Asymmetry");
379 fhConvAsym->SetXTitle("Pair Mass (GeV/c^2)");
380 outputContainer->Add(fhConvAsym) ;
381
382 fhConvPt = new TH2F
383 ("hConvPt","p_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.);
384 fhConvPt->SetYTitle("Pair p_{T} (GeV/c)");
385 fhConvPt->SetXTitle("Pair Mass (GeV/c^2)");
386 outputContainer->Add(fhConvPt) ;
c4a7d28a 387
388 fhConvDistEta = new TH2F
389 ("hConvDistEta","distance to conversion vertex",100,-0.7,0.7,100,0.,5.);
390 fhConvDistEta->SetXTitle("#eta");
391 fhConvDistEta->SetYTitle(" distance (m)");
392 outputContainer->Add(fhConvDistEta) ;
393
394 fhConvDistEn = new TH2F
395 ("hConvDistEn","distance to conversion vertex",nptbins,ptmin,ptmax,100,0.,5.);
396 fhConvDistEn->SetXTitle("E (GeV)");
397 fhConvDistEn->SetYTitle(" distance (m)");
398 outputContainer->Add(fhConvDistEn) ;
399
400 fhConvDistMass = new TH2F
401 ("hConvDistMass","distance to conversion vertex",100,0,fMassCut,100,0.,5.);
402 fhConvDistMass->SetXTitle("m (GeV/c^2)");
403 fhConvDistMass->SetYTitle(" distance (m)");
404 outputContainer->Add(fhConvDistMass) ;
405
406 fhConvDistEtaCutEta = new TH2F
407 ("hConvDistEtaCutEta","distance to conversion vertex, dEta < 0.05",100,-0.7,0.7,100,0.,5.);
408 fhConvDistEtaCutEta->SetXTitle("#eta");
409 fhConvDistEtaCutEta->SetYTitle(" distance (m)");
410 outputContainer->Add(fhConvDistEtaCutEta) ;
411
412 fhConvDistEnCutEta = new TH2F
413 ("hConvDistEnCutEta","distance to conversion vertex, dEta < 0.05",nptbins,ptmin,ptmax,100,0.,5.);
414 fhConvDistEnCutEta->SetXTitle("E (GeV)");
415 fhConvDistEnCutEta->SetYTitle(" distance (m)");
416 outputContainer->Add(fhConvDistEnCutEta) ;
417
418 fhConvDistMassCutEta = new TH2F
419 ("hConvDistMassCutEta","distance to conversion vertex, dEta < 0.05",100,0,fMassCut,100,0.,5.);
420 fhConvDistMassCutEta->SetXTitle("m (GeV/c^2)");
421 fhConvDistMassCutEta->SetYTitle(" distance (m)");
422 outputContainer->Add(fhConvDistMassCutEta) ;
423
424 fhConvDistEtaCutMass = new TH2F
425 ("hConvDistEtaCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",100,-0.7,0.7,100,0.,5.);
426 fhConvDistEtaCutMass->SetXTitle("#eta");
427 fhConvDistEtaCutMass->SetYTitle(" distance (m)");
428 outputContainer->Add(fhConvDistEtaCutMass) ;
429
430 fhConvDistEnCutMass = new TH2F
431 ("hConvDistEnCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",nptbins,ptmin,ptmax,100,0.,5.);
432 fhConvDistEnCutMass->SetXTitle("E (GeV)");
433 fhConvDistEnCutMass->SetYTitle(" distance (m)");
434 outputContainer->Add(fhConvDistEnCutMass) ;
435
436 fhConvDistEtaCutAsy = new TH2F
437 ("hConvDistEtaCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",100,-0.7,0.7,100,0.,5.);
438 fhConvDistEtaCutAsy->SetXTitle("#eta");
439 fhConvDistEtaCutAsy->SetYTitle(" distance (m)");
440 outputContainer->Add(fhConvDistEtaCutAsy) ;
441
442 fhConvDistEnCutAsy = new TH2F
443 ("hConvDistEnCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",nptbins,ptmin,ptmax,100,0.,5.);
444 fhConvDistEnCutAsy->SetXTitle("E (GeV)");
445 fhConvDistEnCutAsy->SetYTitle(" distance (m)");
446 outputContainer->Add(fhConvDistEnCutAsy) ;
447
20218aea 448 }
6175da48 449
477d6cee 450 if(IsDataMC()){
123fc3bd 451 fhDeltaE = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
452 fhDeltaE->SetXTitle("#Delta E (GeV)");
453 outputContainer->Add(fhDeltaE);
454
455 fhDeltaPt = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
456 fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
457 outputContainer->Add(fhDeltaPt);
458
459 fhRatioE = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
460 fhRatioE->SetXTitle("E_{reco}/E_{gen}");
461 outputContainer->Add(fhRatioE);
477d6cee 462
123fc3bd 463 fhRatioPt = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
464 fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
465 outputContainer->Add(fhRatioPt);
466
467 fh2E = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 468 fh2E->SetXTitle("E_{rec} (GeV)");
469 fh2E->SetYTitle("E_{gen} (GeV)");
123fc3bd 470 outputContainer->Add(fh2E);
471
472 fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 473 fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
474 fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
123fc3bd 475 outputContainer->Add(fh2Pt);
476
c8fe2783 477 fhPtMCPhoton = new TH1F("hPtMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
478 fhPtMCPhoton->SetYTitle("N");
479 fhPtMCPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
480 outputContainer->Add(fhPtMCPhoton) ;
481
482 fhPhiMCPhoton = new TH2F
5ae09196 483 ("hPhiMCPhoton","#phi_{#gamma}, #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
c8fe2783 484 fhPhiMCPhoton->SetYTitle("#phi");
485 fhPhiMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
486 outputContainer->Add(fhPhiMCPhoton) ;
487
488 fhEtaMCPhoton = new TH2F
5ae09196 489 ("hEtaMCPhoton","#eta_{#gamma}, #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
c8fe2783 490 fhEtaMCPhoton->SetYTitle("#eta");
491 fhEtaMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
492 outputContainer->Add(fhEtaMCPhoton) ;
493
591cc579 494 fhPtPrompt = new TH1F("hPtMCPrompt","Number of prompt #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 495 fhPtPrompt->SetYTitle("N");
496 fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
497 outputContainer->Add(fhPtPrompt) ;
498
499 fhPhiPrompt = new TH2F
5ae09196 500 ("hPhiMCPrompt","#phi_{#gamma}, prompt #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 501 fhPhiPrompt->SetYTitle("#phi");
502 fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
503 outputContainer->Add(fhPhiPrompt) ;
504
505 fhEtaPrompt = new TH2F
5ae09196 506 ("hEtaMCPrompt","#eta_{#gamma}, prompt #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 507 fhEtaPrompt->SetYTitle("#eta");
508 fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
509 outputContainer->Add(fhEtaPrompt) ;
510
591cc579 511 fhPtFragmentation = new TH1F("hPtMCFragmentation","Number of fragmentation #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 512 fhPtFragmentation->SetYTitle("N");
513 fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
514 outputContainer->Add(fhPtFragmentation) ;
515
516 fhPhiFragmentation = new TH2F
5ae09196 517 ("hPhiMCFragmentation","#phi_{#gamma}, fragmentation #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 518 fhPhiFragmentation->SetYTitle("#phi");
519 fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
520 outputContainer->Add(fhPhiFragmentation) ;
521
522 fhEtaFragmentation = new TH2F
5ae09196 523 ("hEtaMCFragmentation","#eta_{#gamma}, fragmentation #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 524 fhEtaFragmentation->SetYTitle("#eta");
525 fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
526 outputContainer->Add(fhEtaFragmentation) ;
527
a3aebfff 528 fhPtISR = new TH1F("hPtMCISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 529 fhPtISR->SetYTitle("N");
530 fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
531 outputContainer->Add(fhPtISR) ;
532
533 fhPhiISR = new TH2F
a3aebfff 534 ("hPhiMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 535 fhPhiISR->SetYTitle("#phi");
536 fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
537 outputContainer->Add(fhPhiISR) ;
538
539 fhEtaISR = new TH2F
5ae09196 540 ("hEtaMCISR","#eta_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 541 fhEtaISR->SetYTitle("#eta");
542 fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
543 outputContainer->Add(fhEtaISR) ;
544
591cc579 545 fhPtPi0Decay = new TH1F("hPtMCPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 546 fhPtPi0Decay->SetYTitle("N");
547 fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
548 outputContainer->Add(fhPtPi0Decay) ;
549
550 fhPhiPi0Decay = new TH2F
5ae09196 551 ("hPhiMCPi0Decay","#phi_{#gamma}, #pi^{0} decay #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 552 fhPhiPi0Decay->SetYTitle("#phi");
553 fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
554 outputContainer->Add(fhPhiPi0Decay) ;
555
556 fhEtaPi0Decay = new TH2F
5ae09196 557 ("hEtaMCPi0Decay","#eta_{#gamma}, #pi^{0} #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 558 fhEtaPi0Decay->SetYTitle("#eta");
559 fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
560 outputContainer->Add(fhEtaPi0Decay) ;
561
a3aebfff 562 fhPtOtherDecay = new TH1F("hPtMCOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 563 fhPtOtherDecay->SetYTitle("N");
564 fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
565 outputContainer->Add(fhPtOtherDecay) ;
566
567 fhPhiOtherDecay = new TH2F
5ae09196 568 ("hPhiMCOtherDecay","#phi_{#gamma}, other decay #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 569 fhPhiOtherDecay->SetYTitle("#phi");
570 fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
571 outputContainer->Add(fhPhiOtherDecay) ;
572
573 fhEtaOtherDecay = new TH2F
5ae09196 574 ("hEtaMCOtherDecay","#eta_{#gamma}, other decay #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 575 fhEtaOtherDecay->SetYTitle("#eta");
576 fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
577 outputContainer->Add(fhEtaOtherDecay) ;
578
a3aebfff 579 fhPtConversion = new TH1F("hPtMCConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 580 fhPtConversion->SetYTitle("N");
581 fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
582 outputContainer->Add(fhPtConversion) ;
583
584 fhPhiConversion = new TH2F
5ae09196 585 ("hPhiMCConversion","#phi_{#gamma}, conversion #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 586 fhPhiConversion->SetYTitle("#phi");
587 fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
588 outputContainer->Add(fhPhiConversion) ;
589
590 fhEtaConversion = new TH2F
5ae09196 591 ("hEtaMCConversion","#eta_{#gamma}, conversion #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 592 fhEtaConversion->SetYTitle("#eta");
593 fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
594 outputContainer->Add(fhEtaConversion) ;
595
6175da48 596 fhEtaPhiConversion = new TH2F
597 ("hEtaPhiConversion","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
598 fhEtaPhiConversion->SetYTitle("#phi (rad)");
599 fhEtaPhiConversion->SetXTitle("#eta");
600 outputContainer->Add(fhEtaPhiConversion) ;
601
602 fhEtaPhi05Conversion = new TH2F
603 ("hEtaPhi05Conversion","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
604 fhEtaPhi05Conversion->SetYTitle("#phi (rad)");
605 fhEtaPhi05Conversion->SetXTitle("#eta");
606 outputContainer->Add(fhEtaPhi05Conversion) ;
607
608 fhPtAntiNeutron = new TH1F("hPtMCAntiNeutron","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
609 fhPtAntiNeutron->SetYTitle("N");
610 fhPtAntiNeutron->SetXTitle("p_{T #gamma}(GeV/c)");
611 outputContainer->Add(fhPtAntiNeutron) ;
612
613 fhPhiAntiNeutron = new TH2F
614 ("hPhiMCAntiNeutron","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
615 fhPhiAntiNeutron->SetYTitle("#phi");
616 fhPhiAntiNeutron->SetXTitle("p_{T #gamma} (GeV/c)");
617 outputContainer->Add(fhPhiAntiNeutron) ;
618
619 fhEtaAntiNeutron = new TH2F
620 ("hEtaMCAntiNeutron","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
621 fhEtaAntiNeutron->SetYTitle("#eta");
622 fhEtaAntiNeutron->SetXTitle("p_{T #gamma} (GeV/c)");
623 outputContainer->Add(fhEtaAntiNeutron) ;
624
41121cfe 625 fhPtAntiProton = new TH1F("hPtMCAntiProton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
6175da48 626 fhPtAntiProton->SetYTitle("N");
627 fhPtAntiProton->SetXTitle("p_{T #gamma}(GeV/c)");
628 outputContainer->Add(fhPtAntiProton) ;
629
630 fhPhiAntiProton = new TH2F
631 ("hPhiMCAntiProton","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
632 fhPhiAntiProton->SetYTitle("#phi");
633 fhPhiAntiProton->SetXTitle("p_{T #gamma} (GeV/c)");
634 outputContainer->Add(fhPhiAntiProton) ;
635
636 fhEtaAntiProton = new TH2F
637 ("hEtaMCAntiProton","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
638 fhEtaAntiProton->SetYTitle("#eta");
639 fhEtaAntiProton->SetXTitle("p_{T #gamma} (GeV/c)");
640 outputContainer->Add(fhEtaAntiProton) ;
641
a3aebfff 642 fhPtUnknown = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 643 fhPtUnknown->SetYTitle("N");
644 fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
645 outputContainer->Add(fhPtUnknown) ;
646
647 fhPhiUnknown = new TH2F
5ae09196 648 ("hPhiMCUnknown","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 649 fhPhiUnknown->SetYTitle("#phi");
650 fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
651 outputContainer->Add(fhPhiUnknown) ;
652
653 fhEtaUnknown = new TH2F
5ae09196 654 ("hEtaMCUnknown","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 655 fhEtaUnknown->SetYTitle("#eta");
656 fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
657 outputContainer->Add(fhEtaUnknown) ;
a3aebfff 658
20218aea 659 if(fCheckConversion){
660 fhPtConversionTagged = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
661 fhPtConversionTagged->SetYTitle("N");
662 fhPtConversionTagged->SetXTitle("p_{T #gamma}(GeV/c)");
663 outputContainer->Add(fhPtConversionTagged) ;
664
665
666 fhPtAntiNeutronTagged = new TH1F("hPtMCAntiNeutronTagged","Number of AntiNeutron id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
667 fhPtAntiNeutronTagged->SetYTitle("N");
668 fhPtAntiNeutronTagged->SetXTitle("p_{T #gamma}(GeV/c)");
669 outputContainer->Add(fhPtAntiNeutronTagged) ;
670
671 fhPtAntiProtonTagged = new TH1F("hPtMCAntiProtonTagged","Number of AntiProton id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
672 fhPtAntiProtonTagged->SetYTitle("N");
673 fhPtAntiProtonTagged->SetXTitle("p_{T #gamma}(GeV/c)");
674 outputContainer->Add(fhPtAntiProtonTagged) ;
675
676 fhPtUnknownTagged = new TH1F("hPtMCUnknownTagged","Number of Unknown id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
677 fhPtUnknownTagged->SetYTitle("N");
678 fhPtUnknownTagged->SetXTitle("p_{T #gamma}(GeV/c)");
679 outputContainer->Add(fhPtUnknownTagged) ;
680
681 fhConvDeltaEtaMCConversion = new TH2F
682 ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5);
683 fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
684 fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
685 outputContainer->Add(fhConvDeltaEtaMCConversion) ;
686
687 fhConvDeltaPhiMCConversion = new TH2F
688 ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5);
689 fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
690 fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
691 outputContainer->Add(fhConvDeltaPhiMCConversion) ;
692
693 fhConvDeltaEtaPhiMCConversion = new TH2F
694 ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5);
695 fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
696 fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
697 outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
698
699 fhConvAsymMCConversion = new TH2F
700 ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1);
701 fhConvAsymMCConversion->SetYTitle("Asymmetry");
702 fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
703 outputContainer->Add(fhConvAsymMCConversion) ;
704
705 fhConvPtMCConversion = new TH2F
706 ("hConvPtMCConversion","p_{T} of selected conversion pairs from real conversions",100,0,fMassCut,100,0.,10.);
707 fhConvPtMCConversion->SetYTitle("Pair p_{T} (GeV/c)");
708 fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
709 outputContainer->Add(fhConvPtMCConversion) ;
710
711 fhConvDispersionMCConversion = new TH2F
712 ("hConvDispersionMCConversion","p_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.);
713 fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
714 fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
715 outputContainer->Add(fhConvDispersionMCConversion) ;
716
717 fhConvM02MCConversion = new TH2F
718 ("hConvM02MCConversion","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
719 fhConvM02MCConversion->SetYTitle("M02 cluster 1");
720 fhConvM02MCConversion->SetXTitle("M02 cluster 2");
721 outputContainer->Add(fhConvM02MCConversion) ;
722
723 fhConvDeltaEtaMCAntiNeutron = new TH2F
724 ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5);
725 fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
726 fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
727 outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
728
729 fhConvDeltaPhiMCAntiNeutron = new TH2F
730 ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5);
731 fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
732 fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
733 outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
734
735 fhConvDeltaEtaPhiMCAntiNeutron = new TH2F
736 ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
737 fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
738 fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
739 outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;
740
741 fhConvAsymMCAntiNeutron = new TH2F
742 ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1);
743 fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
744 fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
745 outputContainer->Add(fhConvAsymMCAntiNeutron) ;
746
747 fhConvPtMCAntiNeutron = new TH2F
748 ("hConvPtMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0.,10.);
749 fhConvPtMCAntiNeutron->SetYTitle("Pair p_{T} (GeV/c)");
750 fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
751 outputContainer->Add(fhConvPtMCAntiNeutron) ;
752
753 fhConvDispersionMCAntiNeutron = new TH2F
754 ("hConvDispersionMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.);
755 fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
756 fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
757 outputContainer->Add(fhConvDispersionMCAntiNeutron) ;
758
759 fhConvM02MCAntiNeutron = new TH2F
760 ("hConvM02MCAntiNeutron","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
761 fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
762 fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
763 outputContainer->Add(fhConvM02MCAntiNeutron) ;
764
765 fhConvDeltaEtaMCAntiProton = new TH2F
766 ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5);
767 fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
768 fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
769 outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
770
771 fhConvDeltaPhiMCAntiProton = new TH2F
772 ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5);
773 fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
774 fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
775 outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
776
777 fhConvDeltaEtaPhiMCAntiProton = new TH2F
778 ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
779 fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
780 fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
781 outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;
782
783 fhConvAsymMCAntiProton = new TH2F
784 ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1);
785 fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
786 fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
787 outputContainer->Add(fhConvAsymMCAntiProton) ;
788
789 fhConvPtMCAntiProton = new TH2F
790 ("hConvPtMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,100,0.,10.);
791 fhConvPtMCAntiProton->SetYTitle("Pair p_{T} (GeV/c)");
792 fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
793 outputContainer->Add(fhConvPtMCAntiProton) ;
794
795 fhConvDispersionMCAntiProton = new TH2F
796 ("hConvDispersionMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.);
797 fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
798 fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
799 outputContainer->Add(fhConvDispersionMCAntiProton) ;
800
801 fhConvM02MCAntiProton = new TH2F
802 ("hConvM02MCAntiProton","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
803 fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
804 fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
805 outputContainer->Add(fhConvM02MCAntiProton) ;
806
807 fhConvDeltaEtaMCString = new TH2F
808 ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5);
809 fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
810 fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/c^2)");
811 outputContainer->Add(fhConvDeltaEtaMCString) ;
812
813 fhConvDeltaPhiMCString = new TH2F
814 ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5);
815 fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
816 fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/c^2)");
817 outputContainer->Add(fhConvDeltaPhiMCString) ;
818
819 fhConvDeltaEtaPhiMCString = new TH2F
820 ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5);
821 fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
822 fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
823 outputContainer->Add(fhConvDeltaEtaPhiMCString) ;
824
825 fhConvAsymMCString = new TH2F
826 ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1);
827 fhConvAsymMCString->SetYTitle("Asymmetry");
828 fhConvAsymMCString->SetXTitle("Pair Mass (GeV/c^2)");
829 outputContainer->Add(fhConvAsymMCString) ;
830
831 fhConvPtMCString = new TH2F
832 ("hConvPtMCString","p_{T} of selected conversion pairs from string",100,0,fMassCut,100,0.,10.);
833 fhConvPtMCString->SetYTitle("Pair p_{T} (GeV/c)");
834 fhConvPtMCString->SetXTitle("Pair Mass (GeV/c^2)");
835 outputContainer->Add(fhConvPtMCString) ;
836
837 fhConvDispersionMCString = new TH2F
838 ("hConvDispersionMCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
839 fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
840 fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
841 outputContainer->Add(fhConvDispersionMCString) ;
842
843 fhConvM02MCString = new TH2F
844 ("hConvM02MCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
845 fhConvM02MCString->SetYTitle("M02 cluster 1");
846 fhConvM02MCString->SetXTitle("M02 cluster 2");
c4a7d28a 847 outputContainer->Add(fhConvM02MCString) ;
848
849 fhConvDistMCConversion = new TH2F
850 ("hConvDistMCConversion","calculated conversion distance vs real vertes for MC conversion",100,0.,5.,100,0.,5.);
851 fhConvDistMCConversion->SetYTitle("distance");
852 fhConvDistMCConversion->SetXTitle("vertex R");
853 outputContainer->Add(fhConvDistMCConversion) ;
854
855 fhConvDistMCConversionCuts = new TH2F
856 ("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.);
857 fhConvDistMCConversionCuts->SetYTitle("distance");
858 fhConvDistMCConversionCuts->SetXTitle("vertex R");
859 outputContainer->Add(fhConvDistMCConversionCuts) ;
860
20218aea 861 }
6175da48 862
477d6cee 863 }//Histos with MC
0c1383b5 864
d39cba7e 865 //Store calo PID histograms
866 if(fRejectTrackMatch){
867 TList * caloPIDHistos = GetCaloPID()->GetCreateOutputObjects() ;
868 for(Int_t i = 0; i < caloPIDHistos->GetEntries(); i++) {
869 outputContainer->Add(caloPIDHistos->At(i)) ;
870 }
871 delete caloPIDHistos;
872 }
873
477d6cee 874 return outputContainer ;
875
1c5acb87 876}
877
6639984f 878//____________________________________________________________________________
879void AliAnaPhoton::Init()
880{
881
882 //Init
883 //Do some checks
1e86c71e 884 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
591cc579 885 printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 886 abort();
887 }
1e86c71e 888 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
591cc579 889 printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 890 abort();
891 }
892
893}
894
895
1c5acb87 896//____________________________________________________________________________
897void AliAnaPhoton::InitParameters()
898{
899
900 //Initialize the parameters of the analysis.
a3aebfff 901 AddToHistogramsName("AnaPhoton_");
902
6175da48 903 fCalorimeter = "EMCAL" ;
904 fMinDist = 2.;
905 fMinDist2 = 4.;
906 fMinDist3 = 5.;
907 fMassCut = 0.03; //30 MeV
1e86c71e 908
4cf55759 909 fTimeCutMin = -1;
910 fTimeCutMax = 9999999;
6175da48 911 fNCellsCut = 0;
2ac125bf 912
1e86c71e 913 fRejectTrackMatch = kTRUE ;
914 fCheckConversion = kFALSE;
20218aea 915 fRemoveConvertedPair = kFALSE;
1e86c71e 916 fAddConvertedPairsToAOD = kFALSE;
917
1c5acb87 918}
919
920//__________________________________________________________________
921void AliAnaPhoton::MakeAnalysisFillAOD()
922{
f8006433 923 //Do photon analysis and fill aods
f37fa8d2 924
6175da48 925 //Get the vertex
5025c139 926 Double_t v[3] = {0,0,0}; //vertex ;
927 GetReader()->GetVertex(v);
f8006433 928
f37fa8d2 929 //Select the Calorimeter of the photon
c8fe2783 930 TObjArray * pl = 0x0;
477d6cee 931 if(fCalorimeter == "PHOS")
be518ab0 932 pl = GetPHOSClusters();
477d6cee 933 else if (fCalorimeter == "EMCAL")
be518ab0 934 pl = GetEMCALClusters();
5ae09196 935
936 if(!pl) {
937 Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
938 return;
939 }
0ae57829 940
6175da48 941 //Init arrays, variables, get number of clusters
1e86c71e 942 TLorentzVector mom, mom2 ;
943 Int_t nCaloClusters = pl->GetEntriesFast();
6175da48 944 //List to be used in conversion analysis, to tag the cluster as candidate for conversion
20218aea 945 Bool_t * indexConverted = 0x0;
946 if(fCheckConversion){
947 indexConverted = new Bool_t[nCaloClusters];
948 for (Int_t i = 0; i < nCaloClusters; i++)
949 indexConverted[i] = kFALSE;
950 }
951
6175da48 952 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
953
954 //----------------------------------------------------
955 // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
956 //----------------------------------------------------
957 // Loop on clusters
1e86c71e 958 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
959
0ae57829 960 AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
961 //printf("calo %d, %f\n",icalo,calo->E());
c4a7d28a 962
f8006433 963 //Get the index where the cluster comes, to retrieve the corresponding vertex
c8fe2783 964 Int_t evtIndex = 0 ;
965 if (GetMixedEvent()) {
966 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
5025c139 967 //Get the vertex and check it is not too large in z
96539743 968 if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
c8fe2783 969 }
f8006433 970
f37fa8d2 971 //Cluster selection, not charged, with photon id and in fiducial cut
233e0df8 972
f37fa8d2 973 //Input from second AOD?
f8006433 974 //Int_t input = 0;
be518ab0 975 // if (fCalorimeter == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= icalo)
f37fa8d2 976 // input = 1 ;
be518ab0 977 // else if(fCalorimeter == "PHOS" && GetReader()->GetPHOSClustersNormalInputEntries() <= icalo)
f37fa8d2 978 // input = 1;
233e0df8 979
f37fa8d2 980 //Get Momentum vector,
f8006433 981 //if (input == 0)
982 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
983 calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
984 else{
985 Double_t vertex[]={0,0,0};
986 calo->GetMomentum(mom,vertex) ;
987 }
f8006433 988
f37fa8d2 989 // else if(input == 1)
990 // calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
c8fe2783 991
6175da48 992 //--------------------------------------
993 // Cluster selection
994 //--------------------------------------
477d6cee 995
c4a7d28a 996 if(!ClusterSelected(calo,mom)) continue;
6175da48 997
998 //----------------------------
999 //Create AOD for analysis
1000 //----------------------------
1001 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
1002
1003 //...............................................
1004 //Set the indeces of the original caloclusters (MC, ID), and calorimeter
1005 Int_t label = calo->GetLabel();
1006 aodph.SetLabel(label);
c4a7d28a 1007 //aodph.SetInputFileIndex(input);
6175da48 1008 aodph.SetCaloLabel(calo->GetID(),-1);
1009 aodph.SetDetector(fCalorimeter);
c4a7d28a 1010 //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
6175da48 1011
1012 //...............................................
1013 //Set bad channel distance bit
c4a7d28a 1014 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
f37fa8d2 1015 if (distBad > fMinDist3) aodph.SetDistToBad(2) ;
477d6cee 1016 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
f37fa8d2 1017 else aodph.SetDistToBad(0) ;
af7b3903 1018 //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
c8fe2783 1019
6175da48 1020 //...............................................
1021 //Set number of cells in this cluster
1022 //Temporary patch FIXME
1023 aodph.SetBtag(calo->GetNCells());
1024 // MEFIX
1025
1026 //-------------------------------------
f37fa8d2 1027 //PID selection or bit setting
6175da48 1028 //-------------------------------------
1029 // MC
477d6cee 1030 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
f37fa8d2 1031 //Get most probable PID, check PID weights (in MC this option is mandatory)
21a4b1c0 1032 aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
1033 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
f37fa8d2 1034 //If primary is not photon, skip it.
21a4b1c0 1035 if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
6175da48 1036 }
1037 //...............................................
1038 // Data, PID check on
477d6cee 1039 else if(IsCaloPIDOn()){
f37fa8d2 1040 //Get most probable PID, 2 options check PID weights
1041 //or redo PID, recommended option for EMCal.
477d6cee 1042 if(!IsCaloPIDRecalculationOn())
21a4b1c0 1043 aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
477d6cee 1044 else
21a4b1c0 1045 aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
477d6cee 1046
21a4b1c0 1047 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
477d6cee 1048
f37fa8d2 1049 //If cluster does not pass pid, not photon, skip it.
21a4b1c0 1050 if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
477d6cee 1051
1052 }
6175da48 1053 //...............................................
1054 // Data, PID check off
477d6cee 1055 else{
f37fa8d2 1056 //Set PID bits for later selection (AliAnaPi0 for example)
1057 //GetPDG already called in SetPIDBits.
f2ccb5b8 1058 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
a3aebfff 1059 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
477d6cee 1060 }
1061
21a4b1c0 1062 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetIdentifiedParticleType());
477d6cee 1063
6175da48 1064 //--------------------------------------------------------------------------------------
f37fa8d2 1065 //Play with the MC stack if available
6175da48 1066 //--------------------------------------------------------------------------------------
1067
f37fa8d2 1068 //Check origin of the candidates
477d6cee 1069 if(IsDataMC()){
c8fe2783 1070 aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
fb516de2 1071 if(GetDebug() > 0)
1072 printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
477d6cee 1073 }//Work with stack also
1074
6175da48 1075 //--------------------------------------------------------------------------------------
1076 // Conversions pairs analysis
f37fa8d2 1077 // Check if cluster comes from a conversion in the material in front of the calorimeter
1078 // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
6175da48 1079 //--------------------------------------------------------------------------------------
1080
1081 // Do analysis only if there are more than one cluster
20218aea 1082 if( nCaloClusters > 1 && fCheckConversion){
c8fe2783 1083 Bool_t bConverted = kFALSE;
1084 Int_t id2 = -1;
1e86c71e 1085
f37fa8d2 1086 //Check if set previously as converted couple, if so skip its use.
20218aea 1087 if (indexConverted[icalo]) continue;
1e86c71e 1088
6175da48 1089 // Second cluster loop
c8fe2783 1090 for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
f37fa8d2 1091 //Check if set previously as converted couple, if so skip its use.
c8fe2783 1092 if (indexConverted[jcalo]) continue;
f37fa8d2 1093 //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
0ae57829 1094 AliVCluster * calo2 = (AliVCluster*) (pl->At(jcalo)); //Get cluster kinematics
6175da48 1095
c4a7d28a 1096
6175da48 1097 //Mixed event, get index of event
c8fe2783 1098 Int_t evtIndex2 = 0 ;
1099 if (GetMixedEvent()) {
1100 evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
f8006433 1101
6175da48 1102 }
1103
1104 //Get kinematics of second cluster
f8006433 1105 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
8cdc266d 1106 calo2->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
f8006433 1107 else{
1108 Double_t vertex[]={0,0,0};
8cdc266d 1109 calo2->GetMomentum(mom2,vertex) ;
f8006433 1110 }
1111
c4a7d28a 1112 //--------------------------------------
1113 // Cluster selection
1114 //--------------------------------------
1115
1116 if(!ClusterSelected(calo2,mom2)) continue;
1117
1118 //................................................
1119 // Get TOF of each cluster in pair, calculate difference if small,
1120 // take this pair. Way to reject clusters from hadrons (or pileup?)
1121 Double_t t12diff = calo2->GetTOF()-calo->GetTOF()*1e9;
1122 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
c8fe2783 1123
6175da48 1124 //................................................
f37fa8d2 1125 //Get mass of pair, if small, take this pair.
41121cfe 1126 Float_t pairM = (mom+mom2).M();
1127 //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
1128 if(pairM < fMassCut){
1129 aodph.SetTagged(kFALSE);
c8fe2783 1130 id2 = calo2->GetID();
6175da48 1131 indexConverted[icalo]=kTRUE;
c4a7d28a 1132 indexConverted[jcalo]=kTRUE;
41121cfe 1133 Float_t asymmetry = TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E());
1134 Float_t dPhi = mom.Phi()-mom2.Phi();
d39cba7e 1135 Float_t dEta = mom.Eta()-mom2.Eta();
41121cfe 1136
c4a7d28a 1137 //...............................................
1138 //Fill few histograms with kinematics of the pair
1139 //FIXME, move all this to MakeAnalysisFillHistograms ...
1140
1141 fhConvDeltaEta ->Fill( pairM, dPhi );
1142 fhConvDeltaPhi ->Fill( pairM, dEta );
1143 fhConvAsym ->Fill( pairM, asymmetry );
1144 fhConvDeltaEtaPhi->Fill( dEta , dPhi );
1145 fhConvPt ->Fill( pairM, (mom+mom2).Pt());
1146
d39cba7e 1147 //Estimate conversion distance, T. Awes, M. Ivanov
1148 //Under the assumption that the pair has zero mass, and that each electron
1149 //of the pair has the same momentum, they will each have the same bend radius
1150 //given by R=p/(qB) = p / (300 B) with p in [MeV/c], B in [Tesla] and R in [m].
1151 //With nominal ALICE magnet current of 30kA B=0.5T, and so with E_cluster=p,
1152 //R = E/1.5 [cm]. Under these assumptions, the distance from the conversion
1153 //point to the EMCal can be related to the separation distance, L=2y, on the EMCal
1154 //as d = sqrt(R^2 -(R-y)^2) = sqrt(2Ry - y^2). And since R>>y we can write as
1155 //d = sqrt(E*L/1.5) where E is the cluster energy and L is the distance in cm between
1156 //the clusters.
c4a7d28a 1157 Float_t pos1[3];
1158 calo->GetPosition(pos1);
1159 Float_t pos2[3];
1160 calo2->GetPosition(pos2);
1161 Float_t clustDist = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0])+
1162 (pos1[1]-pos2[1])*(pos1[1]-pos2[1])+
1163 (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
1164
1165 Float_t convDist = TMath::Sqrt(mom.E() *clustDist*0.01/0.15);
1166 Float_t convDist2 = TMath::Sqrt(mom2.E()*clustDist*0.01/0.15);
1167 //printf("l = %f, e1 = %f, d1=%f, e2 = %f, d2=%f\n",clustDist,mom.E(),convDist,mom2.E(),convDist2);
8cdc266d 1168 if(GetDebug() > 2)
41121cfe 1169 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",
1170 pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
8cdc266d 1171 calo->GetID(),calo->E(),GetCaloUtils()->GetModuleNumber(calo), mom.Eta(), mom.Phi(),
1172 id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2),mom2.Eta(), mom2.Phi());
6175da48 1173
c4a7d28a 1174 fhConvDistEta ->Fill(mom .Eta(), convDist );
1175 fhConvDistEta ->Fill(mom2.Eta(), convDist2);
1176 fhConvDistEn ->Fill(mom .E(), convDist );
1177 fhConvDistEn ->Fill(mom2.E(), convDist2);
1178 fhConvDistMass->Fill((mom+mom2).M(), convDist );
1179 //dEta cut
1180 if(dEta<0.05){
1181 fhConvDistEtaCutEta ->Fill(mom .Eta(), convDist );
1182 fhConvDistEtaCutEta ->Fill(mom2.Eta(), convDist2);
1183 fhConvDistEnCutEta ->Fill(mom .E(), convDist );
1184 fhConvDistEnCutEta ->Fill(mom2.E(), convDist2);
1185 fhConvDistMassCutEta->Fill((mom+mom2).M(), convDist );
1186 //mass cut
1187 if(pairM<0.01){//10 MeV
1188 fhConvDistEtaCutMass ->Fill(mom .Eta(), convDist );
1189 fhConvDistEtaCutMass ->Fill(mom2.Eta(), convDist2);
1190 fhConvDistEnCutMass ->Fill(mom .E(), convDist );
1191 fhConvDistEnCutMass ->Fill(mom2.E(), convDist2);
1192 // asymmetry cut
1193 if(asymmetry<0.1){
1194 fhConvDistEtaCutAsy ->Fill(mom .Eta(), convDist );
1195 fhConvDistEtaCutAsy ->Fill(mom2.Eta(), convDist2);
1196 fhConvDistEnCutAsy ->Fill(mom .E(), convDist );
1197 fhConvDistEnCutAsy ->Fill(mom2.E(), convDist2);
1198 }//asymmetry cut
1199 }//mass cut
1200 }//dEta cut
6175da48 1201
1202 //...............................................
1203 //Select pairs in a eta-phi window
41121cfe 1204 if(TMath::Abs(dEta) < fConvDEtaCut &&
1205 TMath::Abs(dPhi) < fConvDPhiMaxCut &&
1206 TMath::Abs(dPhi) > fConvDPhiMinCut &&
1207 asymmetry < fConvAsymCut ){
1208 bConverted = kTRUE;
1209 }
1210 //printf("Accepted? %d\n",bConverted);
6175da48 1211 //...........................................
1212 //Fill more histograms, simulated data
1213 //FIXME, move all this to MakeAnalysisFillHistograms ...
1214 if(IsDataMC()){
1215
1216 //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
1217 Int_t ancPDG = 0;
1218 Int_t ancStatus = 0;
1219 TLorentzVector momentum;
c4a7d28a 1220 TVector3 prodVertex;
6175da48 1221 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(calo->GetLabel(), calo2->GetLabel(),
c4a7d28a 1222 GetReader(), ancPDG, ancStatus, momentum, prodVertex);
6175da48 1223
1224 // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1225 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1226
1227 Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(calo2->GetLabels(),calo2->GetNLabels(),GetReader(), 0);
1228 if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCConversion)){
1229 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
41121cfe 1230 fhConvDeltaEtaMCConversion ->Fill( pairM, dEta );
1231 fhConvDeltaPhiMCConversion ->Fill( pairM, dPhi );
1232 fhConvAsymMCConversion ->Fill( pairM, asymmetry );
1233 fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi );
1234 fhConvPtMCConversion ->Fill( pairM, (mom+mom2).Pt());
6175da48 1235 fhConvDispersionMCConversion ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1236 fhConvM02MCConversion ->Fill( calo->GetM02(), calo2->GetM02());
c4a7d28a 1237 fhConvDistMCConversion ->Fill( convDist , prodVertex.Mag() );
1238 fhConvDistMCConversion ->Fill( convDist2, prodVertex.Mag() );
1239
1240 if(dEta<0.05 && pairM<0.01 && asymmetry<0.1){
1241 fhConvDistMCConversionCuts->Fill( convDist , prodVertex.Mag() );
1242 fhConvDistMCConversionCuts->Fill( convDist2, prodVertex.Mag() );
1243 }
6175da48 1244
1245 }
1246 }
1247 else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiNeutron)){
1248 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
41121cfe 1249 fhConvDeltaEtaMCAntiNeutron ->Fill( pairM, dEta );
1250 fhConvDeltaPhiMCAntiNeutron ->Fill( pairM, dPhi );
1251 fhConvAsymMCAntiNeutron ->Fill( pairM, asymmetry );
1252 fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi );
1253 fhConvPtMCAntiNeutron ->Fill( pairM, (mom+mom2).Pt());
6175da48 1254 fhConvDispersionMCAntiNeutron ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1255 fhConvM02MCAntiNeutron ->Fill( calo->GetM02(), calo2->GetM02());
1256 }
1257 }
1258 else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiProton)){
1259 if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
41121cfe 1260 fhConvDeltaEtaMCAntiProton ->Fill( pairM, dEta );
1261 fhConvDeltaPhiMCAntiProton ->Fill( pairM, dPhi );
1262 fhConvAsymMCAntiProton ->Fill( pairM, asymmetry );
1263 fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi );
1264 fhConvPtMCAntiProton ->Fill( pairM, (mom+mom2).Pt());
6175da48 1265 fhConvDispersionMCAntiProton ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1266 fhConvM02MCAntiProton ->Fill( calo->GetM02(), calo2->GetM02());
1267 }
1268 }
1269
1270 //Pairs coming from fragmenting pairs.
1271 if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
41121cfe 1272 fhConvDeltaEtaMCString ->Fill( pairM, dPhi);
1273 fhConvDeltaPhiMCString ->Fill( pairM, dPhi);
1274 fhConvAsymMCString ->Fill( pairM, TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
1275 fhConvDeltaEtaPhiMCString ->Fill( dPhi, dPhi );
1276 fhConvPtMCString ->Fill( pairM, (mom+mom2).Pt());
6175da48 1277 fhConvDispersionMCString ->Fill( calo->GetDispersion(), calo2->GetDispersion());
1278 fhConvM02MCString ->Fill( calo->GetM02(), calo2->GetM02());
1279 }
1280
1281 }// Data MC
1282
c8fe2783 1283 break;
1284 }
1e86c71e 1285
c8fe2783 1286 }//Mass loop
1e86c71e 1287
6175da48 1288 //..........................................................................................................
1289 //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
c8fe2783 1290 if(bConverted){
6175da48 1291 //Add to AOD
c8fe2783 1292 if(fAddConvertedPairsToAOD){
f37fa8d2 1293 //Create AOD of pair analysis
c8fe2783 1294 TLorentzVector mpair = mom+mom2;
1295 AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
1296 aodpair.SetLabel(aodph.GetLabel());
f8006433 1297 //aodpair.SetInputFileIndex(input);
c8fe2783 1298
f37fa8d2 1299 //printf("Index %d, Id %d\n",icalo, calo->GetID());
1300 //Set the indeces of the original caloclusters
c8fe2783 1301 aodpair.SetCaloLabel(calo->GetID(),id2);
1302 aodpair.SetDetector(fCalorimeter);
21a4b1c0 1303 aodpair.SetIdentifiedParticleType(aodph.GetIdentifiedParticleType());
c8fe2783 1304 aodpair.SetTag(aodph.GetTag());
6175da48 1305 aodpair.SetTagged(kTRUE);
f37fa8d2 1306 //Add AOD with pair object to aod branch
c8fe2783 1307 AddAODParticle(aodpair);
f37fa8d2 1308 //printf("\t \t both added pair\n");
c8fe2783 1309 }
1310
f37fa8d2 1311 //Do not add the current calocluster
20218aea 1312 if(fRemoveConvertedPair) continue;
6175da48 1313 else {
1314 //printf("TAGGED\n");
1315 //Tag this cluster as likely conversion
1316 aodph.SetTagged(kTRUE);
1317 }
c8fe2783 1318 }//converted pair
1319 }//check conversion
f37fa8d2 1320 //printf("\t \t added single cluster %d\n",icalo);
1e86c71e 1321
6175da48 1322 //FIXME, this to MakeAnalysisFillHistograms ...
c4a7d28a 1323 fhNCellsE->Fill(aodph.E(),calo->GetNCells());
1324 if(fFillSSHistograms){
1325 if(calo->E()<2){
1326 fhNCellsLam0LowE->Fill(calo->GetNCells(),calo->GetM02());
1327 fhNCellsLam1LowE->Fill(calo->GetNCells(),calo->GetM20());
1328 fhNCellsDispLowE->Fill(calo->GetNCells(),calo->GetDispersion());
1329 fhLam1Lam0LowE ->Fill(calo->GetM20(), calo->GetM02());
1330 }
1331 else {
1332 fhNCellsLam0HighE->Fill(calo->GetNCells(),calo->GetM02());
1333 fhNCellsLam1HighE->Fill(calo->GetNCells(),calo->GetM20());
1334 fhNCellsDispHighE->Fill(calo->GetNCells(),calo->GetDispersion());
1335 fhLam1Lam0HighE ->Fill(calo->GetM20(), calo->GetM02());
1336 }
1337
1338 fhEtaLam0->Fill(aodph.Eta(), calo->GetM02());
1339 fhPhiLam0->Fill(aodph.Phi(), calo->GetM02());
1340 fhLam0E ->Fill(aodph.E(), calo->GetM02());
1341 fhLam1E ->Fill(aodph.E(), calo->GetM20());
1342 }
6175da48 1343
f37fa8d2 1344 //Add AOD with photon object to aod branch
477d6cee 1345 AddAODParticle(aodph);
1346
1347 }//loop
1348
4a745797 1349 delete [] indexConverted;
1350
f37fa8d2 1351 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
477d6cee 1352
1c5acb87 1353}
1354
1355//__________________________________________________________________
1356void AliAnaPhoton::MakeAnalysisFillHistograms()
1357{
6175da48 1358 //Fill histograms
f8006433 1359
6175da48 1360 //-------------------------------------------------------------------
591cc579 1361 // Access MC information in stack if requested, check that it exists.
1362 AliStack * stack = 0x0;
1363 TParticle * primary = 0x0;
1364 TClonesArray * mcparticles0 = 0x0;
f8006433 1365 //TClonesArray * mcparticles1 = 0x0;
591cc579 1366 AliAODMCParticle * aodprimary = 0x0;
1367 if(IsDataMC()){
1368
1369 if(GetReader()->ReadStack()){
1370 stack = GetMCStack() ;
1371 if(!stack) {
1372 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1373 abort();
1374 }
f8006433 1375
591cc579 1376 }
1377 else if(GetReader()->ReadAODMCParticles()){
f8006433 1378
591cc579 1379 //Get the list of MC particles
1380 mcparticles0 = GetReader()->GetAODMCParticles(0);
1381 if(!mcparticles0 && GetDebug() > 0) {
1382 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
1383 }
f8006433 1384 // if(GetReader()->GetSecondInputAODTree()){
1385 // mcparticles1 = GetReader()->GetAODMCParticles(1);
1386 // if(!mcparticles1 && GetDebug() > 0) {
1387 // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
1388 // }
1389 // }
591cc579 1390
1391 }
1392 }// is data and MC
1393
6175da48 1394
1395 // Get vertex
2244659d 1396 Double_t v[3] = {0,0,0}; //vertex ;
1397 GetReader()->GetVertex(v);
6175da48 1398 //fhVertex->Fill(v[0],v[1],v[2]);
1399 if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1400
1401 //----------------------------------
1402 //Loop on stored AOD photons
123fc3bd 1403 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2244659d 1404 fhNtraNclu->Fill(GetReader()->GetTrackMultiplicity(), naod);
123fc3bd 1405 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1406
1407 for(Int_t iaod = 0; iaod < naod ; iaod++){
c8fe2783 1408 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
21a4b1c0 1409 Int_t pdg = ph->GetIdentifiedParticleType();
c8fe2783 1410
1411 if(GetDebug() > 3)
21a4b1c0 1412 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
c8fe2783 1413
1414 //If PID used, fill histos with photons in Calorimeter fCalorimeter
1415 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
1416 if(ph->GetDetector() != fCalorimeter) continue;
1417
1418 if(GetDebug() > 2)
1419 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1420
6175da48 1421 //................................
c8fe2783 1422 //Fill photon histograms
1423 Float_t ptcluster = ph->Pt();
1424 Float_t phicluster = ph->Phi();
1425 Float_t etacluster = ph->Eta();
1426 Float_t ecluster = ph->E();
1427
20218aea 1428 fhEPhoton ->Fill(ecluster);
c8fe2783 1429 fhPtPhoton ->Fill(ptcluster);
1430 fhPhiPhoton ->Fill(ptcluster,phicluster);
20218aea 1431 fhEtaPhoton ->Fill(ptcluster,etacluster);
1432 if(ecluster > 0.5) fhEtaPhiPhoton ->Fill(etacluster, phicluster);
1433 else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
1434
1435 if(fCheckConversion &&ph->IsTagged()){
1436 fhPtPhotonConv->Fill(ptcluster);
1437 if(ecluster > 0.5) fhEtaPhiPhotonConv ->Fill(etacluster, phicluster);
1438 else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
6175da48 1439 }
20218aea 1440
6175da48 1441 //.......................................
c8fe2783 1442 //Play with the MC data if available
1443 if(IsDataMC()){
1444
1445 Int_t tag =ph->GetTag();
fb516de2 1446
c8fe2783 1447 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
f8006433 1448 {
1449 fhPtMCPhoton ->Fill(ptcluster);
1450 fhPhiMCPhoton ->Fill(ptcluster,phicluster);
1451 fhEtaMCPhoton ->Fill(ptcluster,etacluster);
1452
1453 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
1454 {
1455 fhPtConversion ->Fill(ptcluster);
1456 fhPhiConversion ->Fill(ptcluster,phicluster);
1457 fhEtaConversion ->Fill(ptcluster,etacluster);
6175da48 1458 if(ph->IsTagged()) fhPtConversionTagged ->Fill(ptcluster);
1459 if(ptcluster > 0.5)fhEtaPhiConversion ->Fill(etacluster,phicluster);
1460 else fhEtaPhi05Conversion ->Fill(etacluster,phicluster);
f8006433 1461 }
1462
1463 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
1464 fhPtPrompt ->Fill(ptcluster);
1465 fhPhiPrompt ->Fill(ptcluster,phicluster);
1466 fhEtaPrompt ->Fill(ptcluster,etacluster);
1467 }
1468 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
1469 {
1470 fhPtFragmentation ->Fill(ptcluster);
1471 fhPhiFragmentation ->Fill(ptcluster,phicluster);
1472 fhEtaFragmentation ->Fill(ptcluster,etacluster);
1473 }
1474 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
1475 {
1476 fhPtISR ->Fill(ptcluster);
1477 fhPhiISR ->Fill(ptcluster,phicluster);
1478 fhEtaISR ->Fill(ptcluster,etacluster);
1479 }
1480 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
1481 {
1482 fhPtPi0Decay ->Fill(ptcluster);
1483 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
1484 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
1485 }
1486 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
1487 {
1488 fhPtOtherDecay ->Fill(ptcluster);
1489 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
1490 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
1491 }
1492 }
6175da48 1493 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
1494 {
1495 fhPtAntiNeutron ->Fill(ptcluster);
1496 fhPhiAntiNeutron ->Fill(ptcluster,phicluster);
1497 fhEtaAntiNeutron ->Fill(ptcluster,etacluster);
20218aea 1498 if(ph->IsTagged() && fCheckConversion) fhPtAntiNeutronTagged ->Fill(ptcluster);
41121cfe 1499
6175da48 1500 }
1501 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
1502 {
1503 fhPtAntiProton ->Fill(ptcluster);
1504 fhPhiAntiProton ->Fill(ptcluster,phicluster);
1505 fhEtaAntiProton ->Fill(ptcluster,etacluster);
20218aea 1506 if(ph->IsTagged() && fCheckConversion) fhPtAntiProtonTagged ->Fill(ptcluster);
41121cfe 1507
6175da48 1508 }
c8fe2783 1509 else{
1510 fhPtUnknown ->Fill(ptcluster);
1511 fhPhiUnknown ->Fill(ptcluster,phicluster);
1512 fhEtaUnknown ->Fill(ptcluster,etacluster);
20218aea 1513 if(ph->IsTagged() && fCheckConversion) fhPtUnknownTagged ->Fill(ptcluster);
41121cfe 1514
c8fe2783 1515
f8006433 1516 // printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
1517 // ph->GetLabel(),ph->Pt());
1518 // for(Int_t i = 0; i < 20; i++) {
1519 // if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
1520 // }
1521 // printf("\n");
1522
c8fe2783 1523 }
1524
6175da48 1525 //....................................................................
c8fe2783 1526 // Access MC information in stack if requested, check that it exists.
1527 Int_t label =ph->GetLabel();
1528 if(label < 0) {
1529 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
1530 continue;
1531 }
1532
1533 Float_t eprim = 0;
1534 Float_t ptprim = 0;
1535 if(GetReader()->ReadStack()){
1536
1537 if(label >= stack->GetNtrack()) {
f8006433 1538 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
1539 continue ;
c8fe2783 1540 }
1541
1542 primary = stack->Particle(label);
1543 if(!primary){
f8006433 1544 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1545 continue;
c8fe2783 1546 }
1547 eprim = primary->Energy();
1548 ptprim = primary->Pt();
1549
1550 }
1551 else if(GetReader()->ReadAODMCParticles()){
1552 //Check which is the input
1553 if(ph->GetInputFileIndex() == 0){
f8006433 1554 if(!mcparticles0) continue;
1555 if(label >= mcparticles0->GetEntriesFast()) {
1556 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1557 label, mcparticles0->GetEntriesFast());
1558 continue ;
1559 }
1560 //Get the particle
1561 aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
1562
c8fe2783 1563 }
f8006433 1564// else {//Second input
1565// if(!mcparticles1) continue;
1566// if(label >= mcparticles1->GetEntriesFast()) {
1567// if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
1568// label, mcparticles1->GetEntriesFast());
1569// continue ;
1570// }
1571// //Get the particle
1572// aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
1573//
1574// }//second input
c8fe2783 1575
1576 if(!aodprimary){
f8006433 1577 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
1578 continue;
c8fe2783 1579 }
1580
1581 eprim = aodprimary->E();
1582 ptprim = aodprimary->Pt();
1583
1584 }
1585
1586 fh2E ->Fill(ecluster, eprim);
1587 fh2Pt ->Fill(ptcluster, ptprim);
1588 fhDeltaE ->Fill(eprim-ecluster);
1589 fhDeltaPt->Fill(ptprim-ptcluster);
1590 if(eprim > 0) fhRatioE ->Fill(ecluster/eprim);
1591 if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
1592
1593 }//Histograms with MC
1594
1595 }// aod loop
591cc579 1596
1c5acb87 1597}
1598
1599
1600//__________________________________________________________________
1601void AliAnaPhoton::Print(const Option_t * opt) const
1602{
477d6cee 1603 //Print some relevant parameters set for the analysis
1604
1605 if(! opt)
1606 return;
1607
1608 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1609 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 1610
477d6cee 1611 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
1612 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
1613 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1614 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 1615 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1e86c71e 1616 printf("Check Pair Conversion = %d\n",fCheckConversion);
1617 printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD);
1618 printf("Conversion pair mass cut = %f\n",fMassCut);
41121cfe 1619 printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
1620 fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);
4cf55759 1621 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2ac125bf 1622 printf("Number of cells in cluster is > %d \n", fNCellsCut);
477d6cee 1623 printf(" \n") ;
1c5acb87 1624
1625}