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