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