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