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