]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaPi0.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliCalorimeterUtils: Fix to be able to use PHOS bad map and geometry matrices
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPi0.cxx
CommitLineData
1c5acb87 1/**************************************************************************
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 is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15/* $Id: $ */
16
17//_________________________________________________________________________
18// Class to collect two-photon invariant mass distributions for
6175da48 19// extracting raw pi0 yield.
20// Input is produced by AliAnaPhoton (or any other analysis producing output AliAODPWG4Particles),
21// it will do nothing if executed alone
1c5acb87 22//
23//-- Author: Dmitri Peressounko (RRC "KI")
24//-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
25//-- and Gustavo Conesa (INFN-Frascati)
26//_________________________________________________________________________
27
28
29// --- ROOT system ---
30#include "TH3.h"
0333ede6 31#include "TH2F.h"
1c5acb87 32//#include "Riostream.h"
6639984f 33#include "TCanvas.h"
34#include "TPad.h"
35#include "TROOT.h"
477d6cee 36#include "TClonesArray.h"
0c1383b5 37#include "TObjString.h"
6175da48 38#include "TDatabasePDG.h"
1c5acb87 39
40//---- AliRoot system ----
41#include "AliAnaPi0.h"
42#include "AliCaloTrackReader.h"
43#include "AliCaloPID.h"
6639984f 44#include "AliStack.h"
ff45398a 45#include "AliFiducialCut.h"
477d6cee 46#include "TParticle.h"
477d6cee 47#include "AliVEvent.h"
6921fa00 48#include "AliESDCaloCluster.h"
49#include "AliESDEvent.h"
50#include "AliAODEvent.h"
50f39b97 51#include "AliNeutralMesonSelection.h"
c8fe2783 52#include "AliMixedEvent.h"
6175da48 53#include "AliAODMCParticle.h"
591cc579 54
c5693f62 55// --- Detectors ---
56#include "AliPHOSGeoUtils.h"
57#include "AliEMCALGeometry.h"
58
1c5acb87 59ClassImp(AliAnaPi0)
60
61//________________________________________________________________________________________________________________________________________________
62AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
0333ede6 63fDoOwnMix(kFALSE), fEventsList(0x0),
64fCalorimeter(""), fNModules(12),
65fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE), fAngleCut(0), fAngleMaxCut(7.),
66fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
67fNPtCuts(0), fNAsymCuts(0), fNCellNCuts(0), fNPIDBits(0),
68fMakeInvPtPlots(kFALSE), fSameSM(kFALSE), fFillSMCombinations(kFALSE), fCheckConversion(kFALSE),
69fUseTrackMultBins(kFALSE), fUsePhotonMultBins(kFALSE), fUseAverCellEBins(kFALSE), fUseAverClusterEBins(kFALSE),
70fUseAverClusterEDenBins(0), fFillBadDistHisto(kFALSE),
71fhAverTotECluster(0), fhAverTotECell(0), fhAverTotECellvsCluster(0),
72fhEDensityCluster(0), fhEDensityCell(0), fhEDensityCellvsCluster(0),
73//Histograms
74fhReMod(0x0), fhReSameSideEMCALMod(0x0), fhReSameSectorEMCALMod(0x0), fhReDiffPHOSMod(0x0),
75fhMiMod(0x0), fhMiSameSideEMCALMod(0x0), fhMiSameSectorEMCALMod(0x0), fhMiDiffPHOSMod(0x0),
76fhReConv(0x0), fhMiConv(0x0), fhReConv2(0x0), fhMiConv2(0x0),
77fhRe1(0x0), fhMi1(0x0), fhRe2(0x0), fhMi2(0x0),
78fhRe3(0x0), fhMi3(0x0), fhReInvPt1(0x0), fhMiInvPt1(0x0),
79fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
80fhRePtNCellAsymCuts(0x0), fhMiPtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM(),
81fhRePIDBits(0x0), fhRePtMult(0x0), fhReSS(),
82fhRePtAsym(0x0), fhRePtAsymPi0(0x0), fhRePtAsymEta(0x0),
83fhEvents(0x0), fhCentrality(0x0), fhCentralityNoPair(0x0),
84fhEventPlaneAngle(0x0), fhEventPlaneResolution(0x0),
85fhRealOpeningAngle(0x0), fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0), fhMixedCosOpeningAngle(0x0),
86// MC histograms
87fhPrimPi0Pt(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0),
88fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0), fhPrimPi0OpeningAngle(0x0), fhPrimPi0CosOpeningAngle(0x0),
89fhPrimEtaPt(0x0), fhPrimEtaAccPt(0x0), fhPrimEtaY(0x0), fhPrimEtaAccY(0x0),
90fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0), fhPrimPi0PtOrigin(0x0), fhPrimEtaPtOrigin(0x0),
91fhMCOrgMass(), fhMCOrgAsym(), fhMCOrgDeltaEta(), fhMCOrgDeltaPhi(),
92fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(),
93fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec(),
99b8e903 94fhMCPi0PtOrigin(0x0), fhMCEtaPtOrigin(0x0),
95fhReMCFromConversion(0), fhReMCFromNotConversion(0), fhReMCFromMixConversion(0)
1c5acb87 96{
97//Default Ctor
98 InitParameters();
99
100}
1c5acb87 101
102//________________________________________________________________________________________________________________________________________________
103AliAnaPi0::~AliAnaPi0() {
477d6cee 104 // Remove event containers
7e7694bb 105
106 if(fDoOwnMix && fEventsList){
72542aba 107 for(Int_t ic=0; ic<GetNCentrBin(); ic++){
5025c139 108 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
109 for(Int_t irp=0; irp<GetNRPBin(); irp++){
110 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->Delete() ;
111 delete fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] ;
7e7694bb 112 }
477d6cee 113 }
114 }
115 delete[] fEventsList;
477d6cee 116 }
591cc579 117
1c5acb87 118}
119
120//________________________________________________________________________________________________________________________________________________
121void AliAnaPi0::InitParameters()
122{
0333ede6 123 //Init parameters when first called the analysis
124 //Set default parameters
a3aebfff 125 SetInputAODName("PWG4Particle");
126
127 AddToHistogramsName("AnaPi0_");
6921fa00 128 fNModules = 12; // set maximum to maximum number of EMCAL modules
0333ede6 129
477d6cee 130 fCalorimeter = "PHOS";
50f39b97 131 fUseAngleCut = kFALSE;
6175da48 132 fUseAngleEDepCut = kFALSE;
133 fAngleCut = 0.;
134 fAngleMaxCut = TMath::Pi();
0333ede6 135
5ae09196 136 fMultiCutAna = kFALSE;
137
20218aea 138 fNPtCuts = 1;
af7b3903 139 fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
140 for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
5ae09196 141
20218aea 142 fNAsymCuts = 2;
143 fAsymCuts[0] = 1.; fAsymCuts[1] = 0.7; //fAsymCuts[2] = 0.6; // fAsymCuts[3] = 0.1;
af7b3903 144 for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
0333ede6 145
20218aea 146 fNCellNCuts = 1;
af7b3903 147 fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
021c573a 148 for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
0333ede6 149
20218aea 150 fNPIDBits = 1;
af7b3903 151 fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
021c573a 152 for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
0333ede6 153
1c5acb87 154}
1c5acb87 155
0c1383b5 156
157//________________________________________________________________________________________________________________________________________________
158TObjString * AliAnaPi0::GetAnalysisCuts()
159{
af7b3903 160 //Save parameters used for analysis
161 TString parList ; //this will be list of parameters used for this analysis.
162 const Int_t buffersize = 255;
163 char onePar[buffersize] ;
164 snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
165 parList+=onePar ;
72542aba 166 snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",GetNCentrBin()) ;
af7b3903 167 parList+=onePar ;
168 snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
169 parList+=onePar ;
170 snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
171 parList+=onePar ;
72542aba 172 snprintf(onePar,buffersize,"Depth of event buffer: %d \n",GetNMaxEvMix()) ;
af7b3903 173 parList+=onePar ;
20218aea 174 snprintf(onePar,buffersize,"Pair in same Module: %d ; Fill Different SM histos %d; CheckConversions %d; TrackMult as centrality: %d; PhotonMult as centrality: %d; cluster E as centrality: %d; cell as centrality: %d; Fill InvPt histos %d\n",
175 fSameSM, fFillSMCombinations, fCheckConversion, fUseTrackMultBins, fUsePhotonMultBins, fUseAverClusterEBins, fUseAverCellEBins, fMakeInvPtPlots) ;
6175da48 176 parList+=onePar ;
177 snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f,\n",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
af7b3903 178 parList+=onePar ;
179 snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
180 for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
181 parList+=onePar ;
182 snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
183 for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
184 parList+=onePar ;
185 snprintf(onePar,buffersize,"Cuts: \n") ;
186 parList+=onePar ;
187 snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
188 parList+=onePar ;
189 snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
190 parList+=onePar ;
191 snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
192 parList+=onePar ;
db2bf6fd 193 if(fMultiCutAna){
194 snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
195 for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
196 parList+=onePar ;
db2bf6fd 197 snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
198 for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
199 parList+=onePar ;
db2bf6fd 200 }
201
af7b3903 202 return new TObjString(parList) ;
0c1383b5 203}
204
2e557d1c 205//________________________________________________________________________________________________________________________________________________
206TList * AliAnaPi0::GetCreateOutputObjects()
207{
477d6cee 208 // Create histograms to be saved in output file and
209 // store them in fOutputContainer
210
211 //create event containers
72542aba 212 fEventsList = new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ;
1c5acb87 213
72542aba 214 for(Int_t ic=0; ic<GetNCentrBin(); ic++){
5025c139 215 for(Int_t iz=0; iz<GetNZvertBin(); iz++){
216 for(Int_t irp=0; irp<GetNRPBin(); irp++){
217 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
af7b3903 218 fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->SetOwner(kFALSE);
477d6cee 219 }
220 }
221 }
7e7694bb 222
477d6cee 223 TList * outputContainer = new TList() ;
224 outputContainer->SetName(GetName());
6921fa00 225
0333ede6 226 fhReMod = new TH2F*[fNModules] ;
227 fhMiMod = new TH2F*[fNModules] ;
228
8d230fa8 229 if(fCalorimeter == "PHOS"){
0333ede6 230 fhReDiffPHOSMod = new TH2F*[fNModules] ;
231 fhMiDiffPHOSMod = new TH2F*[fNModules] ;
8d230fa8 232 }
233 else{
0333ede6 234 fhReSameSectorEMCALMod = new TH2F*[fNModules/2] ;
235 fhReSameSideEMCALMod = new TH2F*[fNModules-2] ;
236 fhMiSameSectorEMCALMod = new TH2F*[fNModules/2] ;
237 fhMiSameSideEMCALMod = new TH2F*[fNModules-2] ;
8d230fa8 238 }
6175da48 239
af7b3903 240
0333ede6 241 fhRe1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
242 fhMi1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
6175da48 243 if(fFillBadDistHisto){
0333ede6 244 fhRe2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
245 fhRe3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
246 fhMi2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
247 fhMi3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
6175da48 248 }
249 if(fMakeInvPtPlots) {
0333ede6 250 fhReInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
251 fhMiInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
6175da48 252 if(fFillBadDistHisto){
0333ede6 253 fhReInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
254 fhReInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
255 fhMiInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
256 fhMiInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
6175da48 257 }
398c93cc 258 }
6175da48 259
5ae09196 260 const Int_t buffersize = 255;
261 char key[buffersize] ;
262 char title[buffersize] ;
477d6cee 263
5a2dbc3c 264 Int_t nptbins = GetHistoPtBins();
265 Int_t nphibins = GetHistoPhiBins();
266 Int_t netabins = GetHistoEtaBins();
267 Float_t ptmax = GetHistoPtMax();
268 Float_t phimax = GetHistoPhiMax();
269 Float_t etamax = GetHistoEtaMax();
270 Float_t ptmin = GetHistoPtMin();
271 Float_t phimin = GetHistoPhiMin();
272 Float_t etamin = GetHistoEtaMin();
273
274 Int_t nmassbins = GetHistoMassBins();
275 Int_t nasymbins = GetHistoAsymmetryBins();
276 Float_t massmax = GetHistoMassMax();
277 Float_t asymmax = GetHistoAsymmetryMax();
278 Float_t massmin = GetHistoMassMin();
279 Float_t asymmin = GetHistoAsymmetryMin();
af7b3903 280 Int_t ntrmbins = GetHistoTrackMultiplicityBins();
281 Int_t ntrmmax = GetHistoTrackMultiplicityMax();
282 Int_t ntrmmin = GetHistoTrackMultiplicityMin();
6175da48 283
72542aba 284 if(GetNCentrBin() > 1 && (fUseAverCellEBins||fUseAverClusterEBins||fUseAverClusterEDenBins)){
20218aea 285
286 fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
287 fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
288 outputContainer->Add(fhAverTotECluster) ;
289
290 fhAverTotECell = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
291 fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
292 outputContainer->Add(fhAverTotECell) ;
293
294 fhAverTotECellvsCluster = new TH2F("hAverTotECellvsCluster","hAverTotECellvsCluster",200,0,50,200,0,50) ;
295 fhAverTotECellvsCluster->SetYTitle("E_{cell, aver. SM} (GeV)");
296 fhAverTotECellvsCluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
297 outputContainer->Add(fhAverTotECellvsCluster) ;
298
299 fhEDensityCluster = new TH1F("hEDensityCluster","hEDensityCluster",200,0,50) ;
300 fhEDensityCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
301 outputContainer->Add(fhEDensityCluster) ;
302
303 fhEDensityCell = new TH1F("hEDensityCell","hEDensityCell",200,0,50) ;
304 fhEDensityCell->SetXTitle("#Sigma E_{cell} / N_{cell} (GeV)");
305 outputContainer->Add(fhEDensityCell) ;
306
307 fhEDensityCellvsCluster = new TH2F("hEDensityCellvsCluster","hEDensityCellvsCluster",200,0,50,200,0,50) ;
308 fhEDensityCellvsCluster->SetYTitle("#Sigma E_{cell} / N_{cell} (GeV)");
309 fhEDensityCellvsCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
310 outputContainer->Add(fhEDensityCellvsCluster) ;
311
20218aea 312 }//counting and average histograms
156549ae 313
20218aea 314 if(fCheckConversion){
0333ede6 315 fhReConv = new TH2F("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 316 fhReConv->SetXTitle("p_{T} (GeV/c)");
317 fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
318 outputContainer->Add(fhReConv) ;
319
0333ede6 320 fhReConv2 = new TH2F("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 321 fhReConv2->SetXTitle("p_{T} (GeV/c)");
322 fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
323 outputContainer->Add(fhReConv2) ;
324
325 if(fDoOwnMix){
0333ede6 326 fhMiConv = new TH2F("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 327 fhMiConv->SetXTitle("p_{T} (GeV/c)");
328 fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
329 outputContainer->Add(fhMiConv) ;
330
0333ede6 331 fhMiConv2 = new TH2F("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 332 fhMiConv2->SetXTitle("p_{T} (GeV/c)");
333 fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
334 outputContainer->Add(fhMiConv2) ;
335 }
156549ae 336 }
6175da48 337
72542aba 338 for(Int_t ic=0; ic<GetNCentrBin(); ic++){
0333ede6 339 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
340 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
341 Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
342 //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
343 //Distance to bad module 1
344 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
345 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
346 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
347 fhRe1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
348 fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
349 fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
350 //printf("name: %s\n ",fhRe1[index]->GetName());
351 outputContainer->Add(fhRe1[index]) ;
352
353 if(fFillBadDistHisto){
354 //Distance to bad module 2
355 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
356 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
357 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
358 fhRe2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
359 fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
360 fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
361 outputContainer->Add(fhRe2[index]) ;
362
363 //Distance to bad module 3
364 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
365 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
366 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
367 fhRe3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
368 fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
369 fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
370 outputContainer->Add(fhRe3[index]) ;
371 }
372
373 //Inverse pT
374 if(fMakeInvPtPlots){
398c93cc 375 //Distance to bad module 1
0333ede6 376 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
398c93cc 377 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
378 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
0333ede6 379 fhReInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
380 fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
381 fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
382 outputContainer->Add(fhReInvPt1[index]) ;
af7b3903 383
6175da48 384 if(fFillBadDistHisto){
398c93cc 385 //Distance to bad module 2
0333ede6 386 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
6175da48 387 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
398c93cc 388 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
0333ede6 389 fhReInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
390 fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
391 fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
392 outputContainer->Add(fhReInvPt2[index]) ;
398c93cc 393
394 //Distance to bad module 3
0333ede6 395 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
6175da48 396 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
398c93cc 397 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
0333ede6 398 fhReInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
399 fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
400 fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
401 outputContainer->Add(fhReInvPt3[index]) ;
398c93cc 402 }
0333ede6 403 }
404 if(fDoOwnMix){
405 //Distance to bad module 1
406 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
407 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
408 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
409 fhMi1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
410 fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
411 fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
412 outputContainer->Add(fhMi1[index]) ;
413 if(fFillBadDistHisto){
414 //Distance to bad module 2
415 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
416 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
6175da48 417 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
0333ede6 418 fhMi2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
419 fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
420 fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
421 outputContainer->Add(fhMi2[index]) ;
6175da48 422
0333ede6 423 //Distance to bad module 3
424 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
425 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
426 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
427 fhMi3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
428 fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
429 fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
430 outputContainer->Add(fhMi3[index]) ;
6175da48 431 }
0333ede6 432 //Inverse pT
433 if(fMakeInvPtPlots){
6175da48 434 //Distance to bad module 1
0333ede6 435 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
6175da48 436 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
437 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
0333ede6 438 fhMiInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
439 fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
440 fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
441 outputContainer->Add(fhMiInvPt1[index]) ;
6175da48 442 if(fFillBadDistHisto){
443 //Distance to bad module 2
0333ede6 444 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
6175da48 445 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
446 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
0333ede6 447 fhMiInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
448 fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
449 fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
450 outputContainer->Add(fhMiInvPt2[index]) ;
6175da48 451
452 //Distance to bad module 3
0333ede6 453 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
454 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
6175da48 455 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
0333ede6 456 fhMiInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
457 fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
458 fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
459 outputContainer->Add(fhMiInvPt3[index]) ;
6175da48 460 }
0333ede6 461 }
462 }
7e7694bb 463 }
1c5acb87 464 }
0333ede6 465 }
477d6cee 466
0333ede6 467 fhRePtAsym = new TH2F("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 468 fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
469 fhRePtAsym->SetYTitle("Asymmetry");
470 outputContainer->Add(fhRePtAsym);
471
0333ede6 472 fhRePtAsymPi0 = new TH2F("hRePtAsymPi0","Asymmetry vs pt, for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 473 fhRePtAsymPi0->SetXTitle("p_{T} (GeV/c)");
474 fhRePtAsymPi0->SetYTitle("Asymmetry");
475 outputContainer->Add(fhRePtAsymPi0);
0333ede6 476
477 fhRePtAsymEta = new TH2F("hRePtAsymEta","Asymmetry vs pt, for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
af7b3903 478 fhRePtAsymEta->SetXTitle("p_{T} (GeV/c)");
479 fhRePtAsymEta->SetYTitle("Asymmetry");
480 outputContainer->Add(fhRePtAsymEta);
481
5ae09196 482 if(fMultiCutAna){
483
0333ede6 484 fhRePIDBits = new TH2F*[fNPIDBits];
5ae09196 485 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
486 snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
487 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
0333ede6 488 fhRePIDBits[ipid] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
af7b3903 489 fhRePIDBits[ipid]->SetXTitle("p_{T} (GeV/c)");
490 fhRePIDBits[ipid]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
5ae09196 491 outputContainer->Add(fhRePIDBits[ipid]) ;
492 }// pid bit loop
493
0333ede6 494 fhRePtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
495 fhMiPtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
496
9302260a 497 if(fFillSMCombinations){
0333ede6 498 for(Int_t iSM = 0; iSM < fNModules; iSM++) fhRePtNCellAsymCutsSM[iSM] = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
499
9302260a 500 }
501
5ae09196 502 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
503 for(Int_t icell=0; icell<fNCellNCuts; icell++){
504 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
505 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
af7b3903 506 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
5ae09196 507 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
508 //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
0333ede6 509 fhRePtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
af7b3903 510 fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
511 fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
5ae09196 512 outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
6175da48 513
514 snprintf(key, buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
515 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
0333ede6 516 fhMiPtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
6175da48 517 fhMiPtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
518 fhMiPtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
9302260a 519 outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
6175da48 520
0333ede6 521 if(fFillSMCombinations){
522 for(Int_t iSM = 0; iSM < fNModules; iSM++){
523 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM%d",ipt,icell,iasym,iSM) ;
524 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM %d ",
525 fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym],iSM) ;
526 fhRePtNCellAsymCutsSM[iSM][index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
527 fhRePtNCellAsymCutsSM[iSM][index]->SetXTitle("p_{T} (GeV/c)");
528 fhRePtNCellAsymCutsSM[iSM][index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
529 outputContainer->Add(fhRePtNCellAsymCutsSM[iSM][index]) ;
530 }
9302260a 531
9302260a 532 }
5ae09196 533 }
534 }
535 }
821c8090 536
9302260a 537 if(ntrmbins!=0){
0333ede6 538 fhRePtMult = new TH3F*[fNAsymCuts] ;
9302260a 539 for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++){
0333ede6 540 fhRePtMult[iasym] = new TH3F(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
9302260a 541 nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
542 fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
543 fhRePtMult[iasym]->SetYTitle("Track multiplicity");
544 fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
545 outputContainer->Add(fhRePtMult[iasym]) ;
546 }
af7b3903 547 }
5ae09196 548 }// multi cuts analysis
549
0333ede6 550 fhReSS[0] = new TH2F("hRe_SS_Tight"," 0.01 < #lambda_{0}^{2} < 0.4",
551 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
552 fhReSS[0]->SetXTitle("p_{T} (GeV/c)");
553 fhReSS[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
554 outputContainer->Add(fhReSS[0]) ;
555
556
557 fhReSS[1] = new TH2F("hRe_SS_Loose"," #lambda_{0}^{2} > 0.4",
558 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
559 fhReSS[1]->SetXTitle("p_{T} (GeV/c)");
560 fhReSS[1]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
561 outputContainer->Add(fhReSS[1]) ;
562
563
564 fhReSS[2] = new TH2F("hRe_SS_Both"," cluster_{1} #lambda_{0}^{2} > 0.4; cluster_{2} 0.01 < #lambda_{0}^{2} < 0.4",
565 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
566 fhReSS[2]->SetXTitle("p_{T} (GeV/c)");
567 fhReSS[2]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
568 outputContainer->Add(fhReSS[2]) ;
569
570 fhEvents=new TH3F("hEvents","Number of events",GetNCentrBin(),0.,1.*GetNCentrBin(),
5025c139 571 GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
745f04da 572
573 fhEvents->SetXTitle("Centrality bin");
574 fhEvents->SetYTitle("Z vertex bin bin");
575 fhEvents->SetZTitle("RP bin");
477d6cee 576 outputContainer->Add(fhEvents) ;
50f39b97 577
72542aba 578 if(GetNCentrBin()>1){
0333ede6 579 fhCentrality=new TH1F("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
20218aea 580 fhCentrality->SetXTitle("Centrality bin");
581 outputContainer->Add(fhCentrality) ;
582
0333ede6 583 fhCentralityNoPair=new TH1F("hCentralityBinNoPair","Number of events in centrality bin, with no cluster pairs",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
20218aea 584 fhCentralityNoPair->SetXTitle("Centrality bin");
585 outputContainer->Add(fhCentralityNoPair) ;
586 }
587
72542aba 588 if(GetNRPBin() > 1 ){
0333ede6 589
590 fhEventPlaneAngle=new TH1F("hEventPlaneAngleBin","Number of events in centrality bin",100,0.,TMath::TwoPi()) ;
72542aba 591 fhEventPlaneAngle->SetXTitle("EP angle (rad)");
592 outputContainer->Add(fhEventPlaneAngle) ;
593
594 if(GetNCentrBin()>1){
0333ede6 595 fhEventPlaneResolution=new TH2F("hEventPlaneResolution","Event plane resolution",GetNCentrBin(),0,GetNCentrBin(),100,0.,TMath::TwoPi()) ;
72542aba 596 fhEventPlaneResolution->SetYTitle("Resolution");
597 fhEventPlaneResolution->SetXTitle("Centrality Bin");
598 outputContainer->Add(fhEventPlaneResolution) ;
599 }
600 }
601
0333ede6 602 fhRealOpeningAngle = new TH2F
6175da48 603 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
50f39b97 604 fhRealOpeningAngle->SetYTitle("#theta(rad)");
605 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
606 outputContainer->Add(fhRealOpeningAngle) ;
7e7694bb 607
0333ede6 608 fhRealCosOpeningAngle = new TH2F
6175da48 609 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
50f39b97 610 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
611 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
612 outputContainer->Add(fhRealCosOpeningAngle) ;
613
6175da48 614 if(fDoOwnMix){
615
0333ede6 616 fhMixedOpeningAngle = new TH2F
6175da48 617 ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
618 fhMixedOpeningAngle->SetYTitle("#theta(rad)");
619 fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
620 outputContainer->Add(fhMixedOpeningAngle) ;
621
0333ede6 622 fhMixedCosOpeningAngle = new TH2F
6175da48 623 ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
624 fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
625 fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
626 outputContainer->Add(fhMixedCosOpeningAngle) ;
627
628 }
629
477d6cee 630 //Histograms filled only if MC data is requested
0ae57829 631 if(IsDataMC()){
99b8e903 632
633 fhReMCFromConversion = new TH2F("hReMCFromConversion","Invariant mass of 2 clusters originated in conversions",
634 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
635 fhReMCFromConversion->SetXTitle("p_{T} (GeV/c)");
636 fhReMCFromConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
637 outputContainer->Add(fhReMCFromConversion) ;
638
639 fhReMCFromNotConversion = new TH2F("hReMCNotFromConversion","Invariant mass of 2 clusters not originated in conversions",
640 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
641 fhReMCFromNotConversion->SetXTitle("p_{T} (GeV/c)");
642 fhReMCFromNotConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
643 outputContainer->Add(fhReMCFromNotConversion) ;
644
645 fhReMCFromMixConversion = new TH2F("hReMCFromMixConversion","Invariant mass of 2 clusters one from conversion and the other not",
646 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
647 fhReMCFromMixConversion->SetXTitle("p_{T} (GeV/c)");
648 fhReMCFromMixConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
649 outputContainer->Add(fhReMCFromMixConversion) ;
650
6175da48 651 //Pi0
0333ede6 652 fhPrimPi0Pt = new TH1F("hPrimPi0Pt","Primary pi0 pt, Y<1",nptbins,ptmin,ptmax) ;
653 fhPrimPi0AccPt = new TH1F("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
156549ae 654 fhPrimPi0Pt ->SetXTitle("p_{T} (GeV/c)");
655 fhPrimPi0AccPt->SetXTitle("p_{T} (GeV/c)");
6175da48 656 outputContainer->Add(fhPrimPi0Pt) ;
657 outputContainer->Add(fhPrimPi0AccPt) ;
658
b529da89 659 Int_t netabinsopen = TMath::Nint(netabins*4/(etamax-etamin));
0333ede6 660 fhPrimPi0Y = new TH2F("hPrimPi0Rapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
156549ae 661 fhPrimPi0Y ->SetYTitle("Rapidity");
662 fhPrimPi0Y ->SetXTitle("p_{T} (GeV/c)");
6175da48 663 outputContainer->Add(fhPrimPi0Y) ;
664
0333ede6 665 fhPrimPi0AccY = new TH2F("hPrimPi0AccRapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
156549ae 666 fhPrimPi0AccY->SetYTitle("Rapidity");
667 fhPrimPi0AccY->SetXTitle("p_{T} (GeV/c)");
6175da48 668 outputContainer->Add(fhPrimPi0AccY) ;
477d6cee 669
fb516de2 670 Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin));
0333ede6 671 fhPrimPi0Phi = new TH2F("hPrimPi0Phi","Azimuthal of primary pi0, Y<1",nptbins,ptmin,ptmax,nphibinsopen,0,360) ;
156549ae 672 fhPrimPi0Phi->SetYTitle("#phi (deg)");
673 fhPrimPi0Phi->SetXTitle("p_{T} (GeV/c)");
6175da48 674 outputContainer->Add(fhPrimPi0Phi) ;
477d6cee 675
0333ede6 676 fhPrimPi0AccPhi = new TH2F("hPrimPi0AccPhi","Azimuthal of primary pi0 with accepted daughters",nptbins,ptmin,ptmax,
677 nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
156549ae 678 fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
679 fhPrimPi0AccPhi->SetXTitle("p_{T} (GeV/c)");
6175da48 680 outputContainer->Add(fhPrimPi0AccPhi) ;
477d6cee 681
6175da48 682 //Eta
0333ede6 683 fhPrimEtaPt = new TH1F("hPrimEtaPt","Primary eta pt",nptbins,ptmin,ptmax) ;
684 fhPrimEtaAccPt = new TH1F("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
156549ae 685 fhPrimEtaPt ->SetXTitle("p_{T} (GeV/c)");
686 fhPrimEtaAccPt->SetXTitle("p_{T} (GeV/c)");
6175da48 687 outputContainer->Add(fhPrimEtaPt) ;
688 outputContainer->Add(fhPrimEtaAccPt) ;
477d6cee 689
0333ede6 690 fhPrimEtaY = new TH2F("hPrimEtaRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
156549ae 691 fhPrimEtaY->SetYTitle("Rapidity");
692 fhPrimEtaY->SetXTitle("p_{T} (GeV/c)");
6175da48 693 outputContainer->Add(fhPrimEtaY) ;
50f39b97 694
0333ede6 695 fhPrimEtaAccY = new TH2F("hPrimEtaAccRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
156549ae 696 fhPrimEtaAccY->SetYTitle("Rapidity");
697 fhPrimEtaAccY->SetXTitle("p_{T} (GeV/c)");
6175da48 698 outputContainer->Add(fhPrimEtaAccY) ;
50f39b97 699
0333ede6 700 fhPrimEtaPhi = new TH2F("hPrimEtaPhi","Azimuthal of primary eta",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
156549ae 701 fhPrimEtaPhi->SetYTitle("#phi (deg)");
702 fhPrimEtaPhi->SetXTitle("p_{T} (GeV/c)");
6175da48 703 outputContainer->Add(fhPrimEtaPhi) ;
704
0333ede6 705 fhPrimEtaAccPhi = new TH2F("hPrimEtaAccPhi","Azimuthal of primary eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
156549ae 706 fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
707 fhPrimEtaAccPhi->SetXTitle("p_{T} (GeV/c)");
6175da48 708 outputContainer->Add(fhPrimEtaAccPhi) ;
0333ede6 709
50f39b97 710
08a56f5f 711 //Prim origin
712 //Pi0
0333ede6 713 fhPrimPi0PtOrigin = new TH2F("hPrimPi0PtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
08a56f5f 714 fhPrimPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
715 fhPrimPi0PtOrigin->SetYTitle("Origin");
716 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
717 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
718 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
719 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
720 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
721 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
722 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
723 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
724 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
725 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
726 outputContainer->Add(fhPrimPi0PtOrigin) ;
727
0333ede6 728 fhMCPi0PtOrigin = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
08a56f5f 729 fhMCPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
730 fhMCPi0PtOrigin->SetYTitle("Origin");
731 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
732 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
733 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
734 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
735 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
736 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
737 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
738 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
739 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
740 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
741 outputContainer->Add(fhMCPi0PtOrigin) ;
742
743 //Eta
0333ede6 744 fhPrimEtaPtOrigin = new TH2F("hPrimEtaPtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
08a56f5f 745 fhPrimEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
746 fhPrimEtaPtOrigin->SetYTitle("Origin");
747 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
748 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
749 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
750 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
751 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
752 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
0333ede6 753
08a56f5f 754 outputContainer->Add(fhPrimEtaPtOrigin) ;
755
0333ede6 756 fhMCEtaPtOrigin = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
08a56f5f 757 fhMCEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
758 fhMCEtaPtOrigin->SetYTitle("Origin");
759 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
760 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
761 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
762 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
763 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
764 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
0333ede6 765
08a56f5f 766 outputContainer->Add(fhMCEtaPtOrigin) ;
08a56f5f 767
0333ede6 768
769 fhPrimPi0OpeningAngle = new TH2F
6175da48 770 ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
771 fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
772 fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
773 outputContainer->Add(fhPrimPi0OpeningAngle) ;
774
0333ede6 775 fhPrimPi0CosOpeningAngle = new TH2F
6175da48 776 ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
777 fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
778 fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
779 outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
780
781 for(Int_t i = 0; i<13; i++){
0333ede6 782 fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
6175da48 783 fhMCOrgMass[i]->SetXTitle("p_{T} (GeV/c)");
784 fhMCOrgMass[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
785 outputContainer->Add(fhMCOrgMass[i]) ;
786
0333ede6 787 fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
6175da48 788 fhMCOrgAsym[i]->SetXTitle("p_{T} (GeV/c)");
789 fhMCOrgAsym[i]->SetYTitle("A");
790 outputContainer->Add(fhMCOrgAsym[i]) ;
791
0333ede6 792 fhMCOrgDeltaEta[i] = new TH2F(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs pt, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
6175da48 793 fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
156549ae 794 fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
6175da48 795 outputContainer->Add(fhMCOrgDeltaEta[i]) ;
796
0333ede6 797 fhMCOrgDeltaPhi[i]= new TH2F(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs p_{T}, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
6175da48 798 fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
156549ae 799 fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
6175da48 800 outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
801
802 }
50f39b97 803
6175da48 804 if(fMultiCutAnaSim){
0333ede6 805 fhMCPi0MassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
806 fhMCPi0MassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
807 fhMCPi0PtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
808 fhMCEtaMassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
809 fhMCEtaMassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
810 fhMCEtaPtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
6175da48 811 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
812 for(Int_t icell=0; icell<fNCellNCuts; icell++){
813 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
814 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
815
0333ede6 816 fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
6175da48 817 Form("Reconstructed Mass vs reconstructed p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
818 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
819 fhMCPi0MassPtRec[index]->SetXTitle("p_{T, reconstructed} (GeV/c)");
820 fhMCPi0MassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
821 outputContainer->Add(fhMCPi0MassPtRec[index]) ;
822
0333ede6 823 fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
6175da48 824 Form("Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
825 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
826 fhMCPi0MassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
827 fhMCPi0MassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
828 outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
829
0333ede6 830 fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
6175da48 831 Form("Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
832 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
833 fhMCPi0PtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
834 fhMCPi0PtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
835 outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
836
0333ede6 837 fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
6175da48 838 Form("Reconstructed Mass vs reconstructed p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
839 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
840 fhMCEtaMassPtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
841 fhMCEtaMassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
842 outputContainer->Add(fhMCEtaMassPtRec[index]) ;
843
0333ede6 844 fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
6175da48 845 Form("Reconstructed Mass vs generated p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
846 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
847 fhMCEtaMassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
848 fhMCEtaMassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
849 outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
850
0333ede6 851 fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
6175da48 852 Form("Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
853 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
854 fhMCEtaPtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
855 fhMCEtaPtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
856 outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
857 }
858 }
859 }
860 }//multi cut ana
861 else {
0333ede6 862 fhMCPi0MassPtTrue = new TH2F*[1];
863 fhMCPi0PtTruePtRec = new TH2F*[1];
864 fhMCEtaMassPtTrue = new TH2F*[1];
865 fhMCEtaPtTruePtRec = new TH2F*[1];
6175da48 866
0333ede6 867 fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
6175da48 868 fhMCPi0MassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
869 fhMCPi0MassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
870 outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
871
0333ede6 872 fhMCPi0PtTruePtRec[0]= new TH2F("hMCPi0PtTruePtRec","Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
6175da48 873 fhMCPi0PtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
874 fhMCPi0PtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
875 outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
876
0333ede6 877 fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
6175da48 878 fhMCEtaMassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
879 fhMCEtaMassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
880 outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
881
0333ede6 882 fhMCEtaPtTruePtRec[0]= new TH2F("hMCEtaPtTruePtRec","Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
6175da48 883 fhMCEtaPtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
884 fhMCEtaPtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
885 outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
886 }
477d6cee 887 }
50f39b97 888
20218aea 889 if(fFillSMCombinations){
890 TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
891 for(Int_t imod=0; imod<fNModules; imod++){
892 //Module dependent invariant mass
893 snprintf(key, buffersize,"hReMod_%d",imod) ;
894 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
0333ede6 895 fhReMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 896 fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
897 fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
898 outputContainer->Add(fhReMod[imod]) ;
8d230fa8 899 if(fCalorimeter=="PHOS"){
20218aea 900 snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
901 snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
0333ede6 902 fhReDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 903 fhReDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
904 fhReDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
905 outputContainer->Add(fhReDiffPHOSMod[imod]) ;
906 }
8d230fa8 907 else{//EMCAL
908 if(imod<fNModules/2){
20218aea 909 snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
910 snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
0333ede6 911 fhReSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 912 fhReSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
913 fhReSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
914 outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
8d230fa8 915 }
916 if(imod<fNModules-2){
20218aea 917 snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
918 snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
0333ede6 919 fhReSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 920 fhReSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
921 fhReSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
922 outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
8d230fa8 923 }
20218aea 924 }//EMCAL
925
926 if(fDoOwnMix){
927 snprintf(key, buffersize,"hMiMod_%d",imod) ;
928 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
0333ede6 929 fhMiMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 930 fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
931 fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
932 outputContainer->Add(fhMiMod[imod]) ;
933
934 if(fCalorimeter=="PHOS"){
935 snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
936 snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
0333ede6 937 fhMiDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 938 fhMiDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
939 fhMiDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
940 outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
941 }//PHOS
942 else{//EMCAL
943 if(imod<fNModules/2){
944 snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
945 snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
0333ede6 946 fhMiSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 947 fhMiSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
948 fhMiSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
949 outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
950 }
951 if(imod<fNModules-2){
952 snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
953 snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
0333ede6 954 fhMiSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
20218aea 955 fhMiSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
956 fhMiSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
957 outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
958 }
959 }//EMCAL
960 }// own mix
961 }//loop combinations
962 } // SM combinations
0333ede6 963
964 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
965 //
966 // printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
967 //
968 // }
eee5fcf1 969
477d6cee 970 return outputContainer;
1c5acb87 971}
972
973//_________________________________________________________________________________________________________________________________________________
974void AliAnaPi0::Print(const Option_t * /*opt*/) const
975{
477d6cee 976 //Print some relevant parameters set for the analysis
a3aebfff 977 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
477d6cee 978 AliAnaPartCorrBaseClass::Print(" ");
0333ede6 979
72542aba 980 printf("Number of bins in Centrality: %d \n",GetNCentrBin()) ;
5025c139 981 printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
982 printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
72542aba 983 printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
af7b3903 984 printf("Pair in same Module: %d \n",fSameSM) ;
477d6cee 985 printf("Cuts: \n") ;
0333ede6 986 // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
50f39b97 987 printf("Number of modules: %d \n",fNModules) ;
6175da48 988 printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
af7b3903 989 printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
990 printf("\tasymmetry < ");
991 for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
992 printf("\n");
993
994 printf("PID selection bits: n = %d, \n",fNPIDBits) ;
995 printf("\tPID bit = ");
996 for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
997 printf("\n");
998
db2bf6fd 999 if(fMultiCutAna){
1000 printf("pT cuts: n = %d, \n",fNPtCuts) ;
1001 printf("\tpT > ");
1002 for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
1003 printf("GeV/c\n");
1004
1005 printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
1006 printf("\tnCell > ");
1007 for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
1008 printf("\n");
0333ede6 1009
db2bf6fd 1010 }
477d6cee 1011 printf("------------------------------------------------------\n") ;
1c5acb87 1012}
1013
5ae09196 1014//_____________________________________________________________
1015void AliAnaPi0::FillAcceptanceHistograms(){
1016 //Fill acceptance histograms if MC data is available
c8fe2783 1017
6175da48 1018 if(GetReader()->ReadStack()){
5ae09196 1019 AliStack * stack = GetMCStack();
6175da48 1020 if(stack){
fb516de2 1021 for(Int_t i=0 ; i<stack->GetNtrack(); i++){
5ae09196 1022 TParticle * prim = stack->Particle(i) ;
6175da48 1023 Int_t pdg = prim->GetPdgCode();
fb516de2 1024 //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
0333ede6 1025 // prim->GetName(), prim->GetPdgCode());
1026
6175da48 1027 if( pdg == 111 || pdg == 221){
5ae09196 1028 Double_t pi0Pt = prim->Pt() ;
5ae09196 1029 if(prim->Energy() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
1030 Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
1031 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
6175da48 1032 if(pdg == 111){
6eb99ccd 1033 if(TMath::Abs(pi0Y) < 1.0){
91e1ea12 1034 fhPrimPi0Pt ->Fill(pi0Pt) ;
1035 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
6175da48 1036 }
91e1ea12 1037 fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
6175da48 1038 }
1039 else if(pdg == 221){
6eb99ccd 1040 if(TMath::Abs(pi0Y) < 1.0){
fb516de2 1041 fhPrimEtaPt ->Fill(pi0Pt) ;
1042 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
6175da48 1043 }
156549ae 1044 fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
5ae09196 1045 }
08a56f5f 1046
1047 //Origin of meson
1048 Int_t momindex = prim->GetFirstMother();
04131edb 1049 if(momindex >= 0) {
1050 TParticle* mother = stack->Particle(momindex);
1051 Int_t mompdg = TMath::Abs(mother->GetPdgCode());
1052 Int_t momstatus = mother->GetStatusCode();
1053 if(pdg == 111){
1054 if (momstatus == 21)fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1055 else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1056 else if(mompdg > 2100 && mompdg < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1057 else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1058 else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1059 else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1060 else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1061 else if(mompdg >= 310 && mompdg <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1062 else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1063 else if(momstatus == 11 || momstatus == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
91e1ea12 1064 else fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
04131edb 1065 }//pi0
1066 else {
1067 if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1068 else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1069 else if(mompdg > 2100 && mompdg < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1070 else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1071 else if(momstatus == 11 || momstatus == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1072 else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1073 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1074 }
1075 } // pi0 has mother
08a56f5f 1076
5ae09196 1077 //Check if both photons hit Calorimeter
fb516de2 1078 if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
5ae09196 1079 Int_t iphot1=prim->GetFirstDaughter() ;
1080 Int_t iphot2=prim->GetLastDaughter() ;
1081 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
1082 TParticle * phot1 = stack->Particle(iphot1) ;
1083 TParticle * phot2 = stack->Particle(iphot2) ;
1084 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
1085 //printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
1086 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
1087
1088 TLorentzVector lv1, lv2;
1089 phot1->Momentum(lv1);
1090 phot2->Momentum(lv2);
5ae09196 1091 Bool_t inacceptance = kFALSE;
1092 if(fCalorimeter == "PHOS"){
1093 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
1094 Int_t mod ;
1095 Double_t x,z ;
1096 if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
1097 inacceptance = kTRUE;
1098 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1099 }
1100 else{
1101
1102 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1103 inacceptance = kTRUE ;
1104 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1105 }
1106
1107 }
1108 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
1109 if(GetEMCALGeometry()){
156549ae 1110
1111 Int_t absID1=0;
1112 Int_t absID2=0;
1113
1114 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1115 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
0333ede6 1116
156549ae 1117 if( absID1 >= 0 && absID2 >= 0)
5ae09196 1118 inacceptance = kTRUE;
156549ae 1119
0333ede6 1120 // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1121 // inacceptance = kTRUE;
5ae09196 1122 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1123 }
1124 else{
1125 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1126 inacceptance = kTRUE ;
1127 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1128 }
1129 }
1130
1131 if(inacceptance){
6175da48 1132 if(pdg==111){
91e1ea12 1133 fhPrimPi0AccPt ->Fill(pi0Pt) ;
1134 fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1135 fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
6175da48 1136 Double_t angle = lv1.Angle(lv2.Vect());
91e1ea12 1137 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1138 fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
6175da48 1139 }
1140 else if(pdg==221){
156549ae 1141 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1142 fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1143 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
6175da48 1144 }
5ae09196 1145 }//Accepted
1146 }// 2 photons
1147 }//Check daughters exist
156549ae 1148 }// Primary pi0 or eta
5ae09196 1149 }//loop on primaries
1150 }//stack exists and data is MC
1151 }//read stack
1152 else if(GetReader()->ReadAODMCParticles()){
6175da48 1153 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1154 if(mcparticles){
1155 Int_t nprim = mcparticles->GetEntriesFast();
0333ede6 1156
08a56f5f 1157 for(Int_t i=0; i < nprim; i++)
156549ae 1158 {
08a56f5f 1159 AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
04131edb 1160
1161 // Only generator particles, when they come from PYTHIA, PHOJET, HERWIG ...
1162 //if( prim->GetStatus() == 0 && (GetMCAnalysisUtils()->GetMCGenerator()).Length()!=0) break;
0333ede6 1163
6175da48 1164 Int_t pdg = prim->GetPdgCode();
1165 if( pdg == 111 || pdg == 221){
1166 Double_t pi0Pt = prim->Pt() ;
04131edb 1167 //printf("pi0, pt %2.2f, eta %f, phi %f\n",pi0Pt, prim->Eta(), prim->Phi());
fb516de2 1168 if(prim->E() == TMath::Abs(prim->Pz())) continue ; //Protection against floating point exception
1169
6175da48 1170 Double_t pi0Y = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
1171 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
1172 if(pdg == 111){
fb516de2 1173 if(TMath::Abs(pi0Y) < 1){
91e1ea12 1174 fhPrimPi0Pt->Fill(pi0Pt) ;
1175 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
6175da48 1176 }
91e1ea12 1177 fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
6175da48 1178 }
1179 else if(pdg == 221){
fb516de2 1180 if(TMath::Abs(pi0Y) < 1){
1181 fhPrimEtaPt->Fill(pi0Pt) ;
1182 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
6175da48 1183 }
156549ae 1184 fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
6175da48 1185 }
08a56f5f 1186
1187 //Origin of meson
1188 Int_t momindex = prim->GetMother();
04131edb 1189 if(momindex >= 0) {
1190 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1191 Int_t mompdg = TMath::Abs(mother->GetPdgCode());
1192 Int_t momstatus = mother->GetStatus();
08a56f5f 1193 if(pdg == 111){
1194 if (momstatus == 21) fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1195 else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1196 else if(mompdg > 2100 && mompdg < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1197 else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1198 else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1199 else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1200 else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1201 else if(mompdg >= 310 && mompdg <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1202 else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1203 else if(momstatus == 11 || momstatus == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
91e1ea12 1204 else fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
08a56f5f 1205 }//pi0
1206 else {
1207 if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1208 else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1209 else if(mompdg > 2100 && mompdg < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1210 else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1211 else if(momstatus == 11 || momstatus == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1212 else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1213 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1214 }
04131edb 1215 }//pi0 has mother
08a56f5f 1216
6175da48 1217 //Check if both photons hit Calorimeter
fb516de2 1218 if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
6175da48 1219 Int_t iphot1=prim->GetDaughter(0) ;
1220 Int_t iphot2=prim->GetDaughter(1) ;
1221 if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim){
1222 AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);
1223 AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
1224 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
6175da48 1225 TLorentzVector lv1, lv2;
1226 lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1227 lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
0333ede6 1228
6175da48 1229 Bool_t inacceptance = kFALSE;
1230 if(fCalorimeter == "PHOS"){
1231 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
1232 Int_t mod ;
1233 Double_t x,z ;
1234 Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
1235 Double_t vtx2[]={phot2->Xv(),phot2->Yv(),phot2->Zv()};
1236 if(GetPHOSGeometry()->ImpactOnEmc(vtx, phot1->Theta(),phot1->Phi(),mod,z,x) &&
1237 GetPHOSGeometry()->ImpactOnEmc(vtx2,phot2->Theta(),phot2->Phi(),mod,z,x))
1238 inacceptance = kTRUE;
1239 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1240 }
1241 else{
1242
1243 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1244 inacceptance = kTRUE ;
1245 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1246 }
1247
1248 }
1249 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
1250 if(GetEMCALGeometry()){
0333ede6 1251
6175da48 1252 Int_t absID1=0;
6175da48 1253 Int_t absID2=0;
156549ae 1254
156549ae 1255 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1256 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
0333ede6 1257
6175da48 1258 if( absID1 >= 0 && absID2 >= 0)
1259 inacceptance = kTRUE;
156549ae 1260
156549ae 1261
6175da48 1262 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1263 }
1264 else{
1265 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1266 inacceptance = kTRUE ;
1267 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1268 }
1269 }
1270
1271 if(inacceptance){
1272 if(pdg==111){
0333ede6 1273 // printf("ACCEPTED pi0: pt %2.2f, phi %3.2f, eta %1.2f\n",pi0Pt,phi,pi0Y);
91e1ea12 1274 fhPrimPi0AccPt ->Fill(pi0Pt) ;
1275 fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1276 fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
6175da48 1277 Double_t angle = lv1.Angle(lv2.Vect());
91e1ea12 1278 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1279 fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
6175da48 1280 }
1281 else if(pdg==221){
156549ae 1282 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1283 fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1284 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
6175da48 1285 }
1286 }//Accepted
1287 }// 2 photons
1288 }//Check daughters exist
156549ae 1289 }// Primary pi0 or eta
6175da48 1290 }//loop on primaries
1291 }//stack exists and data is MC
1292
1293
1294 } // read AOD MC
5ae09196 1295}
1296
6175da48 1297//_____________________________________________________________
1298void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1, const Int_t index2,
1299 const Float_t pt1, const Float_t pt2,
1300 const Int_t ncell1, const Int_t ncell2,
1301 const Double_t mass, const Double_t pt, const Double_t asym,
1302 const Double_t deta, const Double_t dphi){
1303 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1304 //Adjusted for Pythia, need to see what to do for other generators.
1305 //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1306 // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1307
1308 Int_t ancPDG = 0;
1309 Int_t ancStatus = 0;
1310 TLorentzVector ancMomentum;
c4a7d28a 1311 TVector3 prodVertex;
6175da48 1312 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
c4a7d28a 1313 GetReader(), ancPDG, ancStatus,ancMomentum, prodVertex);
6175da48 1314
08a56f5f 1315 Int_t momindex = -1;
1316 Int_t mompdg = -1;
1317 Int_t momstatus = -1;
6175da48 1318 if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1319 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
0333ede6 1320
1321 if(ancLabel > -1){
1322 if(ancPDG==22){//gamma
1323 fhMCOrgMass[0]->Fill(pt,mass);
1324 fhMCOrgAsym[0]->Fill(pt,asym);
1325 fhMCOrgDeltaEta[0]->Fill(pt,deta);
1326 fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1327 }
1328 else if(TMath::Abs(ancPDG)==11){//e
1329 fhMCOrgMass[1]->Fill(pt,mass);
1330 fhMCOrgAsym[1]->Fill(pt,asym);
1331 fhMCOrgDeltaEta[1]->Fill(pt,deta);
1332 fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1333 }
1334 else if(ancPDG==111){//Pi0
1335 fhMCOrgMass[2]->Fill(pt,mass);
1336 fhMCOrgAsym[2]->Fill(pt,asym);
1337 fhMCOrgDeltaEta[2]->Fill(pt,deta);
1338 fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
1339 if(fMultiCutAnaSim){
1340 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1341 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1342 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1343 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1344 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1345 asym < fAsymCuts[iasym] &&
1346 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1347 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1348 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1349 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1350 }//pass the different cuts
1351 }// pid bit cut loop
1352 }// icell loop
1353 }// pt cut loop
1354 }//Multi cut ana sim
1355 else {
1356 fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1357 if(mass < 0.17 && mass > 0.1) {
1358 fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
08a56f5f 1359
1360 if(GetReader()->ReadStack()){
1361 TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1362 momindex = ancestor->GetFirstMother();
41121cfe 1363 if(momindex < 0) return;
08a56f5f 1364 TParticle* mother = GetMCStack()->Particle(momindex);
1365 mompdg = TMath::Abs(mother->GetPdgCode());
1366 momstatus = mother->GetStatusCode();
1367 }
1368 else {
1369 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1370 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1371 momindex = ancestor->GetMother();
41121cfe 1372 if(momindex < 0) return;
08a56f5f 1373 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1374 mompdg = TMath::Abs(mother->GetPdgCode());
1375 momstatus = mother->GetStatus();
0333ede6 1376 }
1377
1378 if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1379 else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1380 else if(mompdg > 2100 && mompdg < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1381 else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1382 else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1383 else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1384 else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1385 else if(mompdg >= 310 && mompdg <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1386 else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1387 else if(momstatus == 11 || momstatus == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
1388 else fhMCPi0PtOrigin->Fill(pt,7.5);//other?
08a56f5f 1389
0333ede6 1390 }//pi0 mass region
1391
6175da48 1392 }
0333ede6 1393 }
1394 else if(ancPDG==221){//Eta
1395 fhMCOrgMass[3]->Fill(pt,mass);
1396 fhMCOrgAsym[3]->Fill(pt,asym);
1397 fhMCOrgDeltaEta[3]->Fill(pt,deta);
1398 fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1399 if(fMultiCutAnaSim){
1400 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1401 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1402 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1403 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1404 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1405 asym < fAsymCuts[iasym] &&
1406 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1407 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1408 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1409 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1410 }//pass the different cuts
1411 }// pid bit cut loop
1412 }// icell loop
1413 }// pt cut loop
1414 } //Multi cut ana sim
1415 else {
1416 fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1417 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1418
1419 if(GetReader()->ReadStack()){
1420 TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1421 momindex = ancestor->GetFirstMother();
1422 if(momindex < 0) return;
1423 TParticle* mother = GetMCStack()->Particle(momindex);
1424 mompdg = TMath::Abs(mother->GetPdgCode());
1425 momstatus = mother->GetStatusCode();
1426 }
1427 else {
1428 TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1429 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1430 momindex = ancestor->GetMother();
1431 if(momindex < 0) return;
1432 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1433 mompdg = TMath::Abs(mother->GetPdgCode());
1434 momstatus = mother->GetStatus();
1435 }
1436
1437 if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1438 else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1439 else if(mompdg > 2100 && mompdg < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1440 else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1441 else if(momstatus == 11 || momstatus == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1442 else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1443 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1444 }// eta mass region
1445 }
1446 else if(ancPDG==-2212){//AProton
1447 fhMCOrgMass[4]->Fill(pt,mass);
1448 fhMCOrgAsym[4]->Fill(pt,asym);
1449 fhMCOrgDeltaEta[4]->Fill(pt,deta);
1450 fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1451 }
1452 else if(ancPDG==-2112){//ANeutron
1453 fhMCOrgMass[5]->Fill(pt,mass);
1454 fhMCOrgAsym[5]->Fill(pt,asym);
1455 fhMCOrgDeltaEta[5]->Fill(pt,deta);
1456 fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1457 }
1458 else if(TMath::Abs(ancPDG)==13){//muons
1459 fhMCOrgMass[6]->Fill(pt,mass);
1460 fhMCOrgAsym[6]->Fill(pt,asym);
1461 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1462 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1463 }
1464 else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1465 if(ancStatus==1){//Stable particles, converted? not decayed resonances
6175da48 1466 fhMCOrgMass[6]->Fill(pt,mass);
1467 fhMCOrgAsym[6]->Fill(pt,asym);
1468 fhMCOrgDeltaEta[6]->Fill(pt,deta);
0333ede6 1469 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1470 }
1471 else{//resonances and other decays, more hadron conversions?
1472 fhMCOrgMass[7]->Fill(pt,mass);
1473 fhMCOrgAsym[7]->Fill(pt,asym);
1474 fhMCOrgDeltaEta[7]->Fill(pt,deta);
1475 fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
6175da48 1476 }
6175da48 1477 }
0333ede6 1478 else {//Partons, colliding protons, strings, intermediate corrections
1479 if(ancStatus==11 || ancStatus==12){//String fragmentation
1480 fhMCOrgMass[8]->Fill(pt,mass);
1481 fhMCOrgAsym[8]->Fill(pt,asym);
1482 fhMCOrgDeltaEta[8]->Fill(pt,deta);
1483 fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
1484 }
1485 else if (ancStatus==21){
1486 if(ancLabel < 2) {//Colliding protons
1487 fhMCOrgMass[11]->Fill(pt,mass);
1488 fhMCOrgAsym[11]->Fill(pt,asym);
1489 fhMCOrgDeltaEta[11]->Fill(pt,deta);
1490 fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
1491 }//colliding protons
1492 else if(ancLabel < 6){//partonic initial states interactions
1493 fhMCOrgMass[9]->Fill(pt,mass);
1494 fhMCOrgAsym[9]->Fill(pt,asym);
1495 fhMCOrgDeltaEta[9]->Fill(pt,deta);
1496 fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
1497 }
1498 else if(ancLabel < 8){//Final state partons radiations?
1499 fhMCOrgMass[10]->Fill(pt,mass);
1500 fhMCOrgAsym[10]->Fill(pt,asym);
1501 fhMCOrgDeltaEta[10]->Fill(pt,deta);
1502 fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
1503 }
1504 // else {
1505 // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1506 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1507 // }
1508 }//status 21
1509 //else {
1510 // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1511 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1512 // }
1513 }////Partons, colliding protons, strings, intermediate corrections
1514 }//ancLabel > -1
1515 else { //ancLabel <= -1
1516 //printf("Not related at all label = %d\n",ancLabel);
1517 fhMCOrgMass[12]->Fill(pt,mass);
1518 fhMCOrgAsym[12]->Fill(pt,asym);
1519 fhMCOrgDeltaEta[12]->Fill(pt,deta);
1520 fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
1521 }
6175da48 1522}
1523
20218aea 1524//____________________________________________________________________________________________________________________________________________________
1525void AliAnaPi0::CountAndGetAverages(Int_t &nClus,Int_t &nCell, Float_t &eClusTot,Float_t &eCellTot, Float_t &eDenClus,Float_t &eDenCell) {
0333ede6 1526 // Count the number of clusters and cells, deposited energy, and do some averages in case multiplicity bins dependent on such numbers
1527 // are requested
20218aea 1528 if(fCalorimeter=="EMCAL"){
1529 nClus = GetEMCALClusters() ->GetEntriesFast();
1530 nCell = GetEMCALCells()->GetNumberOfCells();
1531 for(Int_t icl=0; icl < nClus; icl++) {
1532 Float_t e1 = ((AliVCluster*)GetEMCALClusters()->At(icl))->E();
1533 eClusTot += e1;
20218aea 1534 }// first cluster
1535
1536 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetEMCALCells()->GetAmplitude(jce);
0333ede6 1537 }
20218aea 1538 else {
1539 nClus = GetPHOSClusters()->GetEntriesFast();
1540 nCell = GetPHOSCells() ->GetNumberOfCells();
1541 for(Int_t icl=0; icl < nClus; icl++) {
1542 Float_t e1 = ((AliVCluster*)GetPHOSClusters()->At(icl))->E();
1543 eClusTot += e1;
20218aea 1544 }// first cluster
1545 for(Int_t jce=0; jce < nCell; jce++) eCellTot += GetPHOSCells()->GetAmplitude(jce);
0333ede6 1546 }
20218aea 1547 if(GetDebug() > 1)
1548 printf("AliAnaPi0::MakeAnalysisFillHistograms() - # Clusters %d, sum cluster E per SM %f,# Cells %d, sum cell E per SM %f\n", nClus,eClusTot,nCell,eCellTot);
0333ede6 1549
1550 //Fill histograms with "energy density", ncell and nclust will be > 0 since there are at least 2 "photons"
1551 eDenClus = eClusTot/nClus;
1552 eDenCell = eCellTot/nCell;
1553 fhEDensityCluster ->Fill(eDenClus);
1554 fhEDensityCell ->Fill(eDenCell);
1555 fhEDensityCellvsCluster->Fill(eDenClus, eDenCell);
1556 //Fill the average number of cells or clusters per SM
1557 eClusTot /=fNModules;
1558 eCellTot /=fNModules;
1559 fhAverTotECluster ->Fill(eClusTot);
1560 fhAverTotECell ->Fill(eCellTot);
1561 fhAverTotECellvsCluster->Fill(eClusTot, eCellTot);
1562 //printf("Average Cluster: E %f, density %f; Average Cell E %f, density %f\n ",eClusTot,eDenClus,eCellTot,eDenCell);
20218aea 1563}
1564
1c5acb87 1565//____________________________________________________________________________________________________________________________________________________
6639984f 1566void AliAnaPi0::MakeAnalysisFillHistograms()
1c5acb87 1567{
477d6cee 1568 //Process one event and extract photons from AOD branch
1569 // filled with AliAnaPhoton and fill histos with invariant mass
1570
6175da48 1571 //In case of simulated data, fill acceptance histograms
1572 if(IsDataMC())FillAcceptanceHistograms();
08a56f5f 1573
1574 //if (GetReader()->GetEventNumber()%10000 == 0)
1575 // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
1576
6175da48 1577 //Init some variables
6175da48 1578 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
1579 Int_t nClus = 0;
1580 Int_t nCell = 0;
1581 Float_t eClusTot = 0;
1582 Float_t eCellTot = 0;
156549ae 1583 Float_t eDenClus = 0;
1584 Float_t eDenCell = 0;
477d6cee 1585
72542aba 1586 if(GetNCentrBin() > 1 && (fUseAverCellEBins||fUseAverClusterEBins||fUseAverClusterEDenBins))
20218aea 1587 CountAndGetAverages(nClus,nCell,eClusTot,eCellTot,eDenClus,eDenCell);
0333ede6 1588
20218aea 1589
156549ae 1590 if(GetDebug() > 1)
1591 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
0333ede6 1592
156549ae 1593 //If less than photon 2 entries in the list, skip this event
20218aea 1594 if(nPhot < 2 ) {
156549ae 1595
20218aea 1596 if(GetDebug() > 2)
1597 printf("AliAnaPi0::MakeAnalysisFillHistograms() - nPhotons %d, cent bin %d continue to next event\n",nPhot, GetEventCentrality());
1598
72542aba 1599 if(GetNCentrBin() > 1) fhCentralityNoPair->Fill(GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt());
20218aea 1600
1601 return ;
6175da48 1602 }
6175da48 1603
1604 //Init variables
1605 Int_t module1 = -1;
1606 Int_t module2 = -1;
1607 Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
1608 Int_t evtIndex1 = 0 ;
1609 Int_t currentEvtIndex = -1;
1610 Int_t curCentrBin = 0 ;
1611 Int_t curRPBin = 0 ;
1612 Int_t curZvertBin = 0 ;
1613
c4a7d28a 1614 //Get shower shape information of clusters
1615 TObjArray *clusters = 0;
474c8976 1616 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1617 else if(fCalorimeter=="PHOS" ) clusters = GetPHOSClusters() ;
c4a7d28a 1618
6175da48 1619 //---------------------------------
1620 //First loop on photons/clusters
1621 //---------------------------------
477d6cee 1622 for(Int_t i1=0; i1<nPhot-1; i1++){
1623 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
6175da48 1624 //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
0333ede6 1625
7e7694bb 1626 // get the event index in the mixed buffer where the photon comes from
1627 // in case of mixing with analysis frame, not own mixing
c8fe2783 1628 evtIndex1 = GetEventIndex(p1, vert) ;
5025c139 1629 //printf("charge = %d\n", track->Charge());
c8fe2783 1630 if ( evtIndex1 == -1 )
1631 return ;
1632 if ( evtIndex1 == -2 )
1633 continue ;
41121cfe 1634
1635 //printf("z vertex %f < %f\n",vert[2],GetZvertexCut());
2244659d 1636 if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
0333ede6 1637
6175da48 1638
1639 //----------------------------------------------------------------------------
1640 // Get the multiplicity bin. Different cases: centrality (PbPb),
1641 // average cluster multiplicity, average cell multiplicity, track multiplicity
1642 // default is centrality bins
1643 //----------------------------------------------------------------------------
c8fe2783 1644 if (evtIndex1 != currentEvtIndex) {
6175da48 1645 if(fUseTrackMultBins){ // Track multiplicity bins
1646 //printf("track mult %d\n",GetTrackMultiplicity());
1647 curCentrBin = (GetTrackMultiplicity()-1)/5;
72542aba 1648 if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
6175da48 1649 //printf("track mult bin %d\n",curCentrBin);
1650 }
1651 else if(fUsePhotonMultBins){ // Photon multiplicity bins
1652 //printf("photon mult %d cluster mult %d\n",nPhot, nClus);
72542aba 1653 curCentrBin = nPhot-2;
1654 if(curCentrBin > GetNCentrBin() -1) curCentrBin=GetNCentrBin()-1;
156549ae 1655 //printf("photon mult bin %d\n",curRPBin);
6175da48 1656 }
156549ae 1657 else if(fUseAverClusterEBins){ // Cluster average energy bins
6175da48 1658 //Bins for pp, if needed can be done in a more general way
72542aba 1659 curCentrBin = (Int_t) eClusTot/10 * GetNCentrBin();
1660 if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
6175da48 1661 //printf("cluster E average %f, bin %d \n",eClusTot,curCentrBin);
1662 }
156549ae 1663 else if(fUseAverCellEBins){ // Cell average energy bins
6175da48 1664 //Bins for pp, if needed can be done in a more general way
72542aba 1665 curCentrBin = (Int_t) eCellTot/10*GetNCentrBin();
1666 if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
6175da48 1667 //printf("cell E average %f, bin %d \n",eCellTot,curCentrBin);
1668 }
156549ae 1669 else if(fUseAverClusterEDenBins){ // Energy density bins
1670 //Bins for pp, if needed can be done in a more general way
72542aba 1671 curCentrBin = (Int_t) eDenClus/10*GetNCentrBin();
1672 if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
156549ae 1673 //printf("cluster Eden average %f, bin %d \n",eDenClus,curCentrBin);
1674 }
6175da48 1675 else { //Event centrality
0333ede6 1676 // Centrality task returns at maximum 10, 20 or 100, depending on option chosen and
1677 // number of bins, the bin has to be corrected
1678 curCentrBin = GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt();
1679 if(GetDebug() > 0 )printf("AliAnaPi0::MakeAnalysisFillHistograms() - curCentrBin %d, centrality %d, n bins %d, max bin from centrality %d\n",
1680 curCentrBin, GetEventCentrality(), GetNCentrBin(), GetReader()->GetCentralityOpt());
6175da48 1681 }
6eb99ccd 1682
72542aba 1683 if (curCentrBin < 0 || curCentrBin >= GetNCentrBin()){
aba0ecc2 1684 if(GetDebug() > 0)
72542aba 1685 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality bin <%d> not expected, n bins <%d> , return\n",curCentrBin,GetNCentrBin());
6eb99ccd 1686 return;
1687 }
0333ede6 1688
156549ae 1689 //Reaction plane bin
c8fe2783 1690 curRPBin = 0 ;
606cbcb1 1691 if(GetNRPBin()>1 && GetEventPlane()){
72542aba 1692 Float_t epAngle = GetEventPlane()->GetEventplane(GetEventPlaneMethod());
1693 fhEventPlaneAngle->Fill(epAngle);
1694 curRPBin = TMath::Nint(epAngle*(GetNRPBin()-1)/TMath::Pi());
1695 if(curRPBin >= GetNRPBin()) printf("RP Bin %d out of range %d\n",curRPBin,GetNRPBin());
1696 //printf("RP: %d, %f, angle %f, n bin %d\n", curRPBin,epAngle*(GetNRPBin()-1)/TMath::Pi(),epAngle,GetNRPBin());
1697 }
1698
156549ae 1699 //Get vertex z bin
5025c139 1700 curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
6175da48 1701
1702 //Fill event bin info
745f04da 1703 fhEvents ->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
72542aba 1704 if(GetNCentrBin() > 1) {
1705 fhCentrality->Fill(curCentrBin);
606cbcb1 1706 if(GetNRPBin() > 1 && GetEventPlane()) fhEventPlaneResolution->Fill(curCentrBin,TMath::Cos(2.*GetEventPlane()->GetQsubRes()));
72542aba 1707 }
c8fe2783 1708 currentEvtIndex = evtIndex1 ;
ca468d44 1709 if(GetDebug() > 1)
6175da48 1710 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d \n",curCentrBin,curRPBin,curZvertBin);
c8fe2783 1711 }
7e7694bb 1712
f8006433 1713 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
af7b3903 1714
6175da48 1715 //Get the momentum of this cluster
477d6cee 1716 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
6175da48 1717
1718 //Get (Super)Module number of this cluster
59b6bd99 1719 module1 = GetModuleNumber(p1);
6175da48 1720
0333ede6 1721 //------------------------------------------
1722 //Get index in VCaloCluster array
c4a7d28a 1723 AliVCluster *cluster1 = 0;
1724 Bool_t bFound1 = kFALSE;
9ab9e937 1725 Int_t caloLabel1 = p1->GetCaloLabel(0);
1726 Bool_t iclus1 =-1;
1727 if(clusters){
1728 for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
1729 AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
1730 if(cluster){
1731 if (cluster->GetID()==caloLabel1) {
1732 bFound1 = kTRUE ;
1733 cluster1 = cluster;
1734 iclus1 = iclus;
1735 }
1736 }
1737 if(bFound1) break;
1738 }
c4a7d28a 1739 }// calorimeter clusters loop
1740
6175da48 1741 //---------------------------------
1742 //Second loop on photons/clusters
1743 //---------------------------------
477d6cee 1744 for(Int_t i2=i1+1; i2<nPhot; i2++){
1745 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
0333ede6 1746
6175da48 1747 //In case of mixing frame, check we are not in the same event as the first cluster
c8fe2783 1748 Int_t evtIndex2 = GetEventIndex(p2, vert) ;
1749 if ( evtIndex2 == -1 )
1750 return ;
1751 if ( evtIndex2 == -2 )
1752 continue ;
1753 if (GetMixedEvent() && (evtIndex1 == evtIndex2))
7e7694bb 1754 continue ;
0333ede6 1755
9ab9e937 1756 //------------------------------------------
0333ede6 1757 //Get index in VCaloCluster array
c4a7d28a 1758 AliVCluster *cluster2 = 0;
1759 Bool_t bFound2 = kFALSE;
1760 Int_t caloLabel2 = p2->GetCaloLabel(0);
9ab9e937 1761 if(clusters){
1762 for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
1763 AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
1764 if(cluster){
1765 if(cluster->GetID()==caloLabel2) {
1766 bFound2 = kTRUE ;
1767 cluster2 = cluster;
1768 }
1769 }
1770 if(bFound2) break;
1771 }// calorimeter clusters loop
1772 }
c4a7d28a 1773
1774 Float_t tof1 = -1;
0333ede6 1775 Float_t l01 = -1;
c4a7d28a 1776 if(cluster1 && bFound1){
1777 tof1 = cluster1->GetTOF()*1e9;
0333ede6 1778 l01 = cluster1->GetM02();
c4a7d28a 1779 }
0333ede6 1780 // else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
1781 // p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
c4a7d28a 1782
1783 Float_t tof2 = -1;
0333ede6 1784 Float_t l02 = -1;
c4a7d28a 1785 if(cluster2 && bFound2){
1786 tof2 = cluster2->GetTOF()*1e9;
0333ede6 1787 l02 = cluster2->GetM02();
1788
c4a7d28a 1789 }
0333ede6 1790 // else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
1791 // p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
6175da48 1792
9ab9e937 1793 if(clusters){
1794 Double_t t12diff = tof1-tof2;
1795 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
1796 }
1797 //------------------------------------------
0333ede6 1798
f8006433 1799 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
0333ede6 1800
6175da48 1801 //Get the momentum of this cluster
477d6cee 1802 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
59b6bd99 1803 //Get module number
6175da48 1804 module2 = GetModuleNumber(p2);
1805
1806 //---------------------------------
1807 // Get pair kinematics
1808 //---------------------------------
1809 Double_t m = (photon1 + photon2).M() ;
1810 Double_t pt = (photon1 + photon2).Pt();
1811 Double_t deta = photon1.Eta() - photon2.Eta();
1812 Double_t dphi = photon1.Phi() - photon2.Phi();
1813 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1814
477d6cee 1815 if(GetDebug() > 2)
6175da48 1816 printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
1817
1818 //--------------------------------
1819 // Opening angle selection
1820 //--------------------------------
50f39b97 1821 //Check if opening angle is too large or too small compared to what is expected
1822 Double_t angle = photon1.Angle(photon2.Vect());
6175da48 1823 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
1824 if(GetDebug() > 2)
1825 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
c8fe2783 1826 continue;
6175da48 1827 }
0333ede6 1828
6175da48 1829 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
1830 if(GetDebug() > 2)
1831 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
1832 continue;
1833 }
0333ede6 1834
6175da48 1835 //-------------------------------------------------------------------------------------------------
af7b3903 1836 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
6175da48 1837 //-------------------------------------------------------------------------------------------------
20218aea 1838 if(a < fAsymCuts[0] && fFillSMCombinations){
af7b3903 1839 if(module1==module2 && module1 >=0 && module1<fNModules)
91e1ea12 1840 fhReMod[module1]->Fill(pt,m) ;
0333ede6 1841
af7b3903 1842 if(fCalorimeter=="EMCAL"){
8d230fa8 1843
1844 // Same sector
1845 Int_t j=0;
1846 for(Int_t i = 0; i < fNModules/2; i++){
1847 j=2*i;
91e1ea12 1848 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
8d230fa8 1849 }
1850
1851 // Same side
1852 for(Int_t i = 0; i < fNModules-2; i++){
91e1ea12 1853 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
8d230fa8 1854 }
1855 }//EMCAL
1856 else {//PHOS
91e1ea12 1857 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
1858 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
1859 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
8d230fa8 1860 }//PHOS
821c8090 1861 }
7e7694bb 1862
af7b3903 1863 //In case we want only pairs in same (super) module, check their origin.
1864 Bool_t ok = kTRUE;
1865 if(fSameSM && module1!=module2) ok=kFALSE;
1866 if(ok){
6175da48 1867
1868 //Check if one of the clusters comes from a conversion
d39cba7e 1869 if(fCheckConversion){
1870 if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
1871 else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
1872 }
6175da48 1873
0333ede6 1874 // Fill shower shape cut histograms
1875 if ( l01 > 0.01 && l01 < 0.4 &&
1876 l02 > 0.01 && l02 < 0.4 ) fhReSS[0]->Fill(pt,m); // Tight
1877 else if( l01 > 0.4 && l02 > 0.4 ) fhReSS[1]->Fill(pt,m); // Loose
1878 else if( l01 > 0.01 && l01 < 0.4 && l02 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
1879 else if( l02 > 0.01 && l02 < 0.4 && l01 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
1880
af7b3903 1881 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
5ae09196 1882 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
af7b3903 1883 if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
1884 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
1885 if(a < fAsymCuts[iasym]){
1886 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
6175da48 1887 //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym);
91e1ea12 1888 fhRe1 [index]->Fill(pt,m);
1889 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
6175da48 1890 if(fFillBadDistHisto){
1891 if(p1->DistToBad()>0 && p2->DistToBad()>0){
91e1ea12 1892 fhRe2 [index]->Fill(pt,m) ;
1893 if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
6175da48 1894 if(p1->DistToBad()>1 && p2->DistToBad()>1){
91e1ea12 1895 fhRe3 [index]->Fill(pt,m) ;
1896 if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
6175da48 1897 }// bad 3
1898 }// bad2
1899 }// Fill bad dist histos
1900 }//assymetry cut
1901 }// asymmetry cut loop
af7b3903 1902 }// bad 1
1903 }// pid bit loop
5ae09196 1904
af7b3903 1905 //Fill histograms with opening angle
91e1ea12 1906 fhRealOpeningAngle ->Fill(pt,angle);
1907 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
af7b3903 1908
1909 //Fill histograms with pair assymmetry
91e1ea12 1910 fhRePtAsym->Fill(pt,a);
1911 if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
1912 if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
af7b3903 1913
6175da48 1914 //-------------------------------------------------------
1915 //Get the number of cells needed for multi cut analysis.
1916 //-------------------------------------------------------
1917 Int_t ncell1 = 0;
1918 Int_t ncell2 = 0;
1919 if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim)){
0333ede6 1920
af7b3903 1921 AliVEvent * event = GetReader()->GetInputEvent();
1922 if(event){
1923 for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
1924 AliVCluster *cluster = event->GetCaloCluster(iclus);
5ae09196 1925
af7b3903 1926 Bool_t is = kFALSE;
1927 if (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
1928 else if(fCalorimeter == "PHOS" && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
5ae09196 1929
af7b3903 1930 if(is){
1931 if (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
1932 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
1933 } // PHOS or EMCAL cluster as requested in analysis
1934
1935 if(ncell2 > 0 && ncell1 > 0) break; // No need to continue the iteration
1936
1937 }
1938 //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
1939 }
6175da48 1940 }
1941
1942 //---------
1943 // MC data
1944 //---------
1945 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
99b8e903 1946 if(IsDataMC()) {
1947
1948 if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
1949 GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
1950 {
1951 fhReMCFromConversion->Fill(pt,m);
1952 }
1953 else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
1954 !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
1955 {
1956 fhReMCFromNotConversion->Fill(pt,m);
1957 }
1958 else
1959 {
1960 fhReMCFromMixConversion->Fill(pt,m);
1961 }
1962
1963 FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
1964 }
6175da48 1965
1966 //-----------------------
1967 //Multi cuts analysis
1968 //-----------------------
1969 if(fMultiCutAna){
1970 //Histograms for different PID bits selection
1971 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1972
1973 if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
91e1ea12 1974 p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ;
6175da48 1975
1976 //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
1977 } // pid bit cut loop
1978
1979 //Several pt,ncell and asymmetry cuts
af7b3903 1980 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1981 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1982 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1983 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
ba1eeb1f 1984 if(p1->E() > fPtCuts[ipt] && p2->E() > fPtCuts[ipt] &&
af7b3903 1985 a < fAsymCuts[iasym] &&
6175da48 1986 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
0333ede6 1987 fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
6175da48 1988 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
9302260a 1989 if(fFillSMCombinations && module1==module2){
0333ede6 1990 fhRePtNCellAsymCutsSM[module1][index]->Fill(pt,m) ;
6175da48 1991 }
1992 }
af7b3903 1993 }// pid bit cut loop
1994 }// icell loop
1995 }// pt cut loop
91e1ea12 1996 if(GetHistoTrackMultiplicityBins()){
1997 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
1998 if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
1999 }
2000 }
af7b3903 2001 }// multiple cuts analysis
2002 }// ok if same sm
7e7694bb 2003 }// second same event particle
2004 }// first cluster
0333ede6 2005
6175da48 2006 //-------------------------------------------------------------
2007 // Mixing
2008 //-------------------------------------------------------------
7e7694bb 2009 if(fDoOwnMix){
156549ae 2010 //printf("Cen bin %d, RP bin %d, e aver %f, mult %d\n",curCentrBin,curRPBin, eClusTot, nPhot);
6175da48 2011 //Recover events in with same characteristics as the current event
5025c139 2012 TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
7e7694bb 2013 Int_t nMixed = evMixList->GetSize() ;
2014 for(Int_t ii=0; ii<nMixed; ii++){
2015 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2016 Int_t nPhot2=ev2->GetEntriesFast() ;
2017 Double_t m = -999;
6175da48 2018 if(GetDebug() > 1)
2019 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, curCentrBin);
7e7694bb 2020
6175da48 2021 //---------------------------------
2022 //First loop on photons/clusters
2023 //---------------------------------
7e7694bb 2024 for(Int_t i1=0; i1<nPhot; i1++){
2025 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
6175da48 2026 if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2027
2028 //Get kinematics of cluster and (super) module of this cluster
7e7694bb 2029 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
af7b3903 2030 module1 = GetModuleNumber(p1);
91e1ea12 2031
6175da48 2032 //---------------------------------
2033 //First loop on photons/clusters
2034 //---------------------------------
7e7694bb 2035 for(Int_t i2=0; i2<nPhot2; i2++){
2036 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
91e1ea12 2037
6175da48 2038 //Get kinematics of second cluster and calculate those of the pair
7e7694bb 2039 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
6175da48 2040 m = (photon1+photon2).M() ;
7e7694bb 2041 Double_t pt = (photon1 + photon2).Pt();
2042 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2043
2044 //Check if opening angle is too large or too small compared to what is expected
2045 Double_t angle = photon1.Angle(photon2.Vect());
6175da48 2046 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
2047 if(GetDebug() > 2)
2048 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2049 continue;
2050 }
2051 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2052 if(GetDebug() > 2)
2053 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
2054 continue;
0333ede6 2055
6175da48 2056 }
7e7694bb 2057
2058 if(GetDebug() > 2)
2059 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
af7b3903 2060 p1->Pt(), p2->Pt(), pt,m,a);
6175da48 2061
af7b3903 2062 //In case we want only pairs in same (super) module, check their origin.
2063 module2 = GetModuleNumber(p2);
6175da48 2064
2065 //-------------------------------------------------------------------------------------------------
2066 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2067 //-------------------------------------------------------------------------------------------------
20218aea 2068 if(a < fAsymCuts[0] && fFillSMCombinations){
6175da48 2069 if(module1==module2 && module1 >=0 && module1<fNModules)
91e1ea12 2070 fhMiMod[module1]->Fill(pt,m) ;
0333ede6 2071
8d230fa8 2072 if(fCalorimeter=="EMCAL"){
2073
2074 // Same sector
2075 Int_t j=0;
2076 for(Int_t i = 0; i < fNModules/2; i++){
2077 j=2*i;
91e1ea12 2078 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
8d230fa8 2079 }
2080
2081 // Same side
2082 for(Int_t i = 0; i < fNModules-2; i++){
91e1ea12 2083 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
8d230fa8 2084 }
2085 }//EMCAL
2086 else {//PHOS
91e1ea12 2087 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2088 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2089 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
8d230fa8 2090 }//PHOS
2091
2092
6175da48 2093 }
2094
af7b3903 2095 Bool_t ok = kTRUE;
2096 if(fSameSM && module1!=module2) ok=kFALSE;
2097 if(ok){
6175da48 2098
2099 //Check if one of the clusters comes from a conversion
d39cba7e 2100 if(fCheckConversion){
2101 if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2102 else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2103 }
6175da48 2104 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
af7b3903 2105 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2106 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
2107 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
2108 if(a < fAsymCuts[iasym]){
2109 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
91e1ea12 2110 fhMi1 [index]->Fill(pt,m) ;
2111 if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
6175da48 2112 if(fFillBadDistHisto){
2113 if(p1->DistToBad()>0 && p2->DistToBad()>0){
2114 fhMi2 [index]->Fill(pt,m) ;
91e1ea12 2115 if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
6175da48 2116 if(p1->DistToBad()>1 && p2->DistToBad()>1){
2117 fhMi3 [index]->Fill(pt,m) ;
91e1ea12 2118 if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
6175da48 2119 }
af7b3903 2120 }
6175da48 2121 }// Fill bad dist histo
af7b3903 2122 }//Asymmetry cut
2123 }// Asymmetry loop
2124 }//PID cut
2125 }//loop for histograms
6175da48 2126
2127 //-----------------------
2128 //Multi cuts analysis
2129 //-----------------------
2130 if(fMultiCutAna){
2131 //Several pt,ncell and asymmetry cuts
0a14e9ae 2132
6175da48 2133 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2134 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2135 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2136 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2137 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
0a14e9ae 2138 a < fAsymCuts[iasym] //&&
2139 //p1->GetBtag() >= fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2140 ){
91e1ea12 2141 fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
6175da48 2142 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2143 }
2144 }// pid bit cut loop
2145 }// icell loop
2146 }// pt cut loop
2147 } // Multi cut ana
2148
2149 //Fill histograms with opening angle
91e1ea12 2150 fhMixedOpeningAngle ->Fill(pt,angle);
2151 fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
af7b3903 2152 }//ok
7e7694bb 2153 }// second cluster loop
2154 }//first cluster loop
2155 }//loop on mixed events
2156
6175da48 2157 //--------------------------------------------------------
2158 //Add the current event to the list of events for mixing
2159 //--------------------------------------------------------
7e7694bb 2160 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
af7b3903 2161 //Add current event to buffer and Remove redundant events
7e7694bb 2162 if(currentEvent->GetEntriesFast()>0){
2163 evMixList->AddFirst(currentEvent) ;
2164 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
72542aba 2165 if(evMixList->GetSize() >= GetNMaxEvMix())
7e7694bb 2166 {
2167 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2168 evMixList->RemoveLast() ;
2169 delete tmp ;
2170 }
2171 }
2172 else{ //empty event
2173 delete currentEvent ;
2174 currentEvent=0 ;
477d6cee 2175 }
7e7694bb 2176 }// DoOwnMix
c8fe2783 2177
1c5acb87 2178}
2179
6639984f 2180//____________________________________________________________________________________________________________________________________________________
c8fe2783 2181Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2182{
f8006433 2183 // retieves the event index and checks the vertex
2184 // in the mixed buffer returns -2 if vertex NOK
2185 // for normal events returns 0 if vertex OK and -1 if vertex NOK
2186
2187 Int_t evtIndex = -1 ;
2188 if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2189
2190 if (GetMixedEvent()){
2191
2192 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2193 GetVertex(vert,evtIndex);
2194
5025c139 2195 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 2196 evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2197 } else {// Single event
2198
2199 GetVertex(vert);
2200
5025c139 2201 if(TMath::Abs(vert[2])> GetZvertexCut())
f8006433 2202 evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2203 else
2204 evtIndex = 0 ;
c8fe2783 2205 }
0ae57829 2206 }//No MC reader
f8006433 2207 else {
2208 evtIndex = 0;
2209 vert[0] = 0. ;
2210 vert[1] = 0. ;
2211 vert[2] = 0. ;
2212 }
0ae57829 2213
f8006433 2214 return evtIndex ;
c8fe2783 2215}
f8006433 2216