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