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