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