1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //_________________________________________________________________________
17 // Class to collect two-photon invariant mass distributions for
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
22 //-- Author: Dmitri Peressounko (RRC "KI")
23 //-- Adapted to CaloTrackCorr frame by Lamia Benhabib (SUBATECH)
24 //-- and Gustavo Conesa (INFN-Frascati)
25 //_________________________________________________________________________
28 // --- ROOT system ---
31 //#include "Riostream.h"
35 #include "TClonesArray.h"
36 #include "TObjString.h"
37 #include "TDatabasePDG.h"
39 //---- AliRoot system ----
40 #include "AliAnaPi0.h"
41 #include "AliCaloTrackReader.h"
42 #include "AliCaloPID.h"
44 #include "AliFiducialCut.h"
45 #include "TParticle.h"
46 #include "AliVEvent.h"
47 #include "AliESDCaloCluster.h"
48 #include "AliESDEvent.h"
49 #include "AliAODEvent.h"
50 #include "AliNeutralMesonSelection.h"
51 #include "AliMixedEvent.h"
52 #include "AliAODMCParticle.h"
55 #include "AliPHOSGeoUtils.h"
56 #include "AliEMCALGeometry.h"
60 //______________________________________________________
61 AliAnaPi0::AliAnaPi0() : AliAnaCaloTrackCorrBaseClass(),
63 fCalorimeter(""), fNModules(22),
64 fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE), fAngleCut(0), fAngleMaxCut(7.),
65 fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
66 fNPtCuts(0), fNAsymCuts(0), fNCellNCuts(0), fNPIDBits(0),
67 fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
68 fFillSMCombinations(kFALSE), fCheckConversion(kFALSE),
69 fFillBadDistHisto(kFALSE), fFillSSCombinations(kFALSE),
70 fFillAngleHisto(kFALSE), fFillAsymmetryHisto(kFALSE), fFillOriginHisto(0), fFillArmenterosThetaStar(0),
71 fCheckAccInSector(kFALSE),
73 fhAverTotECluster(0), fhAverTotECell(0), fhAverTotECellvsCluster(0),
74 fhEDensityCluster(0), fhEDensityCell(0), fhEDensityCellvsCluster(0),
75 fhReMod(0x0), fhReSameSideEMCALMod(0x0), fhReSameSectorEMCALMod(0x0), fhReDiffPHOSMod(0x0),
76 fhMiMod(0x0), fhMiSameSideEMCALMod(0x0), fhMiSameSectorEMCALMod(0x0), fhMiDiffPHOSMod(0x0),
77 fhReConv(0x0), fhMiConv(0x0), fhReConv2(0x0), fhMiConv2(0x0),
78 fhRe1(0x0), fhMi1(0x0), fhRe2(0x0), fhMi2(0x0),
79 fhRe3(0x0), fhMi3(0x0), fhReInvPt1(0x0), fhMiInvPt1(0x0),
80 fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
81 fhRePtNCellAsymCuts(0x0), fhMiPtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM(),
82 fhRePIDBits(0x0), fhRePtMult(0x0), fhReSS(),
83 fhRePtAsym(0x0), fhRePtAsymPi0(0x0), fhRePtAsymEta(0x0),
84 fhEventBin(0), fhEventMixBin(0),
85 fhCentrality(0x0), fhCentralityNoPair(0x0),
86 fhEventPlaneResolution(0x0),
87 fhRealOpeningAngle(0x0), fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0), fhMixedCosOpeningAngle(0x0),
89 fhPrimPi0E(0x0), fhPrimPi0Pt(0x0), fhPrimPi0PtRejected(0x0),
90 fhPrimPi0AccE(0x0), fhPrimPi0AccPt(0x0),
91 fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0),
92 fhPrimPi0Yeta(0x0), fhPrimPi0YetaYcut(0x0), fhPrimPi0AccYeta(0x0),
93 fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
94 fhPrimPi0OpeningAngle(0x0), fhPrimPi0OpeningAngleAsym(0x0),fhPrimPi0CosOpeningAngle(0x0),
95 fhPrimPi0PtCentrality(0), fhPrimPi0PtEventPlane(0),
96 fhPrimPi0AccPtCentrality(0), fhPrimPi0AccPtEventPlane(0),
97 fhPrimEtaE(0x0), fhPrimEtaPt(0x0), fhPrimEtaPtRejected(0x0),
98 fhPrimEtaAccE(0x0), fhPrimEtaAccPt(0x0),
99 fhPrimEtaY(0x0), fhPrimEtaAccY(0x0),
100 fhPrimEtaYeta(0x0), fhPrimEtaYetaYcut(0x0), fhPrimEtaAccYeta(0x0),
101 fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
102 fhPrimEtaOpeningAngle(0x0), fhPrimEtaOpeningAngleAsym(0x0),fhPrimEtaCosOpeningAngle(0x0),
103 fhPrimEtaPtCentrality(0), fhPrimEtaPtEventPlane(0),
104 fhPrimEtaAccPtCentrality(0), fhPrimEtaAccPtEventPlane(0),
105 fhPrimPi0PtOrigin(0x0), fhPrimEtaPtOrigin(0x0),
106 fhMCOrgMass(), fhMCOrgAsym(), fhMCOrgDeltaEta(), fhMCOrgDeltaPhi(),
107 fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(),
108 fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec(),
109 fhMCPi0PtOrigin(0x0), fhMCEtaPtOrigin(0x0),
110 fhReMCFromConversion(0), fhReMCFromNotConversion(0), fhReMCFromMixConversion(0),
111 fhCosThStarPrimPi0(0), fhCosThStarPrimEta(0)//,
117 for(Int_t i = 0; i < 4; i++)
124 //_____________________
125 AliAnaPi0::~AliAnaPi0()
127 // Remove event containers
129 if(DoOwnMix() && fEventsList)
131 for(Int_t ic=0; ic<GetNCentrBin(); ic++)
133 for(Int_t iz=0; iz<GetNZvertBin(); iz++)
135 for(Int_t irp=0; irp<GetNRPBin(); irp++)
137 Int_t bin = GetEventMixBin(ic,iz,irp);
138 fEventsList[bin]->Delete() ;
139 delete fEventsList[bin] ;
143 delete[] fEventsList;
148 //______________________________
149 void AliAnaPi0::InitParameters()
151 //Init parameters when first called the analysis
152 //Set default parameters
153 SetInputAODName("PWG4Particle");
155 AddToHistogramsName("AnaPi0_");
157 fCalorimeter = "PHOS";
158 fUseAngleCut = kFALSE;
159 fUseAngleEDepCut = kFALSE;
161 fAngleMaxCut = TMath::Pi();
163 fMultiCutAna = kFALSE;
166 fPtCuts[0] = 0.; fPtCuts[1] = 0.3; fPtCuts[2] = 0.5;
167 for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
170 fAsymCuts[0] = 1.; fAsymCuts[1] = 0.7; //fAsymCuts[2] = 0.6; // fAsymCuts[3] = 0.1;
171 for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
174 fCellNCuts[0] = 0; fCellNCuts[1] = 1; fCellNCuts[2] = 2;
175 for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i] = 0;
178 fPIDBits[0] = 0; fPIDBits[1] = 2; // fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut, dispersion, neutral, dispersion&&neutral
179 for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
184 //_______________________________________
185 TObjString * AliAnaPi0::GetAnalysisCuts()
187 //Save parameters used for analysis
188 TString parList ; //this will be list of parameters used for this analysis.
189 const Int_t buffersize = 255;
190 char onePar[buffersize] ;
191 snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
193 snprintf(onePar,buffersize,"Number of bins in Centrality: %d \n",GetNCentrBin()) ;
195 snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
197 snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
199 snprintf(onePar,buffersize,"Depth of event buffer: %d \n",GetNMaxEvMix()) ;
201 snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f,\n",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
203 snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
204 for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
206 snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
207 for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
209 snprintf(onePar,buffersize,"Cuts: \n") ;
211 snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
213 snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
215 snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
218 snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
219 for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
221 snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
222 for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
226 return new TObjString(parList) ;
229 //_________________________________________
230 TList * AliAnaPi0::GetCreateOutputObjects()
232 // Create histograms to be saved in output file and
233 // store them in fOutputContainer
235 // Init the number of modules, set in the class AliCalorimeterUtils
236 fNModules = GetCaloUtils()->GetNumberOfSuperModulesUsed();
237 if(fCalorimeter=="PHOS" && fNModules > 4) fNModules = 4;
239 //create event containers
240 fEventsList = new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ;
242 for(Int_t ic=0; ic<GetNCentrBin(); ic++)
244 for(Int_t iz=0; iz<GetNZvertBin(); iz++)
246 for(Int_t irp=0; irp<GetNRPBin(); irp++)
248 Int_t bin = GetEventMixBin(ic,iz,irp);
249 fEventsList[bin] = new TList() ;
250 fEventsList[bin]->SetOwner(kFALSE);
255 TList * outputContainer = new TList() ;
256 outputContainer->SetName(GetName());
258 fhReMod = new TH2F*[fNModules] ;
259 fhMiMod = new TH2F*[fNModules] ;
261 if(fCalorimeter == "PHOS")
263 fhReDiffPHOSMod = new TH2F*[fNModules] ;
264 fhMiDiffPHOSMod = new TH2F*[fNModules] ;
268 fhReSameSectorEMCALMod = new TH2F*[fNModules/2] ;
269 fhReSameSideEMCALMod = new TH2F*[fNModules-2] ;
270 fhMiSameSectorEMCALMod = new TH2F*[fNModules/2] ;
271 fhMiSameSideEMCALMod = new TH2F*[fNModules-2] ;
275 fhRe1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
276 fhMi1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
277 if(fFillBadDistHisto)
279 fhRe2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
280 fhRe3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
281 fhMi2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
282 fhMi3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
286 fhReInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
287 fhMiInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
288 if(fFillBadDistHisto){
289 fhReInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
290 fhReInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
291 fhMiInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
292 fhMiInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
296 const Int_t buffersize = 255;
297 char key[buffersize] ;
298 char title[buffersize] ;
300 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins();
301 Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
302 Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
303 Float_t ptmax = GetHistogramRanges()->GetHistoPtMax();
304 Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
305 Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
306 Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
307 Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
308 Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
310 Int_t nmassbins = GetHistogramRanges()->GetHistoMassBins();
311 Int_t nasymbins = GetHistogramRanges()->GetHistoAsymmetryBins();
312 Float_t massmax = GetHistogramRanges()->GetHistoMassMax();
313 Float_t asymmax = GetHistogramRanges()->GetHistoAsymmetryMax();
314 Float_t massmin = GetHistogramRanges()->GetHistoMassMin();
315 Float_t asymmin = GetHistogramRanges()->GetHistoAsymmetryMin();
316 Int_t ntrmbins = GetHistogramRanges()->GetHistoTrackMultiplicityBins();
317 Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax();
318 Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
322 fhReConv = new TH2F("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
323 fhReConv->SetXTitle("p_{T} (GeV/c)");
324 fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
325 outputContainer->Add(fhReConv) ;
327 fhReConv2 = new TH2F("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
328 fhReConv2->SetXTitle("p_{T} (GeV/c)");
329 fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
330 outputContainer->Add(fhReConv2) ;
334 fhMiConv = new TH2F("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
335 fhMiConv->SetXTitle("p_{T} (GeV/c)");
336 fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
337 outputContainer->Add(fhMiConv) ;
339 fhMiConv2 = new TH2F("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
340 fhMiConv2->SetXTitle("p_{T} (GeV/c)");
341 fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
342 outputContainer->Add(fhMiConv2) ;
346 for(Int_t ic=0; ic<GetNCentrBin(); ic++)
348 for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
350 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
352 Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
353 //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
354 //Distance to bad module 1
355 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
356 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
357 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
358 fhRe1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
359 fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
360 fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
361 //printf("name: %s\n ",fhRe1[index]->GetName());
362 outputContainer->Add(fhRe1[index]) ;
364 if(fFillBadDistHisto)
366 //Distance to bad module 2
367 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
368 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
369 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
370 fhRe2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
371 fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
372 fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
373 outputContainer->Add(fhRe2[index]) ;
375 //Distance to bad module 3
376 snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
377 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
378 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
379 fhRe3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
380 fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
381 fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
382 outputContainer->Add(fhRe3[index]) ;
388 //Distance to bad module 1
389 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
390 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
391 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
392 fhReInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
393 fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
394 fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
395 outputContainer->Add(fhReInvPt1[index]) ;
397 if(fFillBadDistHisto){
398 //Distance to bad module 2
399 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
400 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
401 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
402 fhReInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
403 fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
404 fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
405 outputContainer->Add(fhReInvPt2[index]) ;
407 //Distance to bad module 3
408 snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
409 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
410 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
411 fhReInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
412 fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
413 fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
414 outputContainer->Add(fhReInvPt3[index]) ;
420 //Distance to bad module 1
421 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
422 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
423 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
424 fhMi1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
425 fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
426 fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
427 outputContainer->Add(fhMi1[index]) ;
428 if(fFillBadDistHisto){
429 //Distance to bad module 2
430 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
431 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
432 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
433 fhMi2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
434 fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
435 fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
436 outputContainer->Add(fhMi2[index]) ;
438 //Distance to bad module 3
439 snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
440 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
441 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
442 fhMi3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
443 fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
444 fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
445 outputContainer->Add(fhMi3[index]) ;
451 //Distance to bad module 1
452 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
453 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
454 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
455 fhMiInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
456 fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
457 fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
458 outputContainer->Add(fhMiInvPt1[index]) ;
459 if(fFillBadDistHisto){
460 //Distance to bad module 2
461 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
462 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
463 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
464 fhMiInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
465 fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
466 fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
467 outputContainer->Add(fhMiInvPt2[index]) ;
469 //Distance to bad module 3
470 snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
471 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
472 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
473 fhMiInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
474 fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
475 fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
476 outputContainer->Add(fhMiInvPt3[index]) ;
484 if(fFillAsymmetryHisto)
486 fhRePtAsym = new TH2F("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
487 fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
488 fhRePtAsym->SetYTitle("Asymmetry");
489 outputContainer->Add(fhRePtAsym);
491 fhRePtAsymPi0 = new TH2F("hRePtAsymPi0","Asymmetry vs pt, for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
492 fhRePtAsymPi0->SetXTitle("p_{T} (GeV/c)");
493 fhRePtAsymPi0->SetYTitle("Asymmetry");
494 outputContainer->Add(fhRePtAsymPi0);
496 fhRePtAsymEta = new TH2F("hRePtAsymEta","Asymmetry vs pt, for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
497 fhRePtAsymEta->SetXTitle("p_{T} (GeV/c)");
498 fhRePtAsymEta->SetYTitle("Asymmetry");
499 outputContainer->Add(fhRePtAsymEta);
504 fhRePIDBits = new TH2F*[fNPIDBits];
505 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
506 snprintf(key, buffersize,"hRe_pidbit%d",ipid) ;
507 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
508 fhRePIDBits[ipid] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
509 fhRePIDBits[ipid]->SetXTitle("p_{T} (GeV/c)");
510 fhRePIDBits[ipid]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
511 outputContainer->Add(fhRePIDBits[ipid]) ;
514 fhRePtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
515 fhMiPtNCellAsymCuts = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
517 if(fFillSMCombinations)
518 for(Int_t iSM = 0; iSM < fNModules; iSM++) fhRePtNCellAsymCutsSM[iSM] = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
521 for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
523 for(Int_t icell=0; icell<fNCellNCuts; icell++)
525 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
527 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
528 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
529 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
530 //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
531 fhRePtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
532 fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
533 fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
534 outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
536 snprintf(key, buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
537 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
538 fhMiPtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
539 fhMiPtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
540 fhMiPtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
541 outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
543 if(fFillSMCombinations)
545 for(Int_t iSM = 0; iSM < fNModules; iSM++)
547 snprintf(key, buffersize,"hRe_pt%d_cell%d_asym%d_SM%d",ipt,icell,iasym,iSM) ;
548 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM %d ",
549 fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym],iSM) ;
550 fhRePtNCellAsymCutsSM[iSM][index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
551 fhRePtNCellAsymCutsSM[iSM][index]->SetXTitle("p_{T} (GeV/c)");
552 fhRePtNCellAsymCutsSM[iSM][index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
553 outputContainer->Add(fhRePtNCellAsymCutsSM[iSM][index]) ;
563 fhRePtMult = new TH3F*[fNAsymCuts] ;
564 for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++)
566 fhRePtMult[iasym] = new TH3F(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
567 nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
568 fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
569 fhRePtMult[iasym]->SetYTitle("Track multiplicity");
570 fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
571 outputContainer->Add(fhRePtMult[iasym]) ;
574 }// multi cuts analysis
576 if(fFillSSCombinations)
579 fhReSS[0] = new TH2F("hRe_SS_Tight"," 0.01 < #lambda_{0}^{2} < 0.4",
580 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
581 fhReSS[0]->SetXTitle("p_{T} (GeV/c)");
582 fhReSS[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
583 outputContainer->Add(fhReSS[0]) ;
586 fhReSS[1] = new TH2F("hRe_SS_Loose"," #lambda_{0}^{2} > 0.4",
587 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
588 fhReSS[1]->SetXTitle("p_{T} (GeV/c)");
589 fhReSS[1]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
590 outputContainer->Add(fhReSS[1]) ;
593 fhReSS[2] = new TH2F("hRe_SS_Both"," cluster_{1} #lambda_{0}^{2} > 0.4; cluster_{2} 0.01 < #lambda_{0}^{2} < 0.4",
594 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
595 fhReSS[2]->SetXTitle("p_{T} (GeV/c)");
596 fhReSS[2]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
597 outputContainer->Add(fhReSS[2]) ;
602 fhEventBin=new TH1I("hEventBin","Number of real pairs per bin(cen,vz,rp)",
603 GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
604 GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
605 fhEventBin->SetXTitle("bin");
606 outputContainer->Add(fhEventBin) ;
609 fhEventMixBin=new TH1I("hEventMixBin","Number of mixed pairs per bin(cen,vz,rp)",
610 GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
611 GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
612 fhEventMixBin->SetXTitle("bin");
613 outputContainer->Add(fhEventMixBin) ;
618 fhCentrality=new TH1F("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
619 fhCentrality->SetXTitle("Centrality bin");
620 outputContainer->Add(fhCentrality) ;
622 fhCentralityNoPair=new TH1F("hCentralityBinNoPair","Number of events in centrality bin, with no cluster pairs",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
623 fhCentralityNoPair->SetXTitle("Centrality bin");
624 outputContainer->Add(fhCentralityNoPair) ;
627 if(GetNRPBin() > 1 && GetNCentrBin()>1 )
629 fhEventPlaneResolution=new TH2F("hEventPlaneResolution","Event plane resolution",GetNCentrBin(),0,GetNCentrBin(),100,0.,TMath::TwoPi()) ;
630 fhEventPlaneResolution->SetYTitle("Resolution");
631 fhEventPlaneResolution->SetXTitle("Centrality Bin");
632 outputContainer->Add(fhEventPlaneResolution) ;
637 fhRealOpeningAngle = new TH2F
638 ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
639 fhRealOpeningAngle->SetYTitle("#theta(rad)");
640 fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
641 outputContainer->Add(fhRealOpeningAngle) ;
643 fhRealCosOpeningAngle = new TH2F
644 ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
645 fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
646 fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
647 outputContainer->Add(fhRealCosOpeningAngle) ;
651 fhMixedOpeningAngle = new TH2F
652 ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
653 fhMixedOpeningAngle->SetYTitle("#theta(rad)");
654 fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
655 outputContainer->Add(fhMixedOpeningAngle) ;
657 fhMixedCosOpeningAngle = new TH2F
658 ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
659 fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
660 fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
661 outputContainer->Add(fhMixedCosOpeningAngle) ;
666 //Histograms filled only if MC data is requested
669 fhReMCFromConversion = new TH2F("hReMCFromConversion","Invariant mass of 2 clusters originated in conversions",
670 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
671 fhReMCFromConversion->SetXTitle("p_{T} (GeV/c)");
672 fhReMCFromConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
673 outputContainer->Add(fhReMCFromConversion) ;
675 fhReMCFromNotConversion = new TH2F("hReMCNotFromConversion","Invariant mass of 2 clusters not originated in conversions",
676 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
677 fhReMCFromNotConversion->SetXTitle("p_{T} (GeV/c)");
678 fhReMCFromNotConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
679 outputContainer->Add(fhReMCFromNotConversion) ;
681 fhReMCFromMixConversion = new TH2F("hReMCFromMixConversion","Invariant mass of 2 clusters one from conversion and the other not",
682 nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
683 fhReMCFromMixConversion->SetXTitle("p_{T} (GeV/c)");
684 fhReMCFromMixConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
685 outputContainer->Add(fhReMCFromMixConversion) ;
689 fhPrimPi0E = new TH1F("hPrimPi0E","Primary pi0 E, Y<1",nptbins,ptmin,ptmax) ;
690 fhPrimPi0AccE = new TH1F("hPrimPi0AccE","Primary pi0 E with both photons in acceptance",nptbins,ptmin,ptmax) ;
691 fhPrimPi0E ->SetXTitle("E (GeV)");
692 fhPrimPi0AccE->SetXTitle("E (GeV)");
693 outputContainer->Add(fhPrimPi0E) ;
694 outputContainer->Add(fhPrimPi0AccE) ;
696 fhPrimPi0PtRejected = new TH1F("hPrimPi0PtRejected","Primary pi0 pt",nptbins,ptmin,ptmax) ;
697 fhPrimPi0Pt = new TH1F("hPrimPi0Pt","Primary pi0 pt, Y<1",nptbins,ptmin,ptmax) ;
698 fhPrimPi0AccPt = new TH1F("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
699 fhPrimPi0Pt ->SetXTitle("p_{T} (GeV/c)");
700 fhPrimPi0AccPt->SetXTitle("p_{T} (GeV/c)");
701 outputContainer->Add(fhPrimPi0PtRejected) ;
702 outputContainer->Add(fhPrimPi0Pt) ;
703 outputContainer->Add(fhPrimPi0AccPt) ;
705 Int_t netabinsopen = TMath::Nint(netabins*4/(etamax-etamin));
706 fhPrimPi0Y = new TH2F("hPrimPi0Rapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
707 fhPrimPi0Y ->SetYTitle("Rapidity");
708 fhPrimPi0Y ->SetXTitle("p_{T} (GeV/c)");
709 outputContainer->Add(fhPrimPi0Y) ;
711 fhPrimPi0Yeta = new TH2F("hPrimPi0PseudoRapidity","PseudoRapidity of primary pi0",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
712 fhPrimPi0Yeta ->SetYTitle("#eta");
713 fhPrimPi0Yeta ->SetXTitle("p_{T} (GeV/c)");
714 outputContainer->Add(fhPrimPi0Yeta) ;
716 fhPrimPi0YetaYcut = new TH2F("hPrimPi0PseudoRapidityYcut","PseudoRapidity of primary pi0, |Y|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
717 fhPrimPi0YetaYcut ->SetYTitle("#eta");
718 fhPrimPi0YetaYcut ->SetXTitle("p_{T} (GeV/c)");
719 outputContainer->Add(fhPrimPi0YetaYcut) ;
721 fhPrimPi0AccY = new TH2F("hPrimPi0AccRapidity","Rapidity of primary pi0 with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
722 fhPrimPi0AccY->SetYTitle("Rapidity");
723 fhPrimPi0AccY->SetXTitle("p_{T} (GeV/c)");
724 outputContainer->Add(fhPrimPi0AccY) ;
726 fhPrimPi0AccYeta = new TH2F("hPrimPi0AccPseudoRapidity","PseudoRapidity of primary pi0 with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
727 fhPrimPi0AccYeta ->SetYTitle("#eta");
728 fhPrimPi0AccYeta ->SetXTitle("p_{T} (GeV/c)");
729 outputContainer->Add(fhPrimPi0AccYeta) ;
731 Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin));
732 fhPrimPi0Phi = new TH2F("hPrimPi0Phi","Azimuthal of primary pi0, Y<1",nptbins,ptmin,ptmax,nphibinsopen,0,360) ;
733 fhPrimPi0Phi->SetYTitle("#phi (deg)");
734 fhPrimPi0Phi->SetXTitle("p_{T} (GeV/c)");
735 outputContainer->Add(fhPrimPi0Phi) ;
737 fhPrimPi0AccPhi = new TH2F("hPrimPi0AccPhi","Azimuthal of primary pi0 with accepted daughters",nptbins,ptmin,ptmax,
738 nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
739 fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
740 fhPrimPi0AccPhi->SetXTitle("p_{T} (GeV/c)");
741 outputContainer->Add(fhPrimPi0AccPhi) ;
743 fhPrimPi0PtCentrality = new TH2F("hPrimPi0PtCentrality","Primary pi0 pt vs reco centrality, Y<1",nptbins,ptmin,ptmax, 100, 0, 100) ;
744 fhPrimPi0AccPtCentrality = new TH2F("hPrimPi0AccPtCentrality","Primary pi0 with both photons in acceptance pt vs reco centrality",nptbins,ptmin,ptmax, 100, 0, 100) ;
745 fhPrimPi0PtCentrality ->SetXTitle("p_{T} (GeV/c)");
746 fhPrimPi0AccPtCentrality->SetXTitle("p_{T} (GeV/c)");
747 fhPrimPi0PtCentrality ->SetYTitle("Centrality");
748 fhPrimPi0AccPtCentrality->SetYTitle("Centrality");
749 outputContainer->Add(fhPrimPi0PtCentrality) ;
750 outputContainer->Add(fhPrimPi0AccPtCentrality) ;
752 fhPrimPi0PtEventPlane = new TH2F("hPrimPi0PtEventPlane","Primary pi0 pt vs reco event plane angle, Y<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
753 fhPrimPi0AccPtEventPlane = new TH2F("hPrimPi0AccPtEventPlane","Primary pi0 with both photons in acceptance pt vs reco event plane angle",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
754 fhPrimPi0PtEventPlane ->SetXTitle("p_{T} (GeV/c)");
755 fhPrimPi0AccPtEventPlane->SetXTitle("p_{T} (GeV/c)");
756 fhPrimPi0PtEventPlane ->SetYTitle("Event Plane Angle (rad)");
757 fhPrimPi0AccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
758 outputContainer->Add(fhPrimPi0PtEventPlane) ;
759 outputContainer->Add(fhPrimPi0AccPtEventPlane) ;
763 fhPrimEtaE = new TH1F("hPrimEtaE","Primary eta E",nptbins,ptmin,ptmax) ;
764 fhPrimEtaAccE = new TH1F("hPrimEtaAccE","Primary eta E with both photons in acceptance",nptbins,ptmin,ptmax) ;
765 fhPrimEtaE ->SetXTitle("E (GeV)");
766 fhPrimEtaAccE->SetXTitle("E (GeV)");
767 outputContainer->Add(fhPrimEtaE) ;
768 outputContainer->Add(fhPrimEtaAccE) ;
770 fhPrimEtaPtRejected = new TH1F("hPrimEtaPtRejected","Primary eta pt",nptbins,ptmin,ptmax) ;
771 fhPrimEtaPt = new TH1F("hPrimEtaPt","Primary eta pt",nptbins,ptmin,ptmax) ;
772 fhPrimEtaAccPt = new TH1F("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
773 fhPrimEtaPt ->SetXTitle("p_{T} (GeV/c)");
774 fhPrimEtaAccPt->SetXTitle("p_{T} (GeV/c)");
775 outputContainer->Add(fhPrimEtaPtRejected) ;
776 outputContainer->Add(fhPrimEtaPt) ;
777 outputContainer->Add(fhPrimEtaAccPt) ;
779 fhPrimEtaY = new TH2F("hPrimEtaRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
780 fhPrimEtaY->SetYTitle("Rapidity");
781 fhPrimEtaY->SetXTitle("p_{T} (GeV/c)");
782 outputContainer->Add(fhPrimEtaY) ;
784 fhPrimEtaYeta = new TH2F("hPrimEtaPseudoRapidityEta","PsuedoRapidity of primary eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
785 fhPrimEtaYeta->SetYTitle("Rapidity");
786 fhPrimEtaYeta->SetXTitle("p_{T} (GeV/c)");
787 outputContainer->Add(fhPrimEtaYeta) ;
789 fhPrimEtaYetaYcut = new TH2F("hPrimEtaPseudoRapidityEtaYcut","PseudoRapidity of primary eta, |Y|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
790 fhPrimEtaYetaYcut->SetYTitle("PseudoRapidity");
791 fhPrimEtaYetaYcut->SetXTitle("p_{T} (GeV/c)");
792 outputContainer->Add(fhPrimEtaYetaYcut) ;
794 fhPrimEtaAccY = new TH2F("hPrimEtaAccRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
795 fhPrimEtaAccY->SetYTitle("Rapidity");
796 fhPrimEtaAccY->SetXTitle("p_{T} (GeV/c)");
797 outputContainer->Add(fhPrimEtaAccY) ;
799 fhPrimEtaAccYeta = new TH2F("hPrimEtaAccPseudoRapidity","PseudoRapidity of primary eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
800 fhPrimEtaAccYeta->SetYTitle("PSeudoRapidity");
801 fhPrimEtaAccYeta->SetXTitle("p_{T} (GeV/c)");
802 outputContainer->Add(fhPrimEtaAccYeta) ;
804 fhPrimEtaPhi = new TH2F("hPrimEtaPhi","Azimuthal of primary eta",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
805 fhPrimEtaPhi->SetYTitle("#phi (deg)");
806 fhPrimEtaPhi->SetXTitle("p_{T} (GeV/c)");
807 outputContainer->Add(fhPrimEtaPhi) ;
809 fhPrimEtaAccPhi = new TH2F("hPrimEtaAccPhi","Azimuthal of primary eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
810 fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
811 fhPrimEtaAccPhi->SetXTitle("p_{T} (GeV/c)");
812 outputContainer->Add(fhPrimEtaAccPhi) ;
814 fhPrimEtaPtCentrality = new TH2F("hPrimEtaPtCentrality","Primary eta pt vs reco centrality, Y<1",nptbins,ptmin,ptmax, 100, 0, 100) ;
815 fhPrimEtaAccPtCentrality = new TH2F("hPrimEtaAccPtCentrality","Primary eta with both photons in acceptance pt vs reco centrality",nptbins,ptmin,ptmax, 100, 0, 100) ;
816 fhPrimEtaPtCentrality ->SetXTitle("p_{T} (GeV/c)");
817 fhPrimEtaAccPtCentrality->SetXTitle("p_{T} (GeV/c)");
818 fhPrimEtaPtCentrality ->SetYTitle("Centrality");
819 fhPrimEtaAccPtCentrality->SetYTitle("Centrality");
820 outputContainer->Add(fhPrimEtaPtCentrality) ;
821 outputContainer->Add(fhPrimEtaAccPtCentrality) ;
823 fhPrimEtaPtEventPlane = new TH2F("hPrimEtaPtEventPlane","Primary eta pt vs reco event plane angle, Y<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
824 fhPrimEtaAccPtEventPlane = new TH2F("hPrimEtaAccPtEventPlane","Primary eta with both photons in acceptance pt vs reco event plane angle",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
825 fhPrimEtaPtEventPlane ->SetXTitle("p_{T} (GeV/c)");
826 fhPrimEtaAccPtEventPlane->SetXTitle("p_{T} (GeV/c)");
827 fhPrimEtaPtEventPlane ->SetYTitle("Event Plane Angle (rad)");
828 fhPrimEtaAccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
829 outputContainer->Add(fhPrimEtaPtEventPlane) ;
830 outputContainer->Add(fhPrimEtaAccPtEventPlane) ;
834 fhPrimPi0OpeningAngle = new TH2F
835 ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
836 fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
837 fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
838 outputContainer->Add(fhPrimPi0OpeningAngle) ;
840 fhPrimPi0OpeningAngleAsym = new TH2F
841 ("hPrimPi0OpeningAngleAsym","Angle between all primary #gamma pair vs Asymmetry, p_{T}>5 GeV/c",100,0,1,100,0,0.5);
842 fhPrimPi0OpeningAngleAsym->SetXTitle("|A|=| (E_{1}-E_{2}) / (E_{1}+E_{2}) |");
843 fhPrimPi0OpeningAngleAsym->SetYTitle("#theta(rad)");
844 outputContainer->Add(fhPrimPi0OpeningAngleAsym) ;
846 fhPrimPi0CosOpeningAngle = new TH2F
847 ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
848 fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
849 fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
850 outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
852 fhPrimEtaOpeningAngle = new TH2F
853 ("hPrimEtaOpeningAngle","Angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,0,0.5);
854 fhPrimEtaOpeningAngle->SetYTitle("#theta(rad)");
855 fhPrimEtaOpeningAngle->SetXTitle("E_{#eta} (GeV)");
856 outputContainer->Add(fhPrimEtaOpeningAngle) ;
858 fhPrimEtaOpeningAngleAsym = new TH2F
859 ("hPrimEtaOpeningAngleAsym","Angle between all primary #gamma pair vs Asymmetry, p_{T}>5 GeV/c",100,0,1,100,0,0.5);
860 fhPrimEtaOpeningAngleAsym->SetXTitle("|A|=| (E_{1}-E_{2}) / (E_{1}+E_{2}) |");
861 fhPrimEtaOpeningAngleAsym->SetYTitle("#theta(rad)");
862 outputContainer->Add(fhPrimEtaOpeningAngleAsym) ;
865 fhPrimEtaCosOpeningAngle = new TH2F
866 ("hPrimEtaCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,-1,1);
867 fhPrimEtaCosOpeningAngle->SetYTitle("cos (#theta) ");
868 fhPrimEtaCosOpeningAngle->SetXTitle("E_{ #eta} (GeV)");
869 outputContainer->Add(fhPrimEtaCosOpeningAngle) ;
878 fhPrimPi0PtOrigin = new TH2F("hPrimPi0PtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
879 fhPrimPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
880 fhPrimPi0PtOrigin->SetYTitle("Origin");
881 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
882 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
883 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
884 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
885 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
886 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
887 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
888 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
889 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
890 fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
891 outputContainer->Add(fhPrimPi0PtOrigin) ;
893 fhMCPi0PtOrigin = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
894 fhMCPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
895 fhMCPi0PtOrigin->SetYTitle("Origin");
896 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
897 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
898 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
899 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
900 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
901 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
902 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
903 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
904 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
905 fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
906 outputContainer->Add(fhMCPi0PtOrigin) ;
909 fhPrimEtaPtOrigin = new TH2F("hPrimEtaPtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
910 fhPrimEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
911 fhPrimEtaPtOrigin->SetYTitle("Origin");
912 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
913 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
914 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
915 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
916 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
917 fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
919 outputContainer->Add(fhPrimEtaPtOrigin) ;
921 fhMCEtaPtOrigin = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
922 fhMCEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
923 fhMCEtaPtOrigin->SetYTitle("Origin");
924 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
925 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
926 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
927 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
928 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
929 fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
931 outputContainer->Add(fhMCEtaPtOrigin) ;
933 for(Int_t i = 0; i<13; i++)
935 fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
936 fhMCOrgMass[i]->SetXTitle("p_{T} (GeV/c)");
937 fhMCOrgMass[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
938 outputContainer->Add(fhMCOrgMass[i]) ;
940 fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
941 fhMCOrgAsym[i]->SetXTitle("p_{T} (GeV/c)");
942 fhMCOrgAsym[i]->SetYTitle("A");
943 outputContainer->Add(fhMCOrgAsym[i]) ;
945 fhMCOrgDeltaEta[i] = new TH2F(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs pt, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
946 fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
947 fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
948 outputContainer->Add(fhMCOrgDeltaEta[i]) ;
950 fhMCOrgDeltaPhi[i]= new TH2F(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs p_{T}, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
951 fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
952 fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
953 outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
959 fhMCPi0MassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
960 fhMCPi0MassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
961 fhMCPi0PtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
962 fhMCEtaMassPtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
963 fhMCEtaMassPtTrue = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
964 fhMCEtaPtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
965 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
966 for(Int_t icell=0; icell<fNCellNCuts; icell++){
967 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
968 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
970 fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
971 Form("Reconstructed Mass vs reconstructed p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
972 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
973 fhMCPi0MassPtRec[index]->SetXTitle("p_{T, reconstructed} (GeV/c)");
974 fhMCPi0MassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
975 outputContainer->Add(fhMCPi0MassPtRec[index]) ;
977 fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
978 Form("Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
979 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
980 fhMCPi0MassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
981 fhMCPi0MassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
982 outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
984 fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
985 Form("Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
986 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
987 fhMCPi0PtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
988 fhMCPi0PtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
989 outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
991 fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
992 Form("Reconstructed Mass vs reconstructed p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
993 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
994 fhMCEtaMassPtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
995 fhMCEtaMassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
996 outputContainer->Add(fhMCEtaMassPtRec[index]) ;
998 fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
999 Form("Reconstructed Mass vs generated p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1000 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1001 fhMCEtaMassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
1002 fhMCEtaMassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1003 outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
1005 fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1006 Form("Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1007 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1008 fhMCEtaPtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
1009 fhMCEtaPtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
1010 outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
1017 fhMCPi0MassPtTrue = new TH2F*[1];
1018 fhMCPi0PtTruePtRec = new TH2F*[1];
1019 fhMCEtaMassPtTrue = new TH2F*[1];
1020 fhMCEtaPtTruePtRec = new TH2F*[1];
1022 fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1023 fhMCPi0MassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
1024 fhMCPi0MassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1025 outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
1027 fhMCPi0PtTruePtRec[0]= new TH2F("hMCPi0PtTruePtRec","Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1028 fhMCPi0PtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
1029 fhMCPi0PtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
1030 outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
1032 fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1033 fhMCEtaMassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
1034 fhMCEtaMassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1035 outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
1037 fhMCEtaPtTruePtRec[0]= new TH2F("hMCEtaPtTruePtRec","Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1038 fhMCEtaPtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
1039 fhMCEtaPtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
1040 outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
1045 if(fFillSMCombinations)
1047 TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
1048 for(Int_t imod=0; imod<fNModules; imod++)
1050 //Module dependent invariant mass
1051 snprintf(key, buffersize,"hReMod_%d",imod) ;
1052 snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
1053 fhReMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1054 fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
1055 fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1056 outputContainer->Add(fhReMod[imod]) ;
1057 if(fCalorimeter=="PHOS")
1059 snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
1060 snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1061 fhReDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1062 fhReDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
1063 fhReDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1064 outputContainer->Add(fhReDiffPHOSMod[imod]) ;
1067 if(imod<fNModules/2)
1069 snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
1070 snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1071 fhReSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1072 fhReSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
1073 fhReSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1074 outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
1076 if(imod<fNModules-2)
1078 snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
1079 snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1080 fhReSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1081 fhReSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
1082 fhReSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1083 outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
1089 snprintf(key, buffersize,"hMiMod_%d",imod) ;
1090 snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
1091 fhMiMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1092 fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
1093 fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1094 outputContainer->Add(fhMiMod[imod]) ;
1096 if(fCalorimeter=="PHOS"){
1097 snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
1098 snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1099 fhMiDiffPHOSMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1100 fhMiDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
1101 fhMiDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1102 outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
1105 if(imod<fNModules/2)
1107 snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
1108 snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1109 fhMiSameSectorEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1110 fhMiSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
1111 fhMiSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1112 outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
1114 if(imod<fNModules-2){
1116 snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
1117 snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1118 fhMiSameSideEMCALMod[imod] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1119 fhMiSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
1120 fhMiSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
1121 outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
1125 }//loop combinations
1126 } // SM combinations
1128 if(fFillArmenterosThetaStar && IsDataMC())
1130 TString ebin[] = {"8 < E < 12 GeV","12 < E < 16 GeV", "16 < E < 20 GeV", "E > 20 GeV" };
1131 Int_t narmbins = 400;
1133 Float_t armmax = 0.4;
1135 for(Int_t i = 0; i < 4; i++)
1137 fhArmPrimPi0[i] = new TH2F(Form("hArmenterosPrimPi0EBin%d",i),
1138 Form("Armenteros of primary #pi^{0}, %s",ebin[i].Data()),
1139 200, -1, 1, narmbins,armmin,armmax);
1140 fhArmPrimPi0[i]->SetYTitle("p_{T}^{Arm}");
1141 fhArmPrimPi0[i]->SetXTitle("#alpha^{Arm}");
1142 outputContainer->Add(fhArmPrimPi0[i]) ;
1144 fhArmPrimEta[i] = new TH2F(Form("hArmenterosPrimEtaEBin%d",i),
1145 Form("Armenteros of primary #eta, %s",ebin[i].Data()),
1146 200, -1, 1, narmbins,armmin,armmax);
1147 fhArmPrimEta[i]->SetYTitle("p_{T}^{Arm}");
1148 fhArmPrimEta[i]->SetXTitle("#alpha^{Arm}");
1149 outputContainer->Add(fhArmPrimEta[i]) ;
1153 // Same as asymmetry ...
1154 fhCosThStarPrimPi0 = new TH2F
1155 ("hCosThStarPrimPi0","cos(#theta *) for primary #pi^{0}",nptbins,ptmin,ptmax,200,-1,1);
1156 fhCosThStarPrimPi0->SetYTitle("cos(#theta *)");
1157 fhCosThStarPrimPi0->SetXTitle("E_{ #pi^{0}} (GeV)");
1158 outputContainer->Add(fhCosThStarPrimPi0) ;
1160 fhCosThStarPrimEta = new TH2F
1161 ("hCosThStarPrimEta","cos(#theta *) for primary #eta",nptbins,ptmin,ptmax,200,-1,1);
1162 fhCosThStarPrimEta->SetYTitle("cos(#theta *)");
1163 fhCosThStarPrimEta->SetXTitle("E_{ #eta} (GeV)");
1164 outputContainer->Add(fhCosThStarPrimEta) ;
1168 // for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
1170 // printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
1174 return outputContainer;
1177 //___________________________________________________
1178 void AliAnaPi0::Print(const Option_t * /*opt*/) const
1180 //Print some relevant parameters set for the analysis
1181 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1182 AliAnaCaloTrackCorrBaseClass::Print(" ");
1184 printf("Number of bins in Centrality: %d \n",GetNCentrBin()) ;
1185 printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
1186 printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
1187 printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
1188 printf("Pair in same Module: %d \n",fSameSM) ;
1189 printf("Cuts: \n") ;
1190 // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
1191 printf("Number of modules: %d \n",fNModules) ;
1192 printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
1193 printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
1194 printf("\tasymmetry < ");
1195 for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
1198 printf("PID selection bits: n = %d, \n",fNPIDBits) ;
1199 printf("\tPID bit = ");
1200 for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
1204 printf("pT cuts: n = %d, \n",fNPtCuts) ;
1206 for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
1209 printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
1210 printf("\tnCell > ");
1211 for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
1215 printf("------------------------------------------------------\n") ;
1218 //________________________________________
1219 void AliAnaPi0::FillAcceptanceHistograms()
1221 //Fill acceptance histograms if MC data is available
1223 Float_t cen = GetEventCentrality();
1224 Float_t ep = GetEventPlaneAngle();
1226 if(GetReader()->ReadStack())
1228 AliStack * stack = GetMCStack();
1231 for(Int_t i=0 ; i<stack->GetNtrack(); i++)
1233 if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
1235 TParticle * prim = stack->Particle(i) ;
1236 Int_t pdg = prim->GetPdgCode();
1237 //printf("i %d, %s %d %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
1238 // prim->GetName(), prim->GetPdgCode());
1240 if( pdg == 111 || pdg == 221)
1242 Double_t pi0Pt = prim->Pt() ;
1243 Double_t pi0E = prim->Energy() ;
1244 if(pi0E == TMath::Abs(prim->Pz()))
1246 if( pdg == 111 ) fhPrimPi0PtRejected->Fill(pi0Pt);
1247 else fhPrimEtaPtRejected->Fill(pi0Pt);
1248 continue ; //Protection against floating point exception
1251 Double_t pi0Y = 0.5*TMath::Log((pi0E-prim->Pz())/(pi0E+prim->Pz())) ;
1252 Double_t pi0Yeta = prim->Eta() ;
1253 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
1257 if(TMath::Abs(pi0Y) < 1.0)
1259 fhPrimPi0E ->Fill(pi0E ) ;
1260 fhPrimPi0Pt ->Fill(pi0Pt) ;
1261 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
1262 fhPrimPi0YetaYcut ->Fill(pi0Pt,pi0Yeta) ;
1263 fhPrimPi0PtCentrality->Fill(pi0Pt,cen) ;
1264 fhPrimPi0PtEventPlane->Fill(pi0Pt,ep ) ;
1266 fhPrimPi0Y ->Fill(pi0Pt, pi0Y) ;
1267 fhPrimPi0Yeta->Fill(pi0Pt, pi0Yeta) ;
1271 if(TMath::Abs(pi0Y) < 1.0)
1273 fhPrimEtaE ->Fill(pi0E ) ;
1274 fhPrimEtaPt ->Fill(pi0Pt) ;
1275 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
1276 fhPrimEtaYetaYcut ->Fill(pi0Pt,pi0Yeta) ;
1277 fhPrimEtaPtCentrality->Fill(pi0Pt,cen) ;
1278 fhPrimEtaPtEventPlane->Fill(pi0Pt,ep ) ;
1280 fhPrimEtaY ->Fill(pi0Pt, pi0Y) ;
1281 fhPrimEtaYeta->Fill(pi0Pt, pi0Yeta) ;
1285 if(fFillOriginHisto)
1287 Int_t momindex = prim->GetFirstMother();
1290 TParticle* mother = stack->Particle(momindex);
1291 Int_t mompdg = TMath::Abs(mother->GetPdgCode());
1292 Int_t momstatus = mother->GetStatusCode();
1295 if (momstatus == 21)fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1296 else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1297 else if(mompdg > 2100 && mompdg < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1298 else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1299 else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1300 else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1301 else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1302 else if(mompdg >= 310 && mompdg <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1303 else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1304 else if(momstatus == 11 || momstatus == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1305 else fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1309 if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1310 else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1311 else if(mompdg > 2100 && mompdg < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1312 else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1313 else if(momstatus == 11 || momstatus == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1314 else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1315 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1320 //Check if both photons hit Calorimeter
1321 if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
1322 Int_t iphot1=prim->GetFirstDaughter() ;
1323 Int_t iphot2=prim->GetLastDaughter() ;
1324 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack())
1326 TParticle * phot1 = stack->Particle(iphot1) ;
1327 TParticle * phot2 = stack->Particle(iphot2) ;
1328 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22)
1330 //printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
1331 // phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
1333 TLorentzVector lv1, lv2,lvmeson;
1334 phot1->Momentum(lv1);
1335 phot2->Momentum(lv2);
1336 prim ->Momentum(lvmeson);
1338 if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
1340 Bool_t inacceptance = kFALSE;
1341 if(fCalorimeter == "PHOS")
1343 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
1347 if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
1348 inacceptance = kTRUE;
1349 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1353 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1354 inacceptance = kTRUE ;
1355 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1358 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
1360 if(GetEMCALGeometry())
1365 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1366 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1368 if( absID1 >= 0 && absID2 >= 0)
1369 inacceptance = kTRUE;
1371 if(inacceptance && fCheckAccInSector)
1373 Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
1374 Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
1377 Bool_t sameSector = kFALSE;
1378 for(Int_t isector = 0; isector < fNModules/2; isector++)
1381 if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
1384 if(sm1!=sm2 && !sameSector) inacceptance = kFALSE;
1386 //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
1389 // if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1390 // inacceptance = kTRUE;
1391 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1395 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1396 inacceptance = kTRUE ;
1397 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1403 Float_t asym = TMath::Abs((lv1.E()-lv2.E()) / (lv1.E()+lv2.E()));
1404 Double_t angle = lv1.Angle(lv2.Vect());
1408 fhPrimPi0AccE ->Fill(pi0E) ;
1409 fhPrimPi0AccPt ->Fill(pi0Pt) ;
1410 fhPrimPi0AccPhi ->Fill(pi0Pt, phi) ;
1411 fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
1412 fhPrimPi0AccYeta->Fill(pi0Pt,pi0Yeta) ;
1413 fhPrimPi0AccPtCentrality->Fill(pi0Pt,cen) ;
1414 fhPrimPi0AccPtEventPlane->Fill(pi0Pt,ep ) ;
1418 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1419 if(pi0Pt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle);
1420 fhPrimPi0CosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1425 fhPrimEtaAccE ->Fill(pi0E ) ;
1426 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1427 fhPrimEtaAccPhi ->Fill(pi0Pt, phi) ;
1428 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
1429 fhPrimEtaAccYeta->Fill(pi0Pt, pi0Yeta) ;
1430 fhPrimEtaAccPtCentrality->Fill(pi0Pt,cen) ;
1431 fhPrimEtaAccPtEventPlane->Fill(pi0Pt,ep ) ;
1435 fhPrimEtaOpeningAngle ->Fill(pi0Pt,angle);
1436 if(pi0Pt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle);
1437 fhPrimEtaCosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1442 }//Check daughters exist
1443 }// Primary pi0 or eta
1444 }//loop on primaries
1445 }//stack exists and data is MC
1447 else if(GetReader()->ReadAODMCParticles())
1449 TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1452 Int_t nprim = mcparticles->GetEntriesFast();
1454 for(Int_t i=0; i < nprim; i++)
1456 if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
1458 AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
1460 // Only generator particles, when they come from PYTHIA, PHOJET, HERWIG ...
1461 //if( prim->GetStatus() == 0 && (GetMCAnalysisUtils()->GetMCGenerator()).Length()!=0) break;
1463 Int_t pdg = prim->GetPdgCode();
1464 if( pdg == 111 || pdg == 221)
1466 Double_t pi0Pt = prim->Pt() ;
1467 Double_t pi0E = prim->E() ;
1468 //printf("pi0, pt %2.2f, eta %f, phi %f\n",pi0Pt, prim->Eta(), prim->Phi());
1470 if(pi0E == TMath::Abs(prim->Pz()))
1472 if( pdg == 111 ) fhPrimPi0PtRejected->Fill(pi0Pt);
1473 else fhPrimEtaPtRejected->Fill(pi0Pt);
1474 continue ; //Protection against floating point exception
1477 Double_t pi0Y = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
1478 Double_t pi0Yeta = prim->Eta() ;
1479 Double_t phi = TMath::RadToDeg()*prim->Phi() ;
1483 if(TMath::Abs(pi0Y) < 1)
1485 fhPrimPi0E ->Fill(pi0E ) ;
1486 fhPrimPi0Pt ->Fill(pi0Pt) ;
1487 fhPrimPi0Phi->Fill(pi0Pt, phi) ;
1488 fhPrimPi0YetaYcut ->Fill(pi0Pt,pi0Yeta) ;
1489 fhPrimPi0PtCentrality->Fill(pi0Pt,cen) ;
1490 fhPrimPi0PtEventPlane->Fill(pi0Pt,ep ) ;
1492 fhPrimPi0Y ->Fill(pi0Pt, pi0Y ) ;
1493 fhPrimPi0Yeta->Fill(pi0Pt, pi0Yeta) ;
1497 if(TMath::Abs(pi0Y) < 1)
1499 fhPrimEtaE ->Fill(pi0E ) ;
1500 fhPrimEtaPt ->Fill(pi0Pt) ;
1501 fhPrimEtaPhi->Fill(pi0Pt, phi) ;
1502 fhPrimEtaYetaYcut ->Fill(pi0Pt,pi0Yeta) ;
1503 fhPrimEtaPtCentrality->Fill(pi0Pt,cen) ;
1504 fhPrimEtaPtEventPlane->Fill(pi0Pt,ep ) ;
1506 fhPrimEtaY ->Fill(pi0Pt, pi0Y ) ;
1507 fhPrimEtaYeta->Fill(pi0Pt, pi0Yeta) ;
1511 Int_t momindex = prim->GetMother();
1514 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1515 Int_t mompdg = TMath::Abs(mother->GetPdgCode());
1516 Int_t momstatus = mother->GetStatus();
1517 if(fFillOriginHisto)
1521 if (momstatus == 21) fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1522 else if(mompdg < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1523 else if(mompdg > 2100 && mompdg < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1524 else if(mompdg == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1525 else if(mompdg == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1526 else if(mompdg == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1527 else if(mompdg == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1528 else if(mompdg >= 310 && mompdg <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1529 else if(mompdg == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1530 else if(momstatus == 11 || momstatus == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1531 else fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1535 if (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1536 else if(mompdg < 22 ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1537 else if(mompdg > 2100 && mompdg < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1538 else if(mompdg == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1539 else if(momstatus == 11 || momstatus == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1540 else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1541 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1546 //Check if both photons hit Calorimeter
1547 if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
1548 Int_t iphot1=prim->GetDaughter(0) ;
1549 Int_t iphot2=prim->GetDaughter(1) ;
1550 if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim)
1552 AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);
1553 AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
1554 if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22)
1556 TLorentzVector lv1, lv2,lvmeson;
1557 lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1558 lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1559 lvmeson.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
1561 if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
1563 Bool_t inacceptance = kFALSE;
1564 if(fCalorimeter == "PHOS")
1566 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
1570 Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
1571 Double_t vtx2[]={phot2->Xv(),phot2->Yv(),phot2->Zv()};
1572 if(GetPHOSGeometry()->ImpactOnEmc(vtx, phot1->Theta(),phot1->Phi(),mod,z,x) &&
1573 GetPHOSGeometry()->ImpactOnEmc(vtx2,phot2->Theta(),phot2->Phi(),mod,z,x))
1574 inacceptance = kTRUE;
1575 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1579 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1580 inacceptance = kTRUE ;
1581 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1584 else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
1586 if(GetEMCALGeometry())
1591 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1592 GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1594 if( absID1 >= 0 && absID2 >= 0)
1595 inacceptance = kTRUE;
1597 if(inacceptance && fCheckAccInSector)
1599 Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
1600 Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
1603 Bool_t sameSector = kFALSE;
1604 for(Int_t isector = 0; isector < fNModules/2; isector++)
1607 if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
1610 if(sm1!=sm2 && !sameSector) inacceptance = kFALSE;
1612 //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
1615 if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1619 if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1620 inacceptance = kTRUE ;
1621 if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1627 Float_t asym = TMath::Abs((lv1.E()-lv2.E()) / (lv1.E()+lv2.E()));
1628 Double_t angle = lv1.Angle(lv2.Vect());
1632 // printf("ACCEPTED pi0: pt %2.2f, phi %3.2f, eta %1.2f\n",pi0Pt,phi,pi0Y);
1633 fhPrimPi0AccE ->Fill(pi0E ) ;
1634 fhPrimPi0AccPt ->Fill(pi0Pt) ;
1635 fhPrimPi0AccPhi ->Fill(pi0Pt, phi) ;
1636 fhPrimPi0AccY ->Fill(pi0Pt, pi0Y) ;
1637 fhPrimPi0AccYeta->Fill(pi0Pt, pi0Yeta) ;
1638 fhPrimPi0AccPtCentrality->Fill(pi0Pt,cen) ;
1639 fhPrimPi0AccPtEventPlane->Fill(pi0Pt,ep ) ;
1643 fhPrimPi0OpeningAngle ->Fill(pi0Pt,angle);
1644 if(pi0Pt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle);
1645 fhPrimPi0CosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1650 fhPrimEtaAccE ->Fill(pi0E ) ;
1651 fhPrimEtaAccPt ->Fill(pi0Pt) ;
1652 fhPrimEtaAccPhi ->Fill(pi0Pt, phi) ;
1653 fhPrimEtaAccY ->Fill(pi0Pt, pi0Y) ;
1654 fhPrimEtaAccYeta->Fill(pi0Pt, pi0Yeta) ;
1655 fhPrimEtaAccPtCentrality->Fill(pi0Pt,cen) ;
1656 fhPrimEtaAccPtEventPlane->Fill(pi0Pt,ep ) ;
1660 fhPrimEtaOpeningAngle ->Fill(pi0Pt,angle);
1661 if(pi0Pt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle);
1662 fhPrimEtaCosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1667 }//Check daughters exist
1668 }// Primary pi0 or eta
1669 }//loop on primaries
1670 }//stack exists and data is MC
1676 //__________________________________________________________________________________
1677 void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg, TLorentzVector meson,
1678 TLorentzVector daugh1, TLorentzVector daugh2)
1680 // Fill armenteros plots
1682 // Get pTArm and AlphaArm
1683 Float_t momentumSquaredMother = meson.P()*meson.P();
1684 Float_t momentumDaughter1AlongMother = 0.;
1685 Float_t momentumDaughter2AlongMother = 0.;
1687 if (momentumSquaredMother > 0.)
1689 momentumDaughter1AlongMother = (daugh1.Px()*meson.Px() + daugh1.Py()*meson.Py()+ daugh1.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
1690 momentumDaughter2AlongMother = (daugh2.Px()*meson.Px() + daugh2.Py()*meson.Py()+ daugh2.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
1693 Float_t momentumSquaredDaughter1 = daugh1.P()*daugh1.P();
1694 Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
1697 if (ptArmSquared > 0.)
1698 pTArm = sqrt(ptArmSquared);
1700 Float_t alphaArm = 0.;
1701 if(momentumDaughter1AlongMother +momentumDaughter2AlongMother > 0)
1702 alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
1704 TLorentzVector daugh1Boost = daugh1;
1705 daugh1Boost.Boost(-meson.BoostVector());
1706 Float_t cosThStar=TMath::Cos(daugh1Boost.Vect().Angle(meson.Vect()));
1708 Float_t en = meson.Energy();
1710 if(en > 8 && en <= 12) ebin = 0;
1711 if(en > 12 && en <= 16) ebin = 1;
1712 if(en > 16 && en <= 20) ebin = 2;
1713 if(en > 20) ebin = 3;
1714 if(ebin < 0 || ebin > 3) return ;
1719 fhCosThStarPrimPi0->Fill(en,cosThStar);
1720 fhArmPrimPi0[ebin]->Fill(alphaArm,pTArm);
1724 fhCosThStarPrimEta->Fill(en,cosThStar);
1725 fhArmPrimEta[ebin]->Fill(alphaArm,pTArm);
1731 if(daugh1.E()+daugh2.E() > 0) asym = TMath::Abs(daugh1.E()-daugh2.E())/(daugh1.E()+daugh2.E());
1733 printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
1734 en,alphaArm,pTArm,cosThStar,asym);
1738 //_______________________________________________________________________________________
1739 void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1, Int_t index2,
1740 Float_t pt1, Float_t pt2,
1741 Int_t ncell1, Int_t ncell2,
1742 Double_t mass, Double_t pt, Double_t asym,
1743 Double_t deta, Double_t dphi)
1745 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1746 //Adjusted for Pythia, need to see what to do for other generators.
1747 //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1748 // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1750 if(!fFillOriginHisto) return;
1753 Int_t ancStatus = 0;
1754 TLorentzVector ancMomentum;
1755 TVector3 prodVertex;
1756 Int_t ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1757 GetReader(), ancPDG, ancStatus,ancMomentum, prodVertex);
1759 Int_t momindex = -1;
1761 Int_t momstatus = -1;
1762 if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1763 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1766 if(ancPDG==22){//gamma
1767 fhMCOrgMass[0]->Fill(pt,mass);
1768 fhMCOrgAsym[0]->Fill(pt,asym);
1769 fhMCOrgDeltaEta[0]->Fill(pt,deta);
1770 fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1772 else if(TMath::Abs(ancPDG)==11){//e
1773 fhMCOrgMass[1]->Fill(pt,mass);
1774 fhMCOrgAsym[1]->Fill(pt,asym);
1775 fhMCOrgDeltaEta[1]->Fill(pt,deta);
1776 fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1778 else if(ancPDG==111){//Pi0
1779 fhMCOrgMass[2]->Fill(pt,mass);
1780 fhMCOrgAsym[2]->Fill(pt,asym);
1781 fhMCOrgDeltaEta[2]->Fill(pt,deta);
1782 fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
1783 if(fMultiCutAnaSim){
1784 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1785 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1786 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1787 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1788 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1789 asym < fAsymCuts[iasym] &&
1790 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1791 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1792 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1793 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1794 }//pass the different cuts
1795 }// pid bit cut loop
1798 }//Multi cut ana sim
1800 fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1801 if(mass < 0.17 && mass > 0.1) {
1802 fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1803 if(fFillOriginHisto)
1805 if(GetReader()->ReadStack())
1807 TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1808 momindex = ancestor->GetFirstMother();
1809 if(momindex < 0) return;
1810 TParticle* mother = GetMCStack()->Particle(momindex);
1811 mompdg = TMath::Abs(mother->GetPdgCode());
1812 momstatus = mother->GetStatusCode();
1816 TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1817 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1818 momindex = ancestor->GetMother();
1819 if(momindex < 0) return;
1820 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1821 mompdg = TMath::Abs(mother->GetPdgCode());
1822 momstatus = mother->GetStatus();
1825 if (momstatus == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1826 else if(mompdg < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1827 else if(mompdg > 2100 && mompdg < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1828 else if(mompdg == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1829 else if(mompdg == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1830 else if(mompdg == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1831 else if(mompdg == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1832 else if(mompdg >= 310 && mompdg <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1833 else if(mompdg == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1834 else if(momstatus == 11 || momstatus == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
1835 else fhMCPi0PtOrigin->Fill(pt,7.5);//other?
1841 else if(ancPDG==221){//Eta
1842 fhMCOrgMass[3]->Fill(pt,mass);
1843 fhMCOrgAsym[3]->Fill(pt,asym);
1844 fhMCOrgDeltaEta[3]->Fill(pt,deta);
1845 fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1846 if(fMultiCutAnaSim){
1847 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1848 for(Int_t icell=0; icell<fNCellNCuts; icell++){
1849 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1850 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1851 if(pt1 > fPtCuts[ipt] && pt2 > fPtCuts[ipt] &&
1852 asym < fAsymCuts[iasym] &&
1853 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1854 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1855 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1856 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1857 }//pass the different cuts
1858 }// pid bit cut loop
1861 } //Multi cut ana sim
1863 fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1864 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1866 if(fFillOriginHisto)
1868 if(GetReader()->ReadStack())
1870 TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1871 momindex = ancestor->GetFirstMother();
1872 if(momindex < 0) return;
1873 TParticle* mother = GetMCStack()->Particle(momindex);
1874 mompdg = TMath::Abs(mother->GetPdgCode());
1875 momstatus = mother->GetStatusCode();
1879 TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1880 AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1881 momindex = ancestor->GetMother();
1882 if(momindex < 0) return;
1883 AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1884 mompdg = TMath::Abs(mother->GetPdgCode());
1885 momstatus = mother->GetStatus();
1888 if (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1889 else if(mompdg < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1890 else if(mompdg > 2100 && mompdg < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1891 else if(mompdg == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1892 else if(momstatus == 11 || momstatus == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1893 else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1894 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1898 else if(ancPDG==-2212){//AProton
1899 fhMCOrgMass[4]->Fill(pt,mass);
1900 fhMCOrgAsym[4]->Fill(pt,asym);
1901 fhMCOrgDeltaEta[4]->Fill(pt,deta);
1902 fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1904 else if(ancPDG==-2112){//ANeutron
1905 fhMCOrgMass[5]->Fill(pt,mass);
1906 fhMCOrgAsym[5]->Fill(pt,asym);
1907 fhMCOrgDeltaEta[5]->Fill(pt,deta);
1908 fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1910 else if(TMath::Abs(ancPDG)==13){//muons
1911 fhMCOrgMass[6]->Fill(pt,mass);
1912 fhMCOrgAsym[6]->Fill(pt,asym);
1913 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1914 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1916 else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1917 if(ancStatus==1){//Stable particles, converted? not decayed resonances
1918 fhMCOrgMass[6]->Fill(pt,mass);
1919 fhMCOrgAsym[6]->Fill(pt,asym);
1920 fhMCOrgDeltaEta[6]->Fill(pt,deta);
1921 fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1923 else{//resonances and other decays, more hadron conversions?
1924 fhMCOrgMass[7]->Fill(pt,mass);
1925 fhMCOrgAsym[7]->Fill(pt,asym);
1926 fhMCOrgDeltaEta[7]->Fill(pt,deta);
1927 fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
1930 else {//Partons, colliding protons, strings, intermediate corrections
1931 if(ancStatus==11 || ancStatus==12){//String fragmentation
1932 fhMCOrgMass[8]->Fill(pt,mass);
1933 fhMCOrgAsym[8]->Fill(pt,asym);
1934 fhMCOrgDeltaEta[8]->Fill(pt,deta);
1935 fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
1937 else if (ancStatus==21){
1938 if(ancLabel < 2) {//Colliding protons
1939 fhMCOrgMass[11]->Fill(pt,mass);
1940 fhMCOrgAsym[11]->Fill(pt,asym);
1941 fhMCOrgDeltaEta[11]->Fill(pt,deta);
1942 fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
1943 }//colliding protons
1944 else if(ancLabel < 6){//partonic initial states interactions
1945 fhMCOrgMass[9]->Fill(pt,mass);
1946 fhMCOrgAsym[9]->Fill(pt,asym);
1947 fhMCOrgDeltaEta[9]->Fill(pt,deta);
1948 fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
1950 else if(ancLabel < 8){//Final state partons radiations?
1951 fhMCOrgMass[10]->Fill(pt,mass);
1952 fhMCOrgAsym[10]->Fill(pt,asym);
1953 fhMCOrgDeltaEta[10]->Fill(pt,deta);
1954 fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
1957 // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1958 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1962 // printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1963 // ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1965 }////Partons, colliding protons, strings, intermediate corrections
1967 else { //ancLabel <= -1
1968 //printf("Not related at all label = %d\n",ancLabel);
1969 fhMCOrgMass[12]->Fill(pt,mass);
1970 fhMCOrgAsym[12]->Fill(pt,asym);
1971 fhMCOrgDeltaEta[12]->Fill(pt,deta);
1972 fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
1976 //__________________________________________
1977 void AliAnaPi0::MakeAnalysisFillHistograms()
1979 //Process one event and extract photons from AOD branch
1980 // filled with AliAnaPhoton and fill histos with invariant mass
1982 //In case of simulated data, fill acceptance histograms
1983 if(IsDataMC())FillAcceptanceHistograms();
1985 //if (GetReader()->GetEventNumber()%10000 == 0)
1986 // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
1988 if(!GetInputAODBranch())
1990 printf("AliAnaPi0::MakeAnalysisFillHistograms() - No input aod photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1994 //Init some variables
1995 Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
1998 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
2000 //If less than photon 2 entries in the list, skip this event
2004 printf("AliAnaPi0::MakeAnalysisFillHistograms() - nPhotons %d, cent bin %d continue to next event\n",nPhot, GetEventCentrality());
2006 if(GetNCentrBin() > 1) fhCentralityNoPair->Fill(GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt());
2011 Int_t ncentr = GetNCentrBin();
2012 if(GetNCentrBin()==0) ncentr = 1;
2017 Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
2018 Int_t evtIndex1 = 0 ;
2019 Int_t currentEvtIndex = -1;
2020 Int_t curCentrBin = GetEventCentralityBin();
2021 //Int_t curVzBin = GetEventVzBin();
2022 //Int_t curRPBin = GetEventRPBin();
2023 Int_t eventbin = GetEventMixBin();
2025 if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
2027 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());
2031 //Get shower shape information of clusters
2032 TObjArray *clusters = 0;
2033 if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2034 else if(fCalorimeter=="PHOS" ) clusters = GetPHOSClusters() ;
2036 //---------------------------------
2037 //First loop on photons/clusters
2038 //---------------------------------
2039 for(Int_t i1=0; i1<nPhot-1; i1++)
2041 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2042 //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
2044 // get the event index in the mixed buffer where the photon comes from
2045 // in case of mixing with analysis frame, not own mixing
2046 evtIndex1 = GetEventIndex(p1, vert) ;
2047 //printf("charge = %d\n", track->Charge());
2048 if ( evtIndex1 == -1 )
2050 if ( evtIndex1 == -2 )
2053 //printf("z vertex %f < %f\n",vert[2],GetZvertexCut());
2054 if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ; //vertex cut
2057 if (evtIndex1 != currentEvtIndex)
2059 //Fill event bin info
2060 if(DoOwnMix()) fhEventBin->Fill(eventbin) ;
2061 if(GetNCentrBin() > 1)
2063 fhCentrality->Fill(curCentrBin);
2064 if(GetNRPBin() > 1 && GetEventPlane()) fhEventPlaneResolution->Fill(curCentrBin,TMath::Cos(2.*GetEventPlane()->GetQsubRes()));
2066 currentEvtIndex = evtIndex1 ;
2069 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
2071 //Get the momentum of this cluster
2072 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2074 //Get (Super)Module number of this cluster
2075 module1 = GetModuleNumber(p1);
2077 //------------------------------------------
2078 //Get index in VCaloCluster array
2079 AliVCluster *cluster1 = 0;
2080 Bool_t bFound1 = kFALSE;
2081 Int_t caloLabel1 = p1->GetCaloLabel(0);
2085 for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
2086 AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
2089 if (cluster->GetID()==caloLabel1)
2098 }// calorimeter clusters loop
2100 //---------------------------------
2101 //Second loop on photons/clusters
2102 //---------------------------------
2103 for(Int_t i2=i1+1; i2<nPhot; i2++)
2105 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
2107 //In case of mixing frame, check we are not in the same event as the first cluster
2108 Int_t evtIndex2 = GetEventIndex(p2, vert) ;
2109 if ( evtIndex2 == -1 )
2111 if ( evtIndex2 == -2 )
2113 if (GetMixedEvent() && (evtIndex1 == evtIndex2))
2116 //------------------------------------------
2117 //Get index in VCaloCluster array
2118 AliVCluster *cluster2 = 0;
2119 Bool_t bFound2 = kFALSE;
2120 Int_t caloLabel2 = p2->GetCaloLabel(0);
2122 for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
2123 AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
2125 if(cluster->GetID()==caloLabel2) {
2131 }// calorimeter clusters loop
2136 if(cluster1 && bFound1){
2137 tof1 = cluster1->GetTOF()*1e9;
2138 l01 = cluster1->GetM02();
2140 // else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
2141 // p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
2145 if(cluster2 && bFound2)
2147 tof2 = cluster2->GetTOF()*1e9;
2148 l02 = cluster2->GetM02();
2151 // else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
2152 // p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
2156 Double_t t12diff = tof1-tof2;
2157 if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2159 //------------------------------------------
2161 //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
2163 //Get the momentum of this cluster
2164 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2166 module2 = GetModuleNumber(p2);
2168 //---------------------------------
2169 // Get pair kinematics
2170 //---------------------------------
2171 Double_t m = (photon1 + photon2).M() ;
2172 Double_t pt = (photon1 + photon2).Pt();
2173 Double_t deta = photon1.Eta() - photon2.Eta();
2174 Double_t dphi = photon1.Phi() - photon2.Phi();
2175 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2178 printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
2180 //--------------------------------
2181 // Opening angle selection
2182 //--------------------------------
2183 //Check if opening angle is too large or too small compared to what is expected
2184 Double_t angle = photon1.Angle(photon2.Vect());
2185 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
2187 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2191 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2193 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
2197 //-------------------------------------------------------------------------------------------------
2198 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2199 //-------------------------------------------------------------------------------------------------
2200 if(a < fAsymCuts[0] && fFillSMCombinations)
2202 if(module1==module2 && module1 >=0 && module1<fNModules)
2203 fhReMod[module1]->Fill(pt,m) ;
2205 if(fCalorimeter=="EMCAL")
2209 for(Int_t i = 0; i < fNModules/2; i++)
2212 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
2216 for(Int_t i = 0; i < fNModules-2; i++){
2217 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
2221 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
2222 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
2223 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
2227 //In case we want only pairs in same (super) module, check their origin.
2229 if(fSameSM && module1!=module2) ok=kFALSE;
2232 //Check if one of the clusters comes from a conversion
2233 if(fCheckConversion)
2235 if (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
2236 else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
2239 // Fill shower shape cut histograms
2240 if(fFillSSCombinations)
2242 if ( l01 > 0.01 && l01 < 0.4 &&
2243 l02 > 0.01 && l02 < 0.4 ) fhReSS[0]->Fill(pt,m); // Tight
2244 else if( l01 > 0.4 && l02 > 0.4 ) fhReSS[1]->Fill(pt,m); // Loose
2245 else if( l01 > 0.01 && l01 < 0.4 && l02 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2246 else if( l02 > 0.01 && l02 < 0.4 && l01 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2249 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2250 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2251 if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
2252 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
2253 if(a < fAsymCuts[iasym])
2255 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2256 //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);
2258 if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2260 fhRe1 [index]->Fill(pt,m);
2261 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
2262 if(fFillBadDistHisto){
2263 if(p1->DistToBad()>0 && p2->DistToBad()>0){
2264 fhRe2 [index]->Fill(pt,m) ;
2265 if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
2266 if(p1->DistToBad()>1 && p2->DistToBad()>1){
2267 fhRe3 [index]->Fill(pt,m) ;
2268 if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
2271 }// Fill bad dist histos
2273 }// asymmetry cut loop
2277 //Fill histograms with opening angle
2280 fhRealOpeningAngle ->Fill(pt,angle);
2281 fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2284 //Fill histograms with pair assymmetry
2285 if(fFillAsymmetryHisto)
2287 fhRePtAsym->Fill(pt,a);
2288 if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
2289 if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
2292 //-------------------------------------------------------
2293 //Get the number of cells needed for multi cut analysis.
2294 //-------------------------------------------------------
2297 if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim))
2299 AliVEvent * event = GetReader()->GetInputEvent();
2301 for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++)
2303 AliVCluster *cluster = event->GetCaloCluster(iclus);
2306 if (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
2307 else if(fCalorimeter == "PHOS" && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
2310 if (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
2311 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
2312 } // PHOS or EMCAL cluster as requested in analysis
2314 if(ncell2 > 0 && ncell1 > 0) break; // No need to continue the iteration
2317 //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
2324 //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
2327 if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2328 GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2330 fhReMCFromConversion->Fill(pt,m);
2332 else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2333 !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2335 fhReMCFromNotConversion->Fill(pt,m);
2339 fhReMCFromMixConversion->Fill(pt,m);
2342 FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
2345 //-----------------------
2346 //Multi cuts analysis
2347 //-----------------------
2350 //Histograms for different PID bits selection
2351 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2353 if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton) &&
2354 p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) fhRePIDBits[ipid]->Fill(pt,m) ;
2356 //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
2357 } // pid bit cut loop
2359 //Several pt,ncell and asymmetry cuts
2360 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2361 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2362 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2363 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2364 if(p1->E() > fPtCuts[ipt] && p2->E() > fPtCuts[ipt] &&
2365 a < fAsymCuts[iasym] &&
2366 ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
2367 fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
2368 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2369 if(fFillSMCombinations && module1==module2){
2370 fhRePtNCellAsymCutsSM[module1][index]->Fill(pt,m) ;
2373 }// pid bit cut loop
2376 if(GetHistogramRanges()->GetHistoTrackMultiplicityBins()){
2377 for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
2378 if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
2381 }// multiple cuts analysis
2383 }// second same event particle
2386 //-------------------------------------------------------------
2388 //-------------------------------------------------------------
2391 //Recover events in with same characteristics as the current event
2393 //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2394 if(eventbin < 0) return ;
2396 TList * evMixList=fEventsList[eventbin] ;
2400 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mix event list not available, bin %d \n",eventbin);
2404 Int_t nMixed = evMixList->GetSize() ;
2405 for(Int_t ii=0; ii<nMixed; ii++)
2407 TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2408 Int_t nPhot2=ev2->GetEntriesFast() ;
2411 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, GetEventCentralityBin());
2413 fhEventMixBin->Fill(eventbin) ;
2415 //---------------------------------
2416 //First loop on photons/clusters
2417 //---------------------------------
2418 for(Int_t i1=0; i1<nPhot; i1++){
2419 AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2420 if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2422 //Get kinematics of cluster and (super) module of this cluster
2423 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2424 module1 = GetModuleNumber(p1);
2426 //---------------------------------
2427 //First loop on photons/clusters
2428 //---------------------------------
2429 for(Int_t i2=0; i2<nPhot2; i2++){
2430 AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
2432 //Get kinematics of second cluster and calculate those of the pair
2433 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2434 m = (photon1+photon2).M() ;
2435 Double_t pt = (photon1 + photon2).Pt();
2436 Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2438 //Check if opening angle is too large or too small compared to what is expected
2439 Double_t angle = photon1.Angle(photon2.Vect());
2440 if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
2442 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2445 if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2447 printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
2453 printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
2454 p1->Pt(), p2->Pt(), pt,m,a);
2456 //In case we want only pairs in same (super) module, check their origin.
2457 module2 = GetModuleNumber(p2);
2459 //-------------------------------------------------------------------------------------------------
2460 //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2461 //-------------------------------------------------------------------------------------------------
2462 if(a < fAsymCuts[0] && fFillSMCombinations){
2463 if(module1==module2 && module1 >=0 && module1<fNModules)
2464 fhMiMod[module1]->Fill(pt,m) ;
2466 if(fCalorimeter=="EMCAL"){
2470 for(Int_t i = 0; i < fNModules/2; i++){
2472 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
2476 for(Int_t i = 0; i < fNModules-2; i++){
2477 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
2481 if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2482 if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2483 if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
2490 if(fSameSM && module1!=module2) ok=kFALSE;
2493 //Check if one of the clusters comes from a conversion
2494 if(fCheckConversion){
2495 if (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2496 else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2498 //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2499 for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2500 if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
2502 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2504 if(a < fAsymCuts[iasym])
2506 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2508 if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2510 fhMi1 [index]->Fill(pt,m) ;
2511 if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
2512 if(fFillBadDistHisto)
2514 if(p1->DistToBad()>0 && p2->DistToBad()>0)
2516 fhMi2 [index]->Fill(pt,m) ;
2517 if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
2518 if(p1->DistToBad()>1 && p2->DistToBad()>1)
2520 fhMi3 [index]->Fill(pt,m) ;
2521 if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
2524 }// Fill bad dist histo
2528 }//loop for histograms
2530 //-----------------------
2531 //Multi cuts analysis
2532 //-----------------------
2534 //Several pt,ncell and asymmetry cuts
2536 for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2537 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2538 for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2539 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2540 if(p1->Pt() > fPtCuts[ipt] && p2->Pt() > fPtCuts[ipt] &&
2541 a < fAsymCuts[iasym] //&&
2542 //p1->GetBtag() >= fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2544 fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
2545 //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym, fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2547 }// pid bit cut loop
2552 //Fill histograms with opening angle
2553 if(fFillAngleHisto){
2554 fhMixedOpeningAngle ->Fill(pt,angle);
2555 fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2559 }// second cluster loop
2560 }//first cluster loop
2561 }//loop on mixed events
2563 //--------------------------------------------------------
2564 //Add the current event to the list of events for mixing
2565 //--------------------------------------------------------
2566 TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2567 //Add current event to buffer and Remove redundant events
2568 if(currentEvent->GetEntriesFast()>0){
2569 evMixList->AddFirst(currentEvent) ;
2570 currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2571 if(evMixList->GetSize() >= GetNMaxEvMix())
2573 TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2574 evMixList->RemoveLast() ;
2579 delete currentEvent ;
2586 //________________________________________________________________________
2587 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2589 // retieves the event index and checks the vertex
2590 // in the mixed buffer returns -2 if vertex NOK
2591 // for normal events returns 0 if vertex OK and -1 if vertex NOK
2593 Int_t evtIndex = -1 ;
2594 if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2596 if (GetMixedEvent()){
2598 evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2599 GetVertex(vert,evtIndex);
2601 if(TMath::Abs(vert[2])> GetZvertexCut())
2602 evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2603 } else {// Single event
2607 if(TMath::Abs(vert[2])> GetZvertexCut())
2608 evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)