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