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