]> 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
Implement different event selection based on goodness of the reconstructed vertex...
[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.
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 }
f8006433 432
f37fa8d2 433 // else if(input == 1)
434 // calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
c8fe2783 435
f37fa8d2 436 //If too small or big pt, skip it
477d6cee 437 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
4cf55759 438 Double_t tof = calo->GetTOF()*1e9;
c8fe2783 439
440 if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
2ac125bf 441
0ae57829 442 if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) continue;
2ac125bf 443
0ae57829 444 //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(),
445 // mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
c8fe2783 446
f37fa8d2 447 //Check acceptance selection
ff45398a 448 if(IsFiducialCutOn()){
449 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
477d6cee 450 if(! in ) continue ;
451 }
0ae57829 452 //printf("Fiducial cut passed \n");
c8fe2783 453
f37fa8d2 454 //Create AOD for analysis
477d6cee 455 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
c8fe2783 456 Int_t label = calo->GetLabel();
591cc579 457 aodph.SetLabel(label);
f8006433 458 //aodph.SetInputFileIndex(input);
c8fe2783 459
f37fa8d2 460 //printf("Index %d, Id %d\n",icalo, calo->GetID());
461 //Set the indeces of the original caloclusters
477d6cee 462 aodph.SetCaloLabel(calo->GetID(),-1);
463 aodph.SetDetector(fCalorimeter);
464 if(GetDebug() > 1)
ff45398a 465 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 466
f37fa8d2 467 //Check Distance to Bad channel, set bit.
c8fe2783 468 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
477d6cee 469 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
470 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
471 continue ;
472
a3aebfff 473 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
477d6cee 474
f37fa8d2 475 if (distBad > fMinDist3) aodph.SetDistToBad(2) ;
477d6cee 476 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
f37fa8d2 477 else aodph.SetDistToBad(0) ;
af7b3903 478 //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
479
f37fa8d2 480 //Skip matched clusters with tracks
481 if(fRejectTrackMatch && IsTrackMatched(calo)) continue ;
482 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - TrackMatching cut passed \n");
c8fe2783 483
f37fa8d2 484 //Check PID
485 //PID selection or bit setting
477d6cee 486 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
f37fa8d2 487 //Get most probable PID, check PID weights (in MC this option is mandatory)
c8fe2783 488 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
a3aebfff 489 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
f37fa8d2 490 //If primary is not photon, skip it.
477d6cee 491 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
492 }
493 else if(IsCaloPIDOn()){
c8fe2783 494
f37fa8d2 495 //Get most probable PID, 2 options check PID weights
496 //or redo PID, recommended option for EMCal.
477d6cee 497 if(!IsCaloPIDRecalculationOn())
c8fe2783 498 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
477d6cee 499 else
c8fe2783 500 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
477d6cee 501
a3aebfff 502 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
477d6cee 503
f37fa8d2 504 //If cluster does not pass pid, not photon, skip it.
477d6cee 505 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
506
507 }
508 else{
f37fa8d2 509 //Set PID bits for later selection (AliAnaPi0 for example)
510 //GetPDG already called in SetPIDBits.
f2ccb5b8 511 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
a3aebfff 512 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
477d6cee 513 }
514
a3aebfff 515 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
477d6cee 516
f37fa8d2 517 //Play with the MC stack if available
518 //Check origin of the candidates
477d6cee 519 if(IsDataMC()){
c8fe2783 520
521 aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
591cc579 522 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
477d6cee 523 }//Work with stack also
524
c8fe2783 525
f37fa8d2 526 // Check if cluster comes from a conversion in the material in front of the calorimeter
527 // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
c8fe2783 528
529 if(fCheckConversion && nCaloClusters > 1){
530 Bool_t bConverted = kFALSE;
531 Int_t id2 = -1;
1e86c71e 532
f37fa8d2 533 //Check if set previously as converted couple, if so skip its use.
c8fe2783 534 if (indexConverted[icalo]) continue;
1e86c71e 535
c8fe2783 536 for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
f37fa8d2 537 //Check if set previously as converted couple, if so skip its use.
c8fe2783 538 if (indexConverted[jcalo]) continue;
f37fa8d2 539 //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
0ae57829 540 AliVCluster * calo2 = (AliVCluster*) (pl->At(jcalo)); //Get cluster kinematics
c8fe2783 541 Int_t evtIndex2 = 0 ;
542 if (GetMixedEvent()) {
543 evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
f8006433 544
c8fe2783 545 }
f8006433 546
547 if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
8cdc266d 548 calo2->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
f8006433 549 else{
550 Double_t vertex[]={0,0,0};
8cdc266d 551 calo2->GetMomentum(mom2,vertex) ;
f8006433 552 }
553
f37fa8d2 554 //Check only certain regions
c8fe2783 555 Bool_t in2 = kTRUE;
556 if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
557 if(!in2) continue;
558
f37fa8d2 559 //Get mass of pair, if small, take this pair.
560 //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut);
c8fe2783 561 if((mom+mom2).M() < fMassCut){
562 bConverted = kTRUE;
563 id2 = calo2->GetID();
564 indexConverted[jcalo]=kTRUE;
8cdc266d 565 if(GetDebug() > 2)
566 printf("AliAnaPhoton::MakeAnalysisFillAOD(): Pair with mass %f < %2.4f; \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",
567 (mom+mom2).M(),fMassCut,
568 calo->GetID(),calo->E(),GetCaloUtils()->GetModuleNumber(calo), mom.Eta(), mom.Phi(),
569 id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2),mom2.Eta(), mom2.Phi());
c8fe2783 570 break;
571 }
1e86c71e 572
c8fe2783 573 }//Mass loop
1e86c71e 574
c8fe2783 575 if(bConverted){
576 if(fAddConvertedPairsToAOD){
f37fa8d2 577 //Create AOD of pair analysis
c8fe2783 578 TLorentzVector mpair = mom+mom2;
579 AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
580 aodpair.SetLabel(aodph.GetLabel());
f8006433 581 //aodpair.SetInputFileIndex(input);
c8fe2783 582
f37fa8d2 583 //printf("Index %d, Id %d\n",icalo, calo->GetID());
584 //Set the indeces of the original caloclusters
c8fe2783 585 aodpair.SetCaloLabel(calo->GetID(),id2);
586 aodpair.SetDetector(fCalorimeter);
587 aodpair.SetPdg(aodph.GetPdg());
588 aodpair.SetTag(aodph.GetTag());
589
f37fa8d2 590 //Add AOD with pair object to aod branch
c8fe2783 591 AddAODParticle(aodpair);
f37fa8d2 592 //printf("\t \t both added pair\n");
c8fe2783 593 }
594
f37fa8d2 595 //Do not add the current calocluster
c8fe2783 596 continue;
597 }//converted pair
598 }//check conversion
f37fa8d2 599 //printf("\t \t added single cluster %d\n",icalo);
1e86c71e 600
f37fa8d2 601 //Add AOD with photon object to aod branch
477d6cee 602 AddAODParticle(aodph);
603
604 }//loop
605
4a745797 606 delete [] indexConverted;
607
f37fa8d2 608 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
477d6cee 609
1c5acb87 610}
611
612//__________________________________________________________________
613void AliAnaPhoton::MakeAnalysisFillHistograms()
614{
f8006433 615 //Do analysis and fill histograms
616
591cc579 617 // Access MC information in stack if requested, check that it exists.
618 AliStack * stack = 0x0;
619 TParticle * primary = 0x0;
620 TClonesArray * mcparticles0 = 0x0;
f8006433 621 //TClonesArray * mcparticles1 = 0x0;
591cc579 622 AliAODMCParticle * aodprimary = 0x0;
623 if(IsDataMC()){
624
625 if(GetReader()->ReadStack()){
626 stack = GetMCStack() ;
627 if(!stack) {
628 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
629 abort();
630 }
f8006433 631
591cc579 632 }
633 else if(GetReader()->ReadAODMCParticles()){
f8006433 634
591cc579 635 //Get the list of MC particles
636 mcparticles0 = GetReader()->GetAODMCParticles(0);
637 if(!mcparticles0 && GetDebug() > 0) {
638 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
639 }
f8006433 640 // if(GetReader()->GetSecondInputAODTree()){
641 // mcparticles1 = GetReader()->GetAODMCParticles(1);
642 // if(!mcparticles1 && GetDebug() > 0) {
643 // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
644 // }
645 // }
591cc579 646
647 }
648 }// is data and MC
649
123fc3bd 650 //Loop on stored AOD photons
2244659d 651 Double_t v[3] = {0,0,0}; //vertex ;
652 GetReader()->GetVertex(v);
653 fhVertex->Fill(v[0],v[1],v[2]);
654 if(TMath::Abs(v[2]) > GetZvertexCut()) return ;
123fc3bd 655 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2244659d 656 fhNtraNclu->Fill(GetReader()->GetTrackMultiplicity(), naod);
123fc3bd 657 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
658
659 for(Int_t iaod = 0; iaod < naod ; iaod++){
c8fe2783 660 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
661 Int_t pdg = ph->GetPdg();
662
663 if(GetDebug() > 3)
664 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
665
666 //If PID used, fill histos with photons in Calorimeter fCalorimeter
667 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
668 if(ph->GetDetector() != fCalorimeter) continue;
669
670 if(GetDebug() > 2)
671 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
672
673 //Fill photon histograms
674 Float_t ptcluster = ph->Pt();
675 Float_t phicluster = ph->Phi();
676 Float_t etacluster = ph->Eta();
677 Float_t ecluster = ph->E();
678
679 fhPtPhoton ->Fill(ptcluster);
680 fhPhiPhoton ->Fill(ptcluster,phicluster);
681 fhEtaPhoton ->Fill(ptcluster,etacluster);
682
683 //Play with the MC data if available
684 if(IsDataMC()){
685
686 Int_t tag =ph->GetTag();
687
688 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
f8006433 689 {
690 fhPtMCPhoton ->Fill(ptcluster);
691 fhPhiMCPhoton ->Fill(ptcluster,phicluster);
692 fhEtaMCPhoton ->Fill(ptcluster,etacluster);
693
694 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
695 {
696 fhPtConversion ->Fill(ptcluster);
697 fhPhiConversion ->Fill(ptcluster,phicluster);
698 fhEtaConversion ->Fill(ptcluster,etacluster);
699 }
700
701 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
702 fhPtPrompt ->Fill(ptcluster);
703 fhPhiPrompt ->Fill(ptcluster,phicluster);
704 fhEtaPrompt ->Fill(ptcluster,etacluster);
705 }
706 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
707 {
708 fhPtFragmentation ->Fill(ptcluster);
709 fhPhiFragmentation ->Fill(ptcluster,phicluster);
710 fhEtaFragmentation ->Fill(ptcluster,etacluster);
711 }
712 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
713 {
714 fhPtISR ->Fill(ptcluster);
715 fhPhiISR ->Fill(ptcluster,phicluster);
716 fhEtaISR ->Fill(ptcluster,etacluster);
717 }
718 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
719 {
720 fhPtPi0Decay ->Fill(ptcluster);
721 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
722 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
723 }
724 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
725 {
726 fhPtOtherDecay ->Fill(ptcluster);
727 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
728 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
729 }
730 }
c8fe2783 731 else{
732 fhPtUnknown ->Fill(ptcluster);
733 fhPhiUnknown ->Fill(ptcluster,phicluster);
734 fhEtaUnknown ->Fill(ptcluster,etacluster);
735
f8006433 736 // printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
737 // ph->GetLabel(),ph->Pt());
738 // for(Int_t i = 0; i < 20; i++) {
739 // if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
740 // }
741 // printf("\n");
742
c8fe2783 743 }
744
745
746 // Access MC information in stack if requested, check that it exists.
747 Int_t label =ph->GetLabel();
748 if(label < 0) {
749 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
750 continue;
751 }
752
753 Float_t eprim = 0;
754 Float_t ptprim = 0;
755 if(GetReader()->ReadStack()){
756
757 if(label >= stack->GetNtrack()) {
f8006433 758 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
759 continue ;
c8fe2783 760 }
761
762 primary = stack->Particle(label);
763 if(!primary){
f8006433 764 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
765 continue;
c8fe2783 766 }
767 eprim = primary->Energy();
768 ptprim = primary->Pt();
769
770 }
771 else if(GetReader()->ReadAODMCParticles()){
772 //Check which is the input
773 if(ph->GetInputFileIndex() == 0){
f8006433 774 if(!mcparticles0) continue;
775 if(label >= mcparticles0->GetEntriesFast()) {
776 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
777 label, mcparticles0->GetEntriesFast());
778 continue ;
779 }
780 //Get the particle
781 aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
782
c8fe2783 783 }
f8006433 784// else {//Second input
785// if(!mcparticles1) continue;
786// if(label >= mcparticles1->GetEntriesFast()) {
787// if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
788// label, mcparticles1->GetEntriesFast());
789// continue ;
790// }
791// //Get the particle
792// aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
793//
794// }//second input
c8fe2783 795
796 if(!aodprimary){
f8006433 797 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
798 continue;
c8fe2783 799 }
800
801 eprim = aodprimary->E();
802 ptprim = aodprimary->Pt();
803
804 }
805
806 fh2E ->Fill(ecluster, eprim);
807 fh2Pt ->Fill(ptcluster, ptprim);
808 fhDeltaE ->Fill(eprim-ecluster);
809 fhDeltaPt->Fill(ptprim-ptcluster);
810 if(eprim > 0) fhRatioE ->Fill(ecluster/eprim);
811 if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
812
813 }//Histograms with MC
814
815 }// aod loop
591cc579 816
1c5acb87 817}
818
819
820//__________________________________________________________________
821void AliAnaPhoton::Print(const Option_t * opt) const
822{
477d6cee 823 //Print some relevant parameters set for the analysis
824
825 if(! opt)
826 return;
827
828 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
829 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 830
477d6cee 831 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
832 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
833 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
834 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 835 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1e86c71e 836 printf("Check Pair Conversion = %d\n",fCheckConversion);
837 printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD);
838 printf("Conversion pair mass cut = %f\n",fMassCut);
4cf55759 839 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2ac125bf 840 printf("Number of cells in cluster is > %d \n", fNCellsCut);
477d6cee 841 printf(" \n") ;
1c5acb87 842
843}