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