]> 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
Adding HPTDC parameters that are to be used in MC (hits->digits). The preprocessor...
[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 *
9 * without fee, provided that the above copyright notice appears in all *
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.
22//
23// -- Author: Gustavo Conesa (LNF-INFN)
24//////////////////////////////////////////////////////////////////////////////
25
26
27// --- ROOT system ---
28#include <TH2F.h>
477d6cee 29#include <TClonesArray.h>
30#include <TObjString.h>
31//#include <Riostream.h>
123fc3bd 32#include "TParticle.h"
1c5acb87 33
34// --- Analysis system ---
35#include "AliAnaPhoton.h"
36#include "AliCaloTrackReader.h"
123fc3bd 37#include "AliStack.h"
1c5acb87 38#include "AliCaloPID.h"
6639984f 39#include "AliMCAnalysisUtils.h"
ff45398a 40#include "AliFiducialCut.h"
477d6cee 41#include "AliAODCaloCluster.h"
591cc579 42#include "AliAODMCParticle.h"
1c5acb87 43
44ClassImp(AliAnaPhoton)
45
46//____________________________________________________________________________
47 AliAnaPhoton::AliAnaPhoton() :
48 AliAnaPartCorrBaseClass(), fCalorimeter(""),
a3aebfff 49 fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
1e86c71e 50 fCheckConversion(kFALSE),fAddConvertedPairsToAOD(kFALSE), fMassCut(0),
a3aebfff 51 fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
7175a03a 52 //MC
123fc3bd 53 fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
591cc579 54 fhPtMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0),
1c5acb87 55 fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
56 fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
6639984f 57 fhPtISR(0),fhPhiISR(0),fhEtaISR(0),
1c5acb87 58 fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
59 fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
60 fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
61 fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
62{
63 //default ctor
64
65 //Initialize parameters
66 InitParameters();
67
68}
69
70//____________________________________________________________________________
71AliAnaPhoton::AliAnaPhoton(const AliAnaPhoton & g) :
72 AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter),
1e86c71e 73 fMinDist(g.fMinDist),fMinDist2(g.fMinDist2), fMinDist3(g.fMinDist3),
74 fRejectTrackMatch(g.fRejectTrackMatch),
75 fCheckConversion(g.fCheckConversion),fAddConvertedPairsToAOD(g.fAddConvertedPairsToAOD),
76 fMassCut(g.fMassCut),
77 fhPtPhoton(g.fhPtPhoton),fhPhiPhoton(g.fhPhiPhoton),fhEtaPhoton(g.fhEtaPhoton),
1c5acb87 78 //MC
123fc3bd 79 fhDeltaE(g.fhDeltaE), fhDeltaPt(g.fhDeltaPt),
80 fhRatioE(g.fhRatioE), fhRatioPt(g.fhRatioPt),
81 fh2E(g.fh2E), fh2Pt(g.fh2Pt),
591cc579 82 fhPtMCPhoton(g.fhPtMCPhoton),fhPhiMCPhoton(g.fhPhiMCPhoton),fhEtaMCPhoton(g.fhEtaMCPhoton),
1c5acb87 83 fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
84 fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
6639984f 85 fhPtISR(g.fhPtISR),fhPhiISR(g.fhPhiISR),fhEtaISR(g.fhEtaISR),
1c5acb87 86 fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
87 fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
88 fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
6639984f 89 fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown)
90
1c5acb87 91{
92 // cpy ctor
93
94}
95
96//_________________________________________________________________________
97AliAnaPhoton & AliAnaPhoton::operator = (const AliAnaPhoton & g)
98{
99 // assignment operator
100
101 if(&g == this) return *this;
102
103 fCalorimeter = g.fCalorimeter ;
104 fMinDist = g.fMinDist;
105 fMinDist2 = g.fMinDist2;
106 fMinDist3 = g.fMinDist3;
1e86c71e 107 fMassCut = g.fMassCut;
108
109 fRejectTrackMatch = g.fRejectTrackMatch;
110 fCheckConversion = g.fCheckConversion;
111 fAddConvertedPairsToAOD = g.fAddConvertedPairsToAOD;
112
6639984f 113 fhPtPhoton = g.fhPtPhoton ;
1c5acb87 114 fhPhiPhoton = g.fhPhiPhoton ;
115 fhEtaPhoton = g.fhEtaPhoton ;
123fc3bd 116
117 fhDeltaE = g.fhDeltaE;
118 fhDeltaPt = g.fhDeltaPt;
119 fhRatioE = g.fhRatioE;
120 fhRatioPt = g.fhRatioPt;
121 fh2E = g.fh2E;
122 fh2Pt = g.fh2Pt;
1c5acb87 123
591cc579 124 fhPtMCPhoton = g.fhPtMCPhoton ;
125 fhPhiMCPhoton = g.fhPhiMCPhoton ;
126 fhEtaMCPhoton = g.fhEtaMCPhoton ;
6639984f 127 fhPtPrompt = g.fhPtPrompt;
1c5acb87 128 fhPhiPrompt = g.fhPhiPrompt;
129 fhEtaPrompt = g.fhEtaPrompt;
6639984f 130 fhPtFragmentation = g.fhPtFragmentation;
1c5acb87 131 fhPhiFragmentation = g.fhPhiFragmentation;
132 fhEtaFragmentation = g.fhEtaFragmentation;
6639984f 133 fhPtISR = g.fhPtISR;
134 fhPhiISR = g.fhPhiISR;
135 fhEtaISR = g.fhEtaISR;
136 fhPtPi0Decay = g.fhPtPi0Decay;
1c5acb87 137 fhPhiPi0Decay = g.fhPhiPi0Decay;
138 fhEtaPi0Decay = g.fhEtaPi0Decay;
6639984f 139 fhPtOtherDecay = g.fhPtOtherDecay;
1c5acb87 140 fhPhiOtherDecay = g.fhPhiOtherDecay;
141 fhEtaOtherDecay = g.fhEtaOtherDecay;
6639984f 142 fhPtConversion = g. fhPtConversion;
1c5acb87 143 fhPhiConversion = g.fhPhiConversion;
144 fhEtaConversion = g.fhEtaConversion;
6639984f 145 fhPtUnknown = g.fhPtUnknown;
1c5acb87 146 fhPhiUnknown = g.fhPhiUnknown;
147 fhEtaUnknown = g.fhEtaUnknown;
148
149 return *this;
150
151}
152
153//____________________________________________________________________________
154AliAnaPhoton::~AliAnaPhoton()
155{
156 //dtor
157
158}
159
160
161//________________________________________________________________________
162TList * AliAnaPhoton::GetCreateOutputObjects()
163{
477d6cee 164 // Create histograms to be saved in output file and
165 // store them in outputContainer
166 TList * outputContainer = new TList() ;
167 outputContainer->SetName("PhotonHistos") ;
168
5a2dbc3c 169 Int_t nptbins = GetHistoPtBins();
170 Int_t nphibins = GetHistoPhiBins();
171 Int_t netabins = GetHistoEtaBins();
477d6cee 172 Float_t ptmax = GetHistoPtMax();
173 Float_t phimax = GetHistoPhiMax();
174 Float_t etamax = GetHistoEtaMax();
175 Float_t ptmin = GetHistoPtMin();
176 Float_t phimin = GetHistoPhiMin();
177 Float_t etamin = GetHistoEtaMin();
178
179 //Histograms of highest Photon identified in Event
180 fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
181 fhPtPhoton->SetYTitle("N");
182 fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
183 outputContainer->Add(fhPtPhoton) ;
184
185 fhPhiPhoton = new TH2F
186 ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
187 fhPhiPhoton->SetYTitle("#phi");
188 fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
189 outputContainer->Add(fhPhiPhoton) ;
190
191 fhEtaPhoton = new TH2F
192 ("hEtaPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
193 fhEtaPhoton->SetYTitle("#eta");
194 fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
195 outputContainer->Add(fhEtaPhoton) ;
196
197 if(IsDataMC()){
123fc3bd 198 fhDeltaE = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
199 fhDeltaE->SetXTitle("#Delta E (GeV)");
200 outputContainer->Add(fhDeltaE);
201
202 fhDeltaPt = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
203 fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
204 outputContainer->Add(fhDeltaPt);
205
206 fhRatioE = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
207 fhRatioE->SetXTitle("E_{reco}/E_{gen}");
208 outputContainer->Add(fhRatioE);
477d6cee 209
123fc3bd 210 fhRatioPt = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
211 fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
212 outputContainer->Add(fhRatioPt);
213
214 fh2E = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 215 fh2E->SetXTitle("E_{rec} (GeV)");
216 fh2E->SetYTitle("E_{gen} (GeV)");
123fc3bd 217 outputContainer->Add(fh2E);
218
219 fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 220 fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
221 fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
123fc3bd 222 outputContainer->Add(fh2Pt);
223
591cc579 224 fhPtMCPhoton = new TH1F("hPtMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
225 fhPtMCPhoton->SetYTitle("N");
226 fhPtMCPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
227 outputContainer->Add(fhPtMCPhoton) ;
228
229 fhPhiMCPhoton = new TH2F
230 ("hPhiMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
231 fhPhiMCPhoton->SetYTitle("#phi");
232 fhPhiMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
233 outputContainer->Add(fhPhiMCPhoton) ;
234
235 fhEtaMCPhoton = new TH2F
236 ("hEtaMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
237 fhEtaMCPhoton->SetYTitle("#eta");
238 fhEtaMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
239 outputContainer->Add(fhEtaMCPhoton) ;
240
241 fhPtPrompt = new TH1F("hPtMCPrompt","Number of prompt #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 242 fhPtPrompt->SetYTitle("N");
243 fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
244 outputContainer->Add(fhPtPrompt) ;
245
246 fhPhiPrompt = new TH2F
a3aebfff 247 ("hPhiMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 248 fhPhiPrompt->SetYTitle("#phi");
249 fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
250 outputContainer->Add(fhPhiPrompt) ;
251
252 fhEtaPrompt = new TH2F
a3aebfff 253 ("hEtaMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 254 fhEtaPrompt->SetYTitle("#eta");
255 fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
256 outputContainer->Add(fhEtaPrompt) ;
257
591cc579 258 fhPtFragmentation = new TH1F("hPtMCFragmentation","Number of fragmentation #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 259 fhPtFragmentation->SetYTitle("N");
260 fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
261 outputContainer->Add(fhPtFragmentation) ;
262
263 fhPhiFragmentation = new TH2F
a3aebfff 264 ("hPhiMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 265 fhPhiFragmentation->SetYTitle("#phi");
266 fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
267 outputContainer->Add(fhPhiFragmentation) ;
268
269 fhEtaFragmentation = new TH2F
a3aebfff 270 ("hEtaMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 271 fhEtaFragmentation->SetYTitle("#eta");
272 fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
273 outputContainer->Add(fhEtaFragmentation) ;
274
a3aebfff 275 fhPtISR = new TH1F("hPtMCISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 276 fhPtISR->SetYTitle("N");
277 fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
278 outputContainer->Add(fhPtISR) ;
279
280 fhPhiISR = new TH2F
a3aebfff 281 ("hPhiMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 282 fhPhiISR->SetYTitle("#phi");
283 fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
284 outputContainer->Add(fhPhiISR) ;
285
286 fhEtaISR = new TH2F
a3aebfff 287 ("hEtaMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 288 fhEtaISR->SetYTitle("#eta");
289 fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
290 outputContainer->Add(fhEtaISR) ;
291
591cc579 292 fhPtPi0Decay = new TH1F("hPtMCPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 293 fhPtPi0Decay->SetYTitle("N");
294 fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
295 outputContainer->Add(fhPtPi0Decay) ;
296
297 fhPhiPi0Decay = new TH2F
a3aebfff 298 ("hPhiMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 299 fhPhiPi0Decay->SetYTitle("#phi");
300 fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
301 outputContainer->Add(fhPhiPi0Decay) ;
302
303 fhEtaPi0Decay = new TH2F
a3aebfff 304 ("hEtaMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 305 fhEtaPi0Decay->SetYTitle("#eta");
306 fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
307 outputContainer->Add(fhEtaPi0Decay) ;
308
a3aebfff 309 fhPtOtherDecay = new TH1F("hPtMCOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 310 fhPtOtherDecay->SetYTitle("N");
311 fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
312 outputContainer->Add(fhPtOtherDecay) ;
313
314 fhPhiOtherDecay = new TH2F
a3aebfff 315 ("hPhiMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 316 fhPhiOtherDecay->SetYTitle("#phi");
317 fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
318 outputContainer->Add(fhPhiOtherDecay) ;
319
320 fhEtaOtherDecay = new TH2F
a3aebfff 321 ("hEtaMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 322 fhEtaOtherDecay->SetYTitle("#eta");
323 fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
324 outputContainer->Add(fhEtaOtherDecay) ;
325
a3aebfff 326 fhPtConversion = new TH1F("hPtMCConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 327 fhPtConversion->SetYTitle("N");
328 fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
329 outputContainer->Add(fhPtConversion) ;
330
331 fhPhiConversion = new TH2F
a3aebfff 332 ("hPhiMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 333 fhPhiConversion->SetYTitle("#phi");
334 fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
335 outputContainer->Add(fhPhiConversion) ;
336
337 fhEtaConversion = new TH2F
a3aebfff 338 ("hEtaMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 339 fhEtaConversion->SetYTitle("#eta");
340 fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
341 outputContainer->Add(fhEtaConversion) ;
342
a3aebfff 343 fhPtUnknown = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 344 fhPtUnknown->SetYTitle("N");
345 fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
346 outputContainer->Add(fhPtUnknown) ;
347
348 fhPhiUnknown = new TH2F
a3aebfff 349 ("hPhiMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 350 fhPhiUnknown->SetYTitle("#phi");
351 fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
352 outputContainer->Add(fhPhiUnknown) ;
353
354 fhEtaUnknown = new TH2F
a3aebfff 355 ("hEtaMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 356 fhEtaUnknown->SetYTitle("#eta");
357 fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
358 outputContainer->Add(fhEtaUnknown) ;
a3aebfff 359
477d6cee 360 }//Histos with MC
361
362 //Save parameters used for analysis
363 TString parList ; //this will be list of parameters used for this analysis.
364 char onePar[255] ;
365
366 sprintf(onePar,"--- AliAnaPhoton ---\n") ;
367 parList+=onePar ;
368 sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
369 parList+=onePar ;
370 sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
371 parList+=onePar ;
372 sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
373 parList+=onePar ;
374 sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
375 parList+=onePar ;
a3aebfff 376 sprintf(onePar,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
377 parList+=onePar ;
477d6cee 378
379 //Get parameters set in base class.
380 parList += GetBaseParametersList() ;
381
382 //Get parameters set in PID class.
383 parList += GetCaloPID()->GetPIDParametersList() ;
384
ff45398a 385 //Get parameters set in FiducialCut class (not available yet)
477d6cee 386 //parlist += GetFidCut()->GetFidCutParametersList()
387
388 TObjString *oString= new TObjString(parList) ;
389 outputContainer->Add(oString);
390
391 return outputContainer ;
392
1c5acb87 393}
394
6639984f 395//____________________________________________________________________________
396void AliAnaPhoton::Init()
397{
398
399 //Init
400 //Do some checks
1e86c71e 401 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
591cc579 402 printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 403 abort();
404 }
1e86c71e 405 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
591cc579 406 printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 407 abort();
408 }
409
410}
411
412
1c5acb87 413//____________________________________________________________________________
414void AliAnaPhoton::InitParameters()
415{
416
417 //Initialize the parameters of the analysis.
a3aebfff 418 AddToHistogramsName("AnaPhoton_");
419
1c5acb87 420 fCalorimeter = "PHOS" ;
421 fMinDist = 2.;
422 fMinDist2 = 4.;
423 fMinDist3 = 5.;
1e86c71e 424 fMassCut = 0.03; //30 MeV
425
426 fRejectTrackMatch = kTRUE ;
427 fCheckConversion = kFALSE;
428 fAddConvertedPairsToAOD = kFALSE;
429
1c5acb87 430}
431
432//__________________________________________________________________
433void AliAnaPhoton::MakeAnalysisFillAOD()
434{
477d6cee 435 //Do analysis and fill aods
436 //Search for photons in fCalorimeter
591cc579 437 TObjArray * pl = new TObjArray();
477d6cee 438
439 //Get vertex for photon momentum calculation
233e0df8 440 Double_t vertex[] = {0,0,0} ; //vertex
441 Double_t vertex2[] = {0,0,0} ; //vertex from second input aod
edd59991 442 if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
233e0df8 443 {
444 GetReader()->GetVertex(vertex);
445 if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
446 }
edd59991 447 //printf("Vertex 0: %f,%f,%f\n",vertex[0],vertex[1],vertex[2]);
448 //printf("Vertex 1: %f,%f,%f\n",vertex2[0],vertex2[1],vertex2[2]);
477d6cee 449 //Select the Calorimeter of the photon
450 if(fCalorimeter == "PHOS")
451 pl = GetAODPHOS();
452 else if (fCalorimeter == "EMCAL")
453 pl = GetAODEMCAL();
454
455 //Fill AODCaloClusters and AODParticle with PHOS aods
1e86c71e 456 TLorentzVector mom, mom2 ;
457 Int_t nCaloClusters = pl->GetEntriesFast();
458 Bool_t * indexConverted = new Bool_t[nCaloClusters];
459 for (Int_t i = 0; i < nCaloClusters; i++) indexConverted[i] = kFALSE;
460
461 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
462
463 AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));
233e0df8 464
ff45398a 465 //Cluster selection, not charged, with photon id and in fiducial cut
233e0df8 466
467 //Input from second AOD?
468 Int_t input = 0;
469 if (fCalorimeter == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= icalo) input = 1 ;
470 else if(fCalorimeter == "PHOS" && GetReader()->GetAODPHOSNormalInputEntries() <= icalo) input = 1;
471
477d6cee 472 //Get Momentum vector,
233e0df8 473 if (input == 0) calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
474 else if(input == 1) calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
475
477d6cee 476 //If too small or big pt, skip it
477 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
591cc579 478
479 //printf("AliAnaPhoton::Current Event %d; Current File Name : %s, E %f, pT %f, Ecl %f\n",GetReader()->GetEventNumber(),(GetReader()->GetCurrentFileName()).Data(), mom.E(), mom.Pt(),calo->E());
480
477d6cee 481 //Check acceptance selection
ff45398a 482 if(IsFiducialCutOn()){
483 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
477d6cee 484 if(! in ) continue ;
485 }
591cc579 486
591cc579 487 //Create AOD for analysis
477d6cee 488 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
591cc579 489 Int_t label = calo->GetLabel(0);
591cc579 490 aodph.SetLabel(label);
233e0df8 491 aodph.SetInputFileIndex(input);
7cd4e982 492
477d6cee 493 //printf("Index %d, Id %d\n",icalo, calo->GetID());
494 //Set the indeces of the original caloclusters
495 aodph.SetCaloLabel(calo->GetID(),-1);
496 aodph.SetDetector(fCalorimeter);
497 if(GetDebug() > 1)
ff45398a 498 printf("AliAnaPhoton::MakeAnalysisFillAOD() - Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());
477d6cee 499
500 //Check Distance to Bad channel, set bit.
501 Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
502 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
503 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
504 continue ;
505
a3aebfff 506 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
477d6cee 507
508 if(distBad > fMinDist3) aodph.SetDistToBad(2) ;
509 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
510 else aodph.SetDistToBad(0) ;
511
ece3ffc2 512 //Skip matched clusters with tracks
513 if(fRejectTrackMatch && calo->GetNTracksMatched() > 0) continue ;
514
477d6cee 515 //Check PID
516 //PID selection or bit setting
517 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
518 //Get most probable PID, check PID weights (in MC this option is mandatory)
519 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
a3aebfff 520 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
477d6cee 521 //If primary is not photon, skip it.
522 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
523 }
524 else if(IsCaloPIDOn()){
ece3ffc2 525
477d6cee 526 //Get most probable PID, 2 options check PID weights
527 //or redo PID, recommended option for EMCal.
528 if(!IsCaloPIDRecalculationOn())
529 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
530 else
531 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
532
a3aebfff 533 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
477d6cee 534
535 //If cluster does not pass pid, not photon, skip it.
536 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
537
538 }
539 else{
540 //Set PID bits for later selection (AliAnaPi0 for example)
541 //GetPDG already called in SetPIDBits.
542 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);
a3aebfff 543 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
477d6cee 544 }
545
a3aebfff 546 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
477d6cee 547
548 //Play with the MC stack if available
549 //Check origin of the candidates
550 if(IsDataMC()){
591cc579 551
7cd4e982 552 aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetReader(), aodph.GetInputFileIndex()));
591cc579 553 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
477d6cee 554 }//Work with stack also
555
1e86c71e 556
557 // Check if cluster comes from a conversion in the material in front of the calorimeter
558 // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
559
560 if(fCheckConversion && nCaloClusters > 1){
561 Bool_t bConverted = kFALSE;
562 Int_t id2 = -1;
563
564 //Check if set previously as converted couple, if so skip its use.
565 if (indexConverted[icalo]) continue;
566
567 for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
568 //Check if set previously as converted couple, if so skip its use.
569 if (indexConverted[jcalo]) continue;
570 //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
571 AliAODCaloCluster * calo2 = (AliAODCaloCluster*) (pl->At(jcalo)); //Get cluster kinematics
572 calo2->GetMomentum(mom2,vertex);
573 //Check only certain regions
574 Bool_t in2 = kTRUE;
575 if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
576 if(!in2) continue;
577
578 //Get mass of pair, if small, take this pair.
579 //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut);
580 if((mom+mom2).M() < fMassCut){
581 bConverted = kTRUE;
582 id2 = calo2->GetID();
583 indexConverted[jcalo]=kTRUE;
584 break;
585 }
586
587 }//Mass loop
588
589 if(bConverted){
590 if(fAddConvertedPairsToAOD){
591 //Create AOD of pair analysis
592 TLorentzVector mpair = mom+mom2;
593 AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
594 aodpair.SetLabel(aodph.GetLabel());
595 aodpair.SetInputFileIndex(input);
596
597 //printf("Index %d, Id %d\n",icalo, calo->GetID());
598 //Set the indeces of the original caloclusters
599 aodpair.SetCaloLabel(calo->GetID(),id2);
600 aodpair.SetDetector(fCalorimeter);
601 aodpair.SetPdg(aodph.GetPdg());
602 aodpair.SetTag(aodph.GetTag());
603
604 //Add AOD with pair object to aod branch
605 AddAODParticle(aodpair);
606 //printf("\t \t both added pair\n");
607 }
608
609 //Do not add the current calocluster
610 continue;
611 }//converted pair
612 }//check conversion
613 //printf("\t \t added single cluster %d\n",icalo);
614
477d6cee 615 //Add AOD with photon object to aod branch
616 AddAODParticle(aodph);
617
618 }//loop
619
a3aebfff 620 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs \n");
477d6cee 621
1c5acb87 622}
623
624//__________________________________________________________________
625void AliAnaPhoton::MakeAnalysisFillHistograms()
626{
123fc3bd 627 //Do analysis and fill histograms
628
591cc579 629 // Access MC information in stack if requested, check that it exists.
630 AliStack * stack = 0x0;
631 TParticle * primary = 0x0;
632 TClonesArray * mcparticles0 = 0x0;
633 TClonesArray * mcparticles1 = 0x0;
634 AliAODMCParticle * aodprimary = 0x0;
635 if(IsDataMC()){
636
637 if(GetReader()->ReadStack()){
638 stack = GetMCStack() ;
639 if(!stack) {
640 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
641 abort();
642 }
643
644 }
645 else if(GetReader()->ReadAODMCParticles()){
646
647 //Get the list of MC particles
648 mcparticles0 = GetReader()->GetAODMCParticles(0);
649 if(!mcparticles0 && GetDebug() > 0) {
650 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
651 }
652 if(GetReader()->GetSecondInputAODTree()){
653 mcparticles1 = GetReader()->GetAODMCParticles(1);
654 if(!mcparticles1 && GetDebug() > 0) {
655 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
656 }
657 }
658
659 }
660 }// is data and MC
661
123fc3bd 662 //Loop on stored AOD photons
663 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
664 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
665
666 for(Int_t iaod = 0; iaod < naod ; iaod++){
477d6cee 667 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
668 Int_t pdg = ph->GetPdg();
669
670 if(GetDebug() > 3)
a3aebfff 671 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
477d6cee 672
673 //If PID used, fill histos with photons in Calorimeter fCalorimeter
123fc3bd 674 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
477d6cee 675 if(ph->GetDetector() != fCalorimeter) continue;
676
591cc579 677 if(GetDebug() > 2)
a3aebfff 678 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
477d6cee 679
680 //Fill photon histograms
123fc3bd 681 Float_t ptcluster = ph->Pt();
477d6cee 682 Float_t phicluster = ph->Phi();
683 Float_t etacluster = ph->Eta();
123fc3bd 684 Float_t ecluster = ph->E();
685
686 fhPtPhoton ->Fill(ptcluster);
477d6cee 687 fhPhiPhoton ->Fill(ptcluster,phicluster);
688 fhEtaPhoton ->Fill(ptcluster,etacluster);
123fc3bd 689
591cc579 690 //Play with the MC data if available
477d6cee 691 if(IsDataMC()){
123fc3bd 692
477d6cee 693 Int_t tag =ph->GetTag();
591cc579 694
695 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
696 {
697 fhPtMCPhoton ->Fill(ptcluster);
698 fhPhiMCPhoton ->Fill(ptcluster,phicluster);
699 fhEtaMCPhoton ->Fill(ptcluster,etacluster);
700
701 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
702 {
703 fhPtConversion ->Fill(ptcluster);
704 fhPhiConversion ->Fill(ptcluster,phicluster);
705 fhEtaConversion ->Fill(ptcluster,etacluster);
706 }
707
708 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
709 fhPtPrompt ->Fill(ptcluster);
710 fhPhiPrompt ->Fill(ptcluster,phicluster);
711 fhEtaPrompt ->Fill(ptcluster,etacluster);
712 }
713 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
714 {
715 fhPtFragmentation ->Fill(ptcluster);
716 fhPhiFragmentation ->Fill(ptcluster,phicluster);
717 fhEtaFragmentation ->Fill(ptcluster,etacluster);
718 }
719 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
720 {
721 fhPtISR ->Fill(ptcluster);
722 fhPhiISR ->Fill(ptcluster,phicluster);
723 fhEtaISR ->Fill(ptcluster,etacluster);
724 }
725 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
726 {
727 fhPtPi0Decay ->Fill(ptcluster);
728 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
729 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
730 }
731 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
732 {
733 fhPtOtherDecay ->Fill(ptcluster);
734 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
735 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
736 }
737 }
477d6cee 738 else{
591cc579 739 fhPtUnknown ->Fill(ptcluster);
740 fhPhiUnknown ->Fill(ptcluster,phicluster);
741 fhEtaUnknown ->Fill(ptcluster,etacluster);
742
743// printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
744// ph->GetLabel(),ph->Pt());
745// for(Int_t i = 0; i < 20; i++) {
746// if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
747// }
748// printf("\n");
749
477d6cee 750 }
591cc579 751
123fc3bd 752
591cc579 753 // Access MC information in stack if requested, check that it exists.
123fc3bd 754 Int_t label =ph->GetLabel();
755 if(label < 0) {
8d4bd2e1 756 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
123fc3bd 757 continue;
758 }
759
591cc579 760 Float_t eprim = 0;
761 Float_t ptprim = 0;
762 if(GetReader()->ReadStack()){
763
764 if(label >= stack->GetNtrack()) {
8d4bd2e1 765 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
591cc579 766 continue ;
767 }
768
769 primary = stack->Particle(label);
770 if(!primary){
771 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
772 continue;
773 }
774 eprim = primary->Energy();
775 ptprim = primary->Pt();
123fc3bd 776
123fc3bd 777 }
591cc579 778 else if(GetReader()->ReadAODMCParticles()){
779 //Check which is the input
7cd4e982 780 if(ph->GetInputFileIndex() == 0){
591cc579 781 if(!mcparticles0) continue;
782 if(label >= mcparticles0->GetEntriesFast()) {
783 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
784 label, mcparticles0->GetEntriesFast());
785 continue ;
786 }
787 //Get the particle
788 aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
789
790 }
791 else {//Second input
792 if(!mcparticles1) continue;
591cc579 793 if(label >= mcparticles1->GetEntriesFast()) {
794 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
795 label, mcparticles1->GetEntriesFast());
796 continue ;
797 }
798 //Get the particle
799 aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
800
801 }//second input
802
803 if(!aodprimary){
804 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
805 continue;
806 }
807
808 eprim = aodprimary->E();
809 ptprim = aodprimary->Pt();
123fc3bd 810
591cc579 811 }
812
41e886c8 813 fh2E ->Fill(ecluster, eprim);
814 fh2Pt ->Fill(ptcluster, ptprim);
591cc579 815 fhDeltaE ->Fill(eprim-ecluster);
816 fhDeltaPt->Fill(ptprim-ptcluster);
817 if(eprim > 0) fhRatioE ->Fill(ecluster/eprim);
818 if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
123fc3bd 819
477d6cee 820 }//Histograms with MC
821
822 }// aod loop
823
1c5acb87 824}
825
826
827//__________________________________________________________________
828void AliAnaPhoton::Print(const Option_t * opt) const
829{
477d6cee 830 //Print some relevant parameters set for the analysis
831
832 if(! opt)
833 return;
834
835 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
836 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 837
477d6cee 838 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
839 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
840 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
841 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 842 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1e86c71e 843 printf("Check Pair Conversion = %d\n",fCheckConversion);
844 printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD);
845 printf("Conversion pair mass cut = %f\n",fMassCut);
477d6cee 846 printf(" \n") ;
1c5acb87 847
848}