]> 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
Introduce smaller AliTRDCalDCSFEEv2 object (Frederick)
[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
2244659d 133 fhVertex = new TH3D ("Vertex","vertex position", 100,-50.,50., 100,-50.,50., 100,-50.,50.);
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) ;
477d6cee 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){
548 calo->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
549 else{
550 Double_t vertex[]={0,0,0};
551 calo->GetMomentum(mom2,vertex) ;
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;
565 break;
566 }
1e86c71e 567
c8fe2783 568 }//Mass loop
1e86c71e 569
c8fe2783 570 if(bConverted){
571 if(fAddConvertedPairsToAOD){
f37fa8d2 572 //Create AOD of pair analysis
c8fe2783 573 TLorentzVector mpair = mom+mom2;
574 AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
575 aodpair.SetLabel(aodph.GetLabel());
f8006433 576 //aodpair.SetInputFileIndex(input);
c8fe2783 577
f37fa8d2 578 //printf("Index %d, Id %d\n",icalo, calo->GetID());
579 //Set the indeces of the original caloclusters
c8fe2783 580 aodpair.SetCaloLabel(calo->GetID(),id2);
581 aodpair.SetDetector(fCalorimeter);
582 aodpair.SetPdg(aodph.GetPdg());
583 aodpair.SetTag(aodph.GetTag());
584
f37fa8d2 585 //Add AOD with pair object to aod branch
c8fe2783 586 AddAODParticle(aodpair);
f37fa8d2 587 //printf("\t \t both added pair\n");
c8fe2783 588 }
589
f37fa8d2 590 //Do not add the current calocluster
c8fe2783 591 continue;
592 }//converted pair
593 }//check conversion
f37fa8d2 594 //printf("\t \t added single cluster %d\n",icalo);
1e86c71e 595
f37fa8d2 596 //Add AOD with photon object to aod branch
477d6cee 597 AddAODParticle(aodph);
598
599 }//loop
600
4a745797 601 delete [] indexConverted;
602
f37fa8d2 603 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
477d6cee 604
1c5acb87 605}
606
607//__________________________________________________________________
608void AliAnaPhoton::MakeAnalysisFillHistograms()
609{
f8006433 610 //Do analysis and fill histograms
611
591cc579 612 // Access MC information in stack if requested, check that it exists.
613 AliStack * stack = 0x0;
614 TParticle * primary = 0x0;
615 TClonesArray * mcparticles0 = 0x0;
f8006433 616 //TClonesArray * mcparticles1 = 0x0;
591cc579 617 AliAODMCParticle * aodprimary = 0x0;
618 if(IsDataMC()){
619
620 if(GetReader()->ReadStack()){
621 stack = GetMCStack() ;
622 if(!stack) {
623 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
624 abort();
625 }
f8006433 626
591cc579 627 }
628 else if(GetReader()->ReadAODMCParticles()){
f8006433 629
591cc579 630 //Get the list of MC particles
631 mcparticles0 = GetReader()->GetAODMCParticles(0);
632 if(!mcparticles0 && GetDebug() > 0) {
633 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
634 }
f8006433 635 // if(GetReader()->GetSecondInputAODTree()){
636 // mcparticles1 = GetReader()->GetAODMCParticles(1);
637 // if(!mcparticles1 && GetDebug() > 0) {
638 // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
639 // }
640 // }
591cc579 641
642 }
643 }// is data and MC
644
123fc3bd 645 //Loop on stored AOD photons
2244659d 646 Double_t v[3] = {0,0,0}; //vertex ;
647 GetReader()->GetVertex(v);
648 fhVertex->Fill(v[0],v[1],v[2]);
649 if(TMath::Abs(v[2]) > GetZvertexCut()) return ;
123fc3bd 650 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2244659d 651 fhNtraNclu->Fill(GetReader()->GetTrackMultiplicity(), naod);
123fc3bd 652 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
653
654 for(Int_t iaod = 0; iaod < naod ; iaod++){
c8fe2783 655 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
656 Int_t pdg = ph->GetPdg();
657
658 if(GetDebug() > 3)
659 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
660
661 //If PID used, fill histos with photons in Calorimeter fCalorimeter
662 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
663 if(ph->GetDetector() != fCalorimeter) continue;
664
665 if(GetDebug() > 2)
666 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
667
668 //Fill photon histograms
669 Float_t ptcluster = ph->Pt();
670 Float_t phicluster = ph->Phi();
671 Float_t etacluster = ph->Eta();
672 Float_t ecluster = ph->E();
673
674 fhPtPhoton ->Fill(ptcluster);
675 fhPhiPhoton ->Fill(ptcluster,phicluster);
676 fhEtaPhoton ->Fill(ptcluster,etacluster);
677
678 //Play with the MC data if available
679 if(IsDataMC()){
680
681 Int_t tag =ph->GetTag();
682
683 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
f8006433 684 {
685 fhPtMCPhoton ->Fill(ptcluster);
686 fhPhiMCPhoton ->Fill(ptcluster,phicluster);
687 fhEtaMCPhoton ->Fill(ptcluster,etacluster);
688
689 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
690 {
691 fhPtConversion ->Fill(ptcluster);
692 fhPhiConversion ->Fill(ptcluster,phicluster);
693 fhEtaConversion ->Fill(ptcluster,etacluster);
694 }
695
696 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
697 fhPtPrompt ->Fill(ptcluster);
698 fhPhiPrompt ->Fill(ptcluster,phicluster);
699 fhEtaPrompt ->Fill(ptcluster,etacluster);
700 }
701 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
702 {
703 fhPtFragmentation ->Fill(ptcluster);
704 fhPhiFragmentation ->Fill(ptcluster,phicluster);
705 fhEtaFragmentation ->Fill(ptcluster,etacluster);
706 }
707 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
708 {
709 fhPtISR ->Fill(ptcluster);
710 fhPhiISR ->Fill(ptcluster,phicluster);
711 fhEtaISR ->Fill(ptcluster,etacluster);
712 }
713 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
714 {
715 fhPtPi0Decay ->Fill(ptcluster);
716 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
717 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
718 }
719 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
720 {
721 fhPtOtherDecay ->Fill(ptcluster);
722 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
723 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
724 }
725 }
c8fe2783 726 else{
727 fhPtUnknown ->Fill(ptcluster);
728 fhPhiUnknown ->Fill(ptcluster,phicluster);
729 fhEtaUnknown ->Fill(ptcluster,etacluster);
730
f8006433 731 // printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
732 // ph->GetLabel(),ph->Pt());
733 // for(Int_t i = 0; i < 20; i++) {
734 // if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
735 // }
736 // printf("\n");
737
c8fe2783 738 }
739
740
741 // Access MC information in stack if requested, check that it exists.
742 Int_t label =ph->GetLabel();
743 if(label < 0) {
744 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
745 continue;
746 }
747
748 Float_t eprim = 0;
749 Float_t ptprim = 0;
750 if(GetReader()->ReadStack()){
751
752 if(label >= stack->GetNtrack()) {
f8006433 753 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
754 continue ;
c8fe2783 755 }
756
757 primary = stack->Particle(label);
758 if(!primary){
f8006433 759 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
760 continue;
c8fe2783 761 }
762 eprim = primary->Energy();
763 ptprim = primary->Pt();
764
765 }
766 else if(GetReader()->ReadAODMCParticles()){
767 //Check which is the input
768 if(ph->GetInputFileIndex() == 0){
f8006433 769 if(!mcparticles0) continue;
770 if(label >= mcparticles0->GetEntriesFast()) {
771 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
772 label, mcparticles0->GetEntriesFast());
773 continue ;
774 }
775 //Get the particle
776 aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
777
c8fe2783 778 }
f8006433 779// else {//Second input
780// if(!mcparticles1) continue;
781// if(label >= mcparticles1->GetEntriesFast()) {
782// if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
783// label, mcparticles1->GetEntriesFast());
784// continue ;
785// }
786// //Get the particle
787// aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
788//
789// }//second input
c8fe2783 790
791 if(!aodprimary){
f8006433 792 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
793 continue;
c8fe2783 794 }
795
796 eprim = aodprimary->E();
797 ptprim = aodprimary->Pt();
798
799 }
800
801 fh2E ->Fill(ecluster, eprim);
802 fh2Pt ->Fill(ptcluster, ptprim);
803 fhDeltaE ->Fill(eprim-ecluster);
804 fhDeltaPt->Fill(ptprim-ptcluster);
805 if(eprim > 0) fhRatioE ->Fill(ecluster/eprim);
806 if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
807
808 }//Histograms with MC
809
810 }// aod loop
591cc579 811
1c5acb87 812}
813
814
815//__________________________________________________________________
816void AliAnaPhoton::Print(const Option_t * opt) const
817{
477d6cee 818 //Print some relevant parameters set for the analysis
819
820 if(! opt)
821 return;
822
823 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
824 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 825
477d6cee 826 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
827 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
828 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
829 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 830 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1e86c71e 831 printf("Check Pair Conversion = %d\n",fCheckConversion);
832 printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD);
833 printf("Conversion pair mass cut = %f\n",fMassCut);
4cf55759 834 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2ac125bf 835 printf("Number of cells in cluster is > %d \n", fNCellsCut);
477d6cee 836 printf(" \n") ;
1c5acb87 837
838}