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