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