]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/CaloTrackCorrelations/AliAnaInsideClusterInvariantMass.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Restored Makefile needed for PAR creation. Apparently this disappeared when
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaInsideClusterInvariantMass.cxx
CommitLineData
992b14a7 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//_________________________________________________________________________
17//
18// Split clusters with some criteria and calculate invariant mass
19// to identify them as pi0 or conversion
20//
21//
22//-- Author: Gustavo Conesa (LPSC-Grenoble)
23//_________________________________________________________________________
24
25//////////////////////////////////////////////////////////////////////////////
26
27
28// --- ROOT system ---
29#include <TList.h>
30#include <TClonesArray.h>
31#include <TObjString.h>
32#include <TH3F.h>
29ca9cad 33#include <TCanvas.h>
34#include <TStyle.h>
35#include <TPad.h>
992b14a7 36
37// --- Analysis system ---
38#include "AliAnaInsideClusterInvariantMass.h"
39#include "AliCaloTrackReader.h"
40#include "AliMCAnalysisUtils.h"
41#include "AliStack.h"
42#include "AliFiducialCut.h"
43#include "TParticle.h"
44#include "AliVCluster.h"
45#include "AliAODEvent.h"
46#include "AliAODMCParticle.h"
47#include "AliEMCALGeoParams.h"
48
c5693f62 49// --- Detectors ---
50//#include "AliPHOSGeoUtils.h"
51#include "AliEMCALGeometry.h"
52
992b14a7 53ClassImp(AliAnaInsideClusterInvariantMass)
54
55//__________________________________________________________________
56AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
745913ae 57 AliAnaCaloTrackCorrBaseClass(),
29ca9cad 58 fCalorimeter(""),
71e3889f 59 fM02MaxCut(0), fM02MinCut(0),
60 fMinNCells(0),
243c2909 61 fMassEtaMin(0), fMassEtaMax(0),
62 fMassPi0Min(0), fMassPi0Max(0),
71e3889f 63 fMassConMin(0), fMassConMax(0),
64 fPlotCluster(0)
992b14a7 65{
66 //default ctor
67
68 // Init array of histograms
5c46c992 69 for(Int_t i = 0; i < 7; i++)
70 {
71 for(Int_t j = 0; j < 2; j++)
72 {
73
74 fhMassNLocMax1[i][j] = 0;
75 fhMassNLocMax2[i][j] = 0;
76 fhMassNLocMaxN[i][j] = 0;
77 fhNLocMax[i][j] = 0;
71e3889f 78 fhNLocMaxNLabel[i][j] = 0;
5c46c992 79 fhNLocMaxEMax[i][j] = 0;
80 fhNLocMaxEFrac[i][j] = 0;
81 fhNLocMaxM02Cut[i][j] = 0;
82 fhM02NLocMax1[i][j] = 0;
83 fhM02NLocMax2[i][j] = 0;
84 fhM02NLocMaxN[i][j] = 0;
85 fhNCellNLocMax1[i][j] = 0;
86 fhNCellNLocMax2[i][j] = 0;
87 fhNCellNLocMaxN[i][j] = 0;
88 fhM02Pi0LocMax1[i][j] = 0;
89 fhM02EtaLocMax1[i][j] = 0;
90 fhM02ConLocMax1[i][j] = 0;
91 fhM02Pi0LocMax2[i][j] = 0;
92 fhM02EtaLocMax2[i][j] = 0;
93 fhM02ConLocMax2[i][j] = 0;
94 fhM02Pi0LocMaxN[i][j] = 0;
95 fhM02EtaLocMaxN[i][j] = 0;
96 fhM02ConLocMaxN[i][j] = 0;
97
98 }
99
100 fhTrackMatchedDEtaLocMax1[i] = 0;
101 fhTrackMatchedDPhiLocMax1[i] = 0;
102 fhTrackMatchedDEtaLocMax2[i] = 0;
103 fhTrackMatchedDPhiLocMax2[i] = 0;
104 fhTrackMatchedDEtaLocMaxN[i] = 0;
105 fhTrackMatchedDPhiLocMaxN[i] = 0;
243c2909 106
992b14a7 107 }
108
5c46c992 109 for(Int_t i = 0; i < 2; i++)
110 {
111 fhAnglePairLocMax1 [i] = 0;
112 fhAnglePairLocMax2 [i] = 0;
113 fhAnglePairLocMaxN [i] = 0;
114 fhAnglePairMassLocMax1[i] = 0;
115 fhAnglePairMassLocMax2[i] = 0;
116 fhAnglePairMassLocMaxN[i] = 0;
117 }
118
992b14a7 119 InitParameters();
120
121}
122
123//_______________________________________________________________
124TObjString * AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
125{
126 //Save parameters used for analysis
127 TString parList ; //this will be list of parameters used for this analysis.
128 const Int_t buffersize = 255;
129 char onePar[buffersize] ;
130
131 snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
132 parList+=onePar ;
133
243c2909 134 snprintf(onePar,buffersize,"Calorimeter: %s\n", fCalorimeter.Data()) ;
992b14a7 135 parList+=onePar ;
dbba06ca 136 snprintf(onePar,buffersize,"fLocMaxCutE =%2.2f \n", GetCaloUtils()->GetLocalMaximaCutE()) ;
29ca9cad 137 parList+=onePar ;
dbba06ca 138 snprintf(onePar,buffersize,"fLocMaxCutEDiff =%2.2f \n",GetCaloUtils()->GetLocalMaximaCutEDiff()) ;
29ca9cad 139 parList+=onePar ;
71e3889f 140 snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n", fM02MinCut, fM02MaxCut) ;
992b14a7 141 parList+=onePar ;
e23a0471 142 snprintf(onePar,buffersize,"fMinNCells =%d \n", fMinNCells) ;
992b14a7 143 parList+=onePar ;
243c2909 144 snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f\n", fMassPi0Min,fMassPi0Max);
145 parList+=onePar ;
146 snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f\n", fMassEtaMin,fMassEtaMax);
147 parList+=onePar ;
148 snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f\n", fMassConMin,fMassConMax);
149 parList+=onePar ;
150
992b14a7 151 return new TObjString(parList) ;
152
153}
154
992b14a7 155
156//_____________________________________________________________________________________
157TLorentzVector AliAnaInsideClusterInvariantMass::GetCellMomentum(const Int_t absId,
243c2909 158 Float_t en,
992b14a7 159 AliVCaloCells * cells)
160{
161
162 // Cell momentum calculation for invariant mass
163
164 Double_t cellpos[] = {0, 0, 0};
165 GetEMCALGeometry()->GetGlobal(absId, cellpos);
166
167 if(GetVertex(0)){//calculate direction from vertex
168 cellpos[0]-=GetVertex(0)[0];
169 cellpos[1]-=GetVertex(0)[1];
170 cellpos[2]-=GetVertex(0)[2];
171 }
172
243c2909 173 Double_t r = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]+cellpos[2]*cellpos[2] ) ;
174
175 //If not calculated before, get the energy
29ca9cad 176 if(en <=0 )
243c2909 177 {
178 en = cells->GetCellAmplitude(absId);
dbba06ca 179 GetCaloUtils()->RecalibrateCellAmplitude(en,fCalorimeter,absId);
243c2909 180 }
992b14a7 181
992b14a7 182 TLorentzVector cellMom ;
183 cellMom.SetPxPyPzE( en*cellpos[0]/r, en*cellpos[1]/r, en*cellpos[2]/r, en) ;
184
185 return cellMom;
186
187}
188
189//________________________________________________________________
190TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
191{
192 // Create histograms to be saved in output file and
193 // store them in outputContainer
194 TList * outputContainer = new TList() ;
195 outputContainer->SetName("InsideClusterHistos") ;
196
745913ae 197 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
198 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
199 Int_t mbins = GetHistogramRanges()->GetHistoMassBins(); Float_t mmax = GetHistogramRanges()->GetHistoMassMax(); Float_t mmin = GetHistogramRanges()->GetHistoMassMin();
200 Int_t ncbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t ncmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t ncmin = GetHistogramRanges()->GetHistoNClusterCellMin();
992b14a7 201
5c46c992 202 Int_t nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
203 Float_t resetamax = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
204 Float_t resetamin = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
205 Int_t nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
206 Float_t resphimax = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
207 Float_t resphimin = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
208
992b14a7 209 TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
210 TString pname[] ={"","Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
211
212 Int_t n = 1;
213
214 if(IsDataMC()) n = 7;
215
243c2909 216 Int_t nMaxBins = 10;
217
5c46c992 218 TString sMatched[] = {"","Matched"};
219
220 for(Int_t i = 0; i < n; i++)
221 {
222
223 for(Int_t j = 0; j < 2; j++)
224 {
225
226 fhMassNLocMax1[i][j] = new TH2F(Form("hMassNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
227 Form("Invariant mass of 2 highest energy cells %s %s",ptype[i].Data(),sMatched[j].Data()),
228 nptbins,ptmin,ptmax,mbins,mmin,mmax);
229 fhMassNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
230 fhMassNLocMax1[i][j]->SetXTitle("E (GeV)");
231 outputContainer->Add(fhMassNLocMax1[i][j]) ;
232
233 fhMassNLocMax2[i][j] = new TH2F(Form("hMassNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
234 Form("Invariant mass of 2 local maxima cells %s %s",ptype[i].Data(),sMatched[j].Data()),
235 nptbins,ptmin,ptmax,mbins,mmin,mmax);
236 fhMassNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
237 fhMassNLocMax2[i][j]->SetXTitle("E (GeV)");
238 outputContainer->Add(fhMassNLocMax2[i][j]) ;
239
240 fhMassNLocMaxN[i][j] = new TH2F(Form("hMassNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
241 Form("Invariant mass of N>2 local maxima cells %s %s",ptype[i].Data(),sMatched[j].Data()),
242 nptbins,ptmin,ptmax,mbins,mmin,mmax);
243 fhMassNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
244 fhMassNLocMaxN[i][j]->SetXTitle("E (GeV)");
245 outputContainer->Add(fhMassNLocMaxN[i][j]) ;
246
247
248 fhNLocMax[i][j] = new TH2F(Form("hNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
249 Form("Number of local maxima in cluster %s %s",ptype[i].Data(),sMatched[j].Data()),
250 nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
251 fhNLocMax[i][j] ->SetYTitle("N maxima");
252 fhNLocMax[i][j] ->SetXTitle("E (GeV)");
253 outputContainer->Add(fhNLocMax[i][j]) ;
254
71e3889f 255 if(IsDataMC())
256 {
257 fhNLocMaxNLabel[i][j] = new TH2F(Form("hNLocMaxNLabel%s%s",pname[i].Data(),sMatched[j].Data()),
258 Form("Number of local maxima in cluster vs number of MC labels %s %s",ptype[i].Data(),sMatched[j].Data()),
259 nMaxBins,0,nMaxBins,nMaxBins,0,nMaxBins);
260 fhNLocMaxNLabel[i][j] ->SetYTitle("N maxima");
261 fhNLocMaxNLabel[i][j] ->SetXTitle("N MC labels");
262 outputContainer->Add(fhNLocMaxNLabel[i][j]) ;
263 }
264
5c46c992 265 fhNLocMaxEMax[i][j] = new TH2F(Form("hNLocMaxEMax%s%s",pname[i].Data(),sMatched[j].Data()),
266 Form("Number of local maxima in cluster vs energy of maxima %s %s",ptype[i].Data(),sMatched[j].Data()),
267 nptbins*10,ptmin,ptmax,nMaxBins,0,nMaxBins);
268 fhNLocMaxEMax[i][j] ->SetYTitle("N maxima");
269 fhNLocMaxEMax[i][j] ->SetXTitle("E of maxima (GeV)");
270 outputContainer->Add(fhNLocMaxEMax[i][j]) ;
271
272 fhNLocMaxEFrac[i][j] = new TH2F(Form("hNLocMaxEFrac%s%s",pname[i].Data(),sMatched[j].Data()),
273 Form("Number of local maxima in cluster vs fraction of cluster energy of maxima %s %s",ptype[i].Data(),sMatched[j].Data()),
274 100,0,1,nMaxBins,0,nMaxBins);
275 fhNLocMaxEFrac[i][j] ->SetYTitle("N maxima");
276 fhNLocMaxEFrac[i][j] ->SetXTitle("E maxima / E cluster");
277 outputContainer->Add(fhNLocMaxEFrac[i][j]) ;
278
279 fhNLocMaxM02Cut[i][j] = new TH2F(Form("hNLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
71e3889f 280 Form("Number of local maxima in cluster %s for %2.2f < M02 < %2.2f",ptype[i].Data(),fM02MinCut,fM02MaxCut),
5c46c992 281 nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
282 fhNLocMaxM02Cut[i][j]->SetYTitle("N maxima");
283 fhNLocMaxM02Cut[i][j]->SetXTitle("E (GeV)");
284 outputContainer->Add(fhNLocMaxM02Cut[i][j]) ;
285
286
287 fhM02NLocMax1[i][j] = new TH2F(Form("hM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
288 Form("#lambda_{0}^{2} vs E for N max = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
289 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
290 fhM02NLocMax1[i][j] ->SetYTitle("#lambda_{0}^{2}");
291 fhM02NLocMax1[i][j] ->SetXTitle("E (GeV)");
292 outputContainer->Add(fhM02NLocMax1[i][j]) ;
293
294 fhM02NLocMax2[i][j] = new TH2F(Form("hM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
295 Form("#lambda_{0}^{2} vs E for N max = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
296 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
297 fhM02NLocMax2[i][j] ->SetYTitle("#lambda_{0}^{2}");
298 fhM02NLocMax2[i][j] ->SetXTitle("E (GeV)");
299 outputContainer->Add(fhM02NLocMax2[i][j]) ;
300
301
302 fhM02NLocMaxN[i][j] = new TH2F(Form("hM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
303 Form("#lambda_{0}^{2} vs E for N max > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
304 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
305 fhM02NLocMaxN[i][j] ->SetYTitle("#lambda_{0}^{2}");
306 fhM02NLocMaxN[i][j] ->SetXTitle("E (GeV)");
307 outputContainer->Add(fhM02NLocMaxN[i][j]) ;
308
309
310 fhNCellNLocMax1[i][j] = new TH2F(Form("hNCellNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
311 Form("#lambda_{0}^{2} vs E for N max = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
312 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
313 fhNCellNLocMax1[i][j] ->SetYTitle("N cells");
314 fhNCellNLocMax1[i][j] ->SetXTitle("E (GeV)");
315 outputContainer->Add(fhNCellNLocMax1[i][j]) ;
316
317 fhNCellNLocMax2[i][j] = new TH2F(Form("hNCellNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
318 Form("#lambda_{0}^{2} vs E for N max = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
319 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
320 fhNCellNLocMax2[i][j] ->SetYTitle("N cells");
321 fhNCellNLocMax2[i][j] ->SetXTitle("E (GeV)");
322 outputContainer->Add(fhNCellNLocMax2[i][j]) ;
323
324
325 fhNCellNLocMaxN[i][j] = new TH2F(Form("hNCellNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
326 Form("#lambda_{0}^{2} vs E for N max > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
327 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
328 fhNCellNLocMaxN[i][j] ->SetYTitle("N cells");
329 fhNCellNLocMaxN[i][j] ->SetXTitle("E (GeV)");
330 outputContainer->Add(fhNCellNLocMaxN[i][j]) ;
331
332
333 fhM02Pi0LocMax1[i][j] = new TH2F(Form("hM02Pi0LocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
334 Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 1",fMassPi0Min,fMassPi0Max,ptype[i].Data()),
335 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
336 fhM02Pi0LocMax1[i][j] ->SetYTitle("#lambda_{0}^{2}");
337 fhM02Pi0LocMax1[i][j] ->SetXTitle("E (GeV)");
338 outputContainer->Add(fhM02Pi0LocMax1[i][j]) ;
339
340 fhM02EtaLocMax1[i][j] = new TH2F(Form("hM02EtaLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
341 Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",fMassEtaMin,fMassEtaMax,ptype[i].Data()),
342 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
343 fhM02EtaLocMax1[i][j] ->SetYTitle("#lambda_{0}^{2}");
344 fhM02EtaLocMax1[i][j] ->SetXTitle("E (GeV)");
345 outputContainer->Add(fhM02EtaLocMax1[i][j]) ;
346
347 fhM02ConLocMax1[i][j] = new TH2F(Form("hM02ConLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
348 Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",fMassConMin,fMassConMax,ptype[i].Data()),
349 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
350 fhM02ConLocMax1[i][j] ->SetYTitle("#lambda_{0}^{2}");
351 fhM02ConLocMax1[i][j] ->SetXTitle("E (GeV)");
352 outputContainer->Add(fhM02ConLocMax1[i][j]) ;
353
354 fhM02Pi0LocMax2[i][j] = new TH2F(Form("hM02Pi0LocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
355 Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 2",fMassPi0Min,fMassPi0Max,ptype[i].Data()),
356 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
357 fhM02Pi0LocMax2[i][j] ->SetYTitle("#lambda_{0}^{2}");
358 fhM02Pi0LocMax2[i][j] ->SetXTitle("E (GeV)");
359 outputContainer->Add(fhM02Pi0LocMax2[i][j]) ;
360
361 fhM02EtaLocMax2[i][j] = new TH2F(Form("hM02EtaLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
362 Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",fMassEtaMin,fMassEtaMax,ptype[i].Data()),
363 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
364 fhM02EtaLocMax2[i][j] ->SetYTitle("#lambda_{0}^{2}");
365 fhM02EtaLocMax2[i][j] ->SetXTitle("E (GeV)");
366 outputContainer->Add(fhM02EtaLocMax2[i][j]) ;
367
368 fhM02ConLocMax2[i][j] = new TH2F(Form("hM02ConLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
369 Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",fMassConMin,fMassConMax,ptype[i].Data()),
370 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
371 fhM02ConLocMax2[i][j] ->SetYTitle("#lambda_{0}^{2}");
372 fhM02ConLocMax2[i][j] ->SetXTitle("E (GeV)");
373 outputContainer->Add(fhM02ConLocMax2[i][j]) ;
374
375 fhM02Pi0LocMaxN[i][j] = new TH2F(Form("hM02Pi0LocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
376 Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max > 2",fMassPi0Min,fMassPi0Max,ptype[i].Data()),
377 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
378 fhM02Pi0LocMaxN[i][j] ->SetYTitle("#lambda_{0}^{2}");
379 fhM02Pi0LocMaxN[i][j] ->SetXTitle("E (GeV)");
380 outputContainer->Add(fhM02Pi0LocMaxN[i][j]) ;
381
382 fhM02EtaLocMaxN[i][j] = new TH2F(Form("hM02EtaLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
383 Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max > 2", fMassEtaMin,fMassEtaMax,ptype[i].Data()),
384 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
385 fhM02EtaLocMaxN[i][j] ->SetYTitle("#lambda_{0}^{2}");
386 fhM02EtaLocMaxN[i][j] ->SetXTitle("E (GeV)");
387 outputContainer->Add(fhM02EtaLocMaxN[i][j]) ;
388
389 fhM02ConLocMaxN[i][j] = new TH2F(Form("hM02ConLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
390 Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",fMassConMin,fMassConMax,ptype[i].Data()),
391 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
392 fhM02ConLocMaxN[i][j] ->SetYTitle("#lambda_{0}^{2}");
393 fhM02ConLocMaxN[i][j] ->SetXTitle("E (GeV)");
394 outputContainer->Add(fhM02ConLocMaxN[i][j]) ;
395
396 } // matched, not matched
397
398 } // MC particle list
399
400 for(Int_t i = 0; i < n; i++)
401 {
402 fhTrackMatchedDEtaLocMax1[i] = new TH2F
403 (Form("hTrackMatchedDEtaLocMax1%s",pname[i].Data()),
404 Form("d#eta of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
405 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
406 fhTrackMatchedDEtaLocMax1[i]->SetYTitle("d#eta");
407 fhTrackMatchedDEtaLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
408
409 fhTrackMatchedDPhiLocMax1[i] = new TH2F
410 (Form("hTrackMatchedDPhiLocMax1%s",pname[i].Data()),
411 Form("d#phi of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
412 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
413 fhTrackMatchedDPhiLocMax1[i]->SetYTitle("d#phi (rad)");
414 fhTrackMatchedDPhiLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
415
416 outputContainer->Add(fhTrackMatchedDEtaLocMax1[i]) ;
417 outputContainer->Add(fhTrackMatchedDPhiLocMax1[i]) ;
418
419 fhTrackMatchedDEtaLocMax2[i] = new TH2F
420 (Form("hTrackMatchedDEtaLocMax2%s",pname[i].Data()),
421 Form("d#eta of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
422 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
423 fhTrackMatchedDEtaLocMax2[i]->SetYTitle("d#eta");
424 fhTrackMatchedDEtaLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
425
426 fhTrackMatchedDPhiLocMax2[i] = new TH2F
427 (Form("hTrackMatchedDPhiLocMax2%s",pname[i].Data()),
428 Form("d#phi of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
429 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
430 fhTrackMatchedDPhiLocMax2[i]->SetYTitle("d#phi (rad)");
431 fhTrackMatchedDPhiLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
432
433 outputContainer->Add(fhTrackMatchedDEtaLocMax2[i]) ;
434 outputContainer->Add(fhTrackMatchedDPhiLocMax2[i]) ;
435
436 fhTrackMatchedDEtaLocMaxN[i] = new TH2F
437 (Form("hTrackMatchedDEtaLocMaxN%s",pname[i].Data()),
438 Form("d#eta of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
439 nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
440 fhTrackMatchedDEtaLocMaxN[i]->SetYTitle("d#eta");
441 fhTrackMatchedDEtaLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
442
443 fhTrackMatchedDPhiLocMaxN[i] = new TH2F
444 (Form("hTrackMatchedDPhiLocMaxN%s",pname[i].Data()),
445 Form("d#phi of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
446 nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
447 fhTrackMatchedDPhiLocMaxN[i]->SetYTitle("d#phi (rad)");
448 fhTrackMatchedDPhiLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
449
450 outputContainer->Add(fhTrackMatchedDEtaLocMaxN[i]) ;
451 outputContainer->Add(fhTrackMatchedDPhiLocMaxN[i]) ;
452
453 }
454
455 for(Int_t j = 0; j < 2; j++)
456 {
457
458 fhAnglePairLocMax1[j] = new TH2F(Form("hAnglePairLocMax1%s",sMatched[j].Data()),
459 Form("Opening angle of 2 highest energy cells vs pair Energy, %s",sMatched[j].Data()),
460 nptbins,ptmin,ptmax,200,0,0.2);
461 fhAnglePairLocMax1[j]->SetYTitle("#alpha (rad)");
462 fhAnglePairLocMax1[j]->SetXTitle("E (GeV)");
463 outputContainer->Add(fhAnglePairLocMax1[j]) ;
464
465 fhAnglePairLocMax2[j] = new TH2F(Form("hAnglePairLocMax2%s",sMatched[j].Data()),
466 Form("Opening angle of 2 local maxima cells vs Energy, %s",sMatched[j].Data()),
467 nptbins,ptmin,ptmax,200,0,0.2);
468 fhAnglePairLocMax2[j]->SetYTitle("#alpha (rad)");
469 fhAnglePairLocMax2[j]->SetXTitle("E (GeV)");
470 outputContainer->Add(fhAnglePairLocMax2[j]) ;
471
472 fhAnglePairLocMaxN[j] = new TH2F(Form("hAnglePairLocMaxN%s",sMatched[j].Data()),
473 Form("Opening angle of N>2 local maxima cells vs Energy, %s",sMatched[j].Data()),
474 nptbins,ptmin,ptmax,200,0,0.2);
475 fhAnglePairLocMaxN[j]->SetYTitle("#alpha (rad)");
476 fhAnglePairLocMaxN[j]->SetXTitle("E (GeV)");
477 outputContainer->Add(fhAnglePairLocMaxN[j]) ;
478
479 fhAnglePairMassLocMax1[j] = new TH2F(Form("hAnglePairMassLocMax1%s",sMatched[j].Data()),
480 Form("Opening angle of 2 highest energy cells vs Mass for E > 5 GeV, %s",sMatched[j].Data()),
481 mbins,mmin,mmax,200,0,0.2);
482 fhAnglePairMassLocMax1[j]->SetXTitle("M (GeV/c^{2})");
483 fhAnglePairMassLocMax1[j]->SetYTitle("#alpha (rad)");
484 outputContainer->Add(fhAnglePairMassLocMax1[j]) ;
485
486 fhAnglePairMassLocMax2[j] = new TH2F(Form("hAnglePairMassLocMax2%s",sMatched[j].Data()),
487 Form("Opening angle of 2 local maxima cells vs Mass for E > 5 GeV, %s",sMatched[j].Data()),
488 mbins,mmin,mmax,200,0,0.2);
489 fhAnglePairMassLocMax2[j]->SetXTitle("M (GeV/c^{2})");
490 fhAnglePairMassLocMax2[j]->SetYTitle("#alpha (rad)");
491 outputContainer->Add(fhAnglePairMassLocMax2[j]) ;
492
493 fhAnglePairMassLocMaxN[j] = new TH2F(Form("hAnglePairMassLocMaxN%s",sMatched[j].Data()),
494 Form("Opening angle of N>2 local maxima cells vs Mass for E > 5 GeV, %s",sMatched[j].Data()),
495 mbins,mmin,mmax,200,0,0.2);
496 fhAnglePairMassLocMaxN[j]->SetXTitle("M (GeV/c^{2})");
497 fhAnglePairMassLocMaxN[j]->SetYTitle("#alpha (rad)");
498 outputContainer->Add(fhAnglePairMassLocMaxN[j]) ;
992b14a7 499
500 }
501
502 return outputContainer ;
503
504}
505
992b14a7 506//___________________________________________
507void AliAnaInsideClusterInvariantMass::Init()
508{
509 //Init
510 //Do some checks
511 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
512 printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
513 abort();
514 }
515 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
516 printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
517 abort();
518 }
519
520 if( GetReader()->GetDataType() == AliCaloTrackReader::kMC ){
521 printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
522 abort();
523
524 }
525
526}
527
528//_____________________________________________________
529void AliAnaInsideClusterInvariantMass::InitParameters()
530{
531 //Initialize the parameters of the analysis.
532 AddToHistogramsName("AnaPi0InsideClusterInvariantMass_");
533
534 fCalorimeter = "EMCAL" ;
29ca9cad 535
71e3889f 536 fM02MinCut = 0.26 ;
537 fM02MaxCut = 10 ;
538
992b14a7 539 fMinNCells = 4 ;
243c2909 540
29ca9cad 541 fMassEtaMin = 0.4;
542 fMassEtaMax = 0.6;
243c2909 543
29ca9cad 544 fMassPi0Min = 0.08;
545 fMassPi0Max = 0.20;
243c2909 546
29ca9cad 547 fMassConMin = 0.0;
548 fMassConMax = 0.05;
243c2909 549
992b14a7 550}
551
552
553//__________________________________________________________________
554void AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
555{
556 //Search for pi0 in fCalorimeter with shower shape analysis
557
558 TObjArray * pl = 0x0;
559 AliVCaloCells* cells = 0x0;
560
561 //Select the Calorimeter of the photon
562 if(fCalorimeter == "PHOS"){
563 pl = GetPHOSClusters();
564 cells = GetPHOSCells();
565 }
566 else if (fCalorimeter == "EMCAL"){
567 pl = GetEMCALClusters();
568 cells = GetEMCALCells();
569 }
570
571 if(!pl || !cells) {
572 Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
573 return;
574 }
575
576 if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
577
578 for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++){
579 AliVCluster * cluster = (AliVCluster*) (pl->At(icluster));
580
581 // Study clusters with large shape parameter
582 Float_t en = cluster->E();
583 Float_t l0 = cluster->GetM02();
584 Int_t nc = cluster->GetNCells();
243c2909 585
992b14a7 586 //If too small or big E or low number of cells, skip it
e23a0471 587 if( en < GetMinEnergy() || en > GetMaxEnergy() || nc < fMinNCells) continue ;
992b14a7 588
589 Int_t absId1 = -1; Int_t absId2 = -1;
590 Int_t *absIdList = new Int_t [nc];
591 Float_t *maxEList = new Float_t[nc];
dbba06ca 592 Int_t nMax = GetCaloUtils()->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
5c46c992 593 Bool_t matched = IsTrackMatched(cluster,GetReader()->GetInputEvent());
992b14a7 594
5c46c992 595 if (nMax <= 0)
596 {
29ca9cad 597 printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found!\n");
5c46c992 598
cab316a8 599 /*
29ca9cad 600 for(Int_t iDigit = 0; iDigit < cluster->GetNCells(); iDigit++ ) {
29ca9cad 601 Float_t ec = cells->GetCellAmplitude(cluster->GetCellsAbsId()[iDigit]);
dbba06ca 602 GetCaloUtils()->RecalibrateCellAmplitude(ec,fCalorimeter,cluster->GetCellsAbsId()[iDigit]);
29ca9cad 603 printf("iDigit %d, absId %d, Ecell %f\n",iDigit,cluster->GetCellsAbsId()[iDigit], ec);
29ca9cad 604 }
cab316a8 605 */
5c46c992 606
c5693f62 607 delete [] absIdList ;
608 delete [] maxEList ;
992b14a7 609 return;
610 }
611
5c46c992 612 fhNLocMax[0][matched]->Fill(en,nMax);
29ca9cad 613 for(Int_t imax = 0; imax < nMax; imax++)
614 {
5c46c992 615 fhNLocMaxEMax [0][matched]->Fill(maxEList[imax] ,nMax);
616 fhNLocMaxEFrac[0][matched]->Fill(maxEList[imax]/en,nMax);
29ca9cad 617 }
618
992b14a7 619
5c46c992 620 if ( nMax == 1 ) { fhM02NLocMax1[0][matched]->Fill(en,l0) ; fhNCellNLocMax1[0][matched]->Fill(en,nc) ; }
621 else if( nMax == 2 ) { fhM02NLocMax2[0][matched]->Fill(en,l0) ; fhNCellNLocMax2[0][matched]->Fill(en,nc) ; }
622 else if( nMax >= 3 ) { fhM02NLocMaxN[0][matched]->Fill(en,l0) ; fhNCellNLocMaxN[0][matched]->Fill(en,nc) ; }
992b14a7 623 else printf("N max smaller than 1 -> %d \n",nMax);
5c46c992 624
625 Float_t dZ = cluster->GetTrackDz();
626 Float_t dR = cluster->GetTrackDx();
627
628 if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
629 {
630 dR = 2000., dZ = 2000.;
631 GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
632 }
633 //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
634
635 if(TMath::Abs(dR) < 999)
636 {
637 if ( nMax == 1 ) { fhTrackMatchedDEtaLocMax1[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMax1[0]->Fill(en,dR); }
638 else if( nMax == 2 ) { fhTrackMatchedDEtaLocMax2[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMax2[0]->Fill(en,dR); }
639 else if( nMax >= 3 ) { fhTrackMatchedDEtaLocMaxN[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMaxN[0]->Fill(en,dR); }
640 }
992b14a7 641
642 // Play with the MC stack if available
643 // Check origin of the candidates
644 Int_t mcindex = -1;
645 if(IsDataMC()){
646
647 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader(), 0);
648
c5693f62 649 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ) mcindex = kmcPi0;
650 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mcindex = kmcEta;
992b14a7 651 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
c5693f62 652 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPhoton;
992b14a7 653 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
c5693f62 654 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
655 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)) mcindex = kmcElectron;
656 else mcindex = kmcHadron;
992b14a7 657
29ca9cad 658 //GetMCAnalysisUtils()->PrintMCTag(tag);
31ae6d59 659 //printf("\t MC index Assigned %d \n",mcindex);
71e3889f 660 fhNLocMaxNLabel[0] [matched]->Fill(en,nMax);
661 fhNLocMaxNLabel[mcindex][matched]->Fill(cluster->GetNLabels(),nMax);
662
5c46c992 663 fhNLocMax[mcindex][matched]->Fill(en,nMax);
29ca9cad 664 for(Int_t imax = 0; imax < nMax; imax++)
665 {
5c46c992 666 fhNLocMaxEMax [mcindex][matched]->Fill(maxEList[imax] ,nMax);
667 fhNLocMaxEFrac[mcindex][matched]->Fill(maxEList[imax]/en,nMax);
668 }
669 if (nMax == 1 ) { fhM02NLocMax1[mcindex][matched]->Fill(en,l0) ; fhNCellNLocMax1[mcindex][matched]->Fill(en,nc) ; }
670 else if(nMax == 2 ) { fhM02NLocMax2[mcindex][matched]->Fill(en,l0) ; fhNCellNLocMax2[mcindex][matched]->Fill(en,nc) ; }
671 else if(nMax >= 3 ) { fhM02NLocMaxN[mcindex][matched]->Fill(en,l0) ; fhNCellNLocMaxN[mcindex][matched]->Fill(en,nc) ; }
672
673 if(TMath::Abs(dR) < 999)
674 {
675 if ( nMax == 1 ) { fhTrackMatchedDEtaLocMax1[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMax1[mcindex]->Fill(en,dR); }
676 else if( nMax == 2 ) { fhTrackMatchedDEtaLocMax2[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMax2[mcindex]->Fill(en,dR); }
677 else if( nMax >= 3 ) { fhTrackMatchedDEtaLocMaxN[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMaxN[mcindex]->Fill(en,dR); }
29ca9cad 678 }
992b14a7 679
680 }
681
29ca9cad 682 //---------------------------------------------------------------------
683 // From here only if M02 is large, fill histograms or split the cluster
684 //---------------------------------------------------------------------
685
71e3889f 686 if( l0 < fM02MinCut || l0 > fM02MaxCut )
29ca9cad 687 {
37f53c3e 688 delete [] absIdList ;
689 delete [] maxEList ;
690 continue;
691 }
992b14a7 692
5c46c992 693 fhNLocMaxM02Cut[0][matched]->Fill(en,nMax);
694 if(IsDataMC()) fhNLocMaxM02Cut[mcindex][matched]->Fill(en,nMax);
243c2909 695
29ca9cad 696 //---------------------------------------------------------------------
992b14a7 697 // Get the 2 max indeces and do inv mass
29ca9cad 698 //---------------------------------------------------------------------
699
243c2909 700 if ( nMax == 2 ) {
701 absId1 = absIdList[0];
702 absId2 = absIdList[1];
703 }
704 else if( nMax == 1 )
705 {
706
707 absId1 = absIdList[0];
992b14a7 708
992b14a7 709 //Find second highest energy cell
243c2909 710
04c5e581 711 Float_t enmax = 0 ;
992b14a7 712 for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++){
713 Int_t absId = cluster->GetCellsAbsId()[iDigit];
714 if( absId == absId1 ) continue ;
715 Float_t endig = cells->GetCellAmplitude(absId);
dbba06ca 716 GetCaloUtils()->RecalibrateCellAmplitude(endig,fCalorimeter,absId);
992b14a7 717 if(endig > enmax) {
718 enmax = endig ;
719 absId2 = absId ;
720 }
992b14a7 721 }// cell loop
243c2909 722 }// 1 maxima
723 else { // n max > 2
724 // loop on maxima, find 2 highest
992b14a7 725
243c2909 726 // First max
727 Float_t enmax = 0 ;
992b14a7 728 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
729 Float_t endig = maxEList[iDigit];
730 if(endig > enmax) {
243c2909 731 enmax = endig ;
732 absId1 = absIdList[iDigit];
733 }
734 }// first maxima loop
735
736 // Second max
737 Float_t enmax2 = 0;
738 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
739 if(absIdList[iDigit]==absId1) continue;
740 Float_t endig = maxEList[iDigit];
741 if(endig > enmax2) {
742 enmax2 = endig ;
992b14a7 743 absId2 = absIdList[iDigit];
744 }
243c2909 745 }// second maxima loop
746
747 } // n local maxima > 2
29ca9cad 748
749 //---------------------------------------------------------------------
750 // Split the cluster energy in 2, around the highest 2 local maxima
751 //---------------------------------------------------------------------
752
5c46c992 753 //Float_t en1 = 0, en2 = 0;
754 //SplitEnergy(absId1,absId2,cluster, cells, en1, en2, nMax /*absIdList, maxEList,*/);
755
756 AliAODCaloCluster *cluster1 = new AliAODCaloCluster(0, 0,NULL,0.,NULL,NULL,1,0);
757 AliAODCaloCluster *cluster2 = new AliAODCaloCluster(1, 0,NULL,0.,NULL,NULL,1,0);
758
759 SplitEnergy(absId1,absId2,cluster, cells, cluster1, cluster2, nMax /*absIdList, maxEList,*/);
760
29ca9cad 761 //---------------------------------------------------------------------
762 // Get mass of pair of clusters
763 //---------------------------------------------------------------------
764
765 // First set position of cluster as local maxima position,
766 // assign splitted energy to calculate momentum
992b14a7 767
5c46c992 768 //TLorentzVector cellMom1 = GetCellMomentum(absId1, en1, cells);
769 //TLorentzVector cellMom2 = GetCellMomentum(absId2, en2, cells);
29ca9cad 770
5c46c992 771 TLorentzVector cellMom1;
772 TLorentzVector cellMom2;
773
774 cluster1->GetMomentum(cellMom1,GetVertex(0));
775 cluster2->GetMomentum(cellMom2,GetVertex(0));
992b14a7 776
5c46c992 777 Float_t mass = (cellMom1+cellMom2).M();
778 Double_t angle = cellMom2.Angle(cellMom1.Vect());
779
243c2909 780 if (nMax==1)
781 {
5c46c992 782 fhAnglePairLocMax1[matched]->Fill(en,angle);
783 if( en > 5 )
784 fhAnglePairMassLocMax1[matched]->Fill(mass,angle);
785
786 fhMassNLocMax1[0][matched]->Fill(en,mass);
787 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax1[0][matched]->Fill(en,l0);
788 else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax1[0][matched]->Fill(en,l0);
789 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax1[0][matched]->Fill(en,l0);
243c2909 790 }
5c46c992 791 else if(nMax==2)
792 {
793 fhAnglePairLocMax2[matched]->Fill(en,angle);
794 if( en > 5 )
795 fhAnglePairMassLocMax2[matched]->Fill(mass,angle);
796
797 fhMassNLocMax2[0][matched]->Fill(en,mass);
798 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax2[0][matched]->Fill(en,l0);
799 else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax2[0][matched]->Fill(en,l0);
800 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax2[0][matched]->Fill(en,l0);
243c2909 801 }
5c46c992 802 else if(nMax >2)
803 {
804 fhAnglePairLocMaxN[matched]->Fill(en,angle);
805 if( en > 5 )
806 fhAnglePairMassLocMaxN[matched]->Fill(mass,angle);
807
808 fhMassNLocMaxN[0][matched]->Fill(en,mass);
809 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMaxN[0][matched]->Fill(en,l0);
810 else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMaxN[0][matched]->Fill(en,l0);
811 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMaxN[0][matched]->Fill(en,l0);
243c2909 812 }
813
992b14a7 814 if(IsDataMC()){
815
243c2909 816 if (nMax==1)
817 {
5c46c992 818 fhMassNLocMax1[mcindex][matched]->Fill(en,mass);
819 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax1[mcindex][matched]->Fill(en,l0);
820 else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax1[mcindex][matched]->Fill(en,l0);
821 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax1[mcindex][matched]->Fill(en,l0);
243c2909 822 }
823 else if(nMax==2) {
5c46c992 824 fhMassNLocMax2[mcindex][matched]->Fill(en,mass);
825 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax2[mcindex][matched]->Fill(en,l0);
826 else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax2[mcindex][matched]->Fill(en,l0);
827 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax2[mcindex][matched]->Fill(en,l0);
243c2909 828 }
829 else if(nMax >2) {
5c46c992 830 fhMassNLocMaxN[mcindex][matched]->Fill(en,mass);
831 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMaxN[mcindex][matched]->Fill(en,l0);
832 else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMaxN[mcindex][matched]->Fill(en,l0);
833 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMaxN[mcindex][matched]->Fill(en,l0);
243c2909 834 }
992b14a7 835
836 }//Work with MC truth first
5c46c992 837
838 delete cluster1 ;
839 delete cluster2 ;
c5693f62 840 delete [] absIdList ;
841 delete [] maxEList ;
992b14a7 842
843 }//loop
844
845 if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");
e23a0471 846
992b14a7 847}
848
849//______________________________________________________________________
850void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
851{
852 //Print some relevant parameters set for the analysis
853 if(! opt)
854 return;
855
856 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 857 AliAnaCaloTrackCorrBaseClass::Print("");
243c2909 858 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
dbba06ca 859 printf("Loc. Max. E > %2.2f\n", GetCaloUtils()->GetLocalMaximaCutE());
860 printf("Loc. Max. E Diff > %2.2f\n", GetCaloUtils()->GetLocalMaximaCutEDiff());
71e3889f 861 printf("%2.2f < lambda_0^2 <%2.2f \n",fM02MinCut,fM02MaxCut);
243c2909 862 printf("pi0 : %2.2f<m<%2.2f \n", fMassPi0Min,fMassPi0Max);
863 printf("eta : %2.2f<m<%2.2f \n", fMassEtaMin,fMassEtaMax);
864 printf("conv: %2.2f<m<%2.2f \n", fMassConMin,fMassConMax);
865
992b14a7 866 printf(" \n") ;
867
868}
869
dbba06ca 870
992b14a7 871
c5693f62 872//________________________________________________________________________________________
873void AliAnaInsideClusterInvariantMass::SplitEnergy(const Int_t absId1, const Int_t absId2,
243c2909 874 AliVCluster* cluster,
875 AliVCaloCells* cells,
5c46c992 876 //Float_t & e1, Float_t & e2,
877 AliAODCaloCluster* cluster1,
878 AliAODCaloCluster* cluster2,
29ca9cad 879 const Int_t nMax)
992b14a7 880{
5c46c992 881
29ca9cad 882 // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
883 // maxima are too close and have common cells, split the energy between the 2
884
71e3889f 885 TH2F* hClusterMap = 0 ;
886 TH2F* hClusterLocMax = 0 ;
887 TH2F* hCluster1 = 0 ;
888 TH2F* hCluster2 = 0 ;
889
890 if(fPlotCluster)
891 {
892 hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
893 hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
894 hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
895 hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
896 hClusterMap ->SetXTitle("column");
897 hClusterMap ->SetYTitle("row");
898 hClusterLocMax ->SetXTitle("column");
899 hClusterLocMax ->SetYTitle("row");
900 hCluster1 ->SetXTitle("column");
901 hCluster1 ->SetYTitle("row");
902 hCluster2 ->SetXTitle("column");
903 hCluster2 ->SetYTitle("row");
904 }
29ca9cad 905
243c2909 906 const Int_t ncells = cluster->GetNCells();
29ca9cad 907 Int_t absIdList[ncells];
908
5c46c992 909 Float_t e1 = 0, e2 = 0 ;
71e3889f 910 Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
29ca9cad 911 Float_t eCluster = 0;
71e3889f 912 Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
243c2909 913 for(Int_t iDigit = 0; iDigit < ncells; iDigit++ ) {
914 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
243c2909 915
71e3889f 916
29ca9cad 917 Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
dbba06ca 918 GetCaloUtils()->RecalibrateCellAmplitude(ec,fCalorimeter, absIdList[iDigit]);
29ca9cad 919 eCluster+=ec;
920
71e3889f 921 if(fPlotCluster)
922 {
923 //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
924 sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList[iDigit], fCalorimeter, icol, irow, iRCU) ;
925 if(icol > maxCol) maxCol = icol;
926 if(icol < minCol) minCol = icol;
927 if(irow > maxRow) maxRow = irow;
928 if(irow < minRow) minRow = irow;
929 hClusterMap->Fill(icol,irow,ec);
930 }
29ca9cad 931
932 }
933
243c2909 934 // Init counters and variables
5c46c992 935 Int_t ncells1 = 1 ;
936 UShort_t absIdList1[9] ;
937 Double_t fracList1 [9] ;
938 absIdList1[0] = absId1 ;
939 fracList1 [0] = 1. ;
243c2909 940
941 Float_t ecell1 = cells->GetCellAmplitude(absId1);
dbba06ca 942 GetCaloUtils()->RecalibrateCellAmplitude(ecell1, fCalorimeter, absId1);
243c2909 943 e1 = ecell1;
944
5c46c992 945 Int_t ncells2 = 1 ;
946 UShort_t absIdList2[9] ;
947 Double_t fracList2 [9] ;
948 absIdList2[0] = absId2 ;
949 fracList2 [0] = 1. ;
950
243c2909 951 Float_t ecell2 = cells->GetCellAmplitude(absId2);
dbba06ca 952 GetCaloUtils()->RecalibrateCellAmplitude(ecell2, fCalorimeter, absId2);
243c2909 953 e2 = ecell2;
29ca9cad 954
71e3889f 955 if(fPlotCluster)
956 {
957 Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
958 sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId1, fCalorimeter, icol1, irow1, iRCU) ;
959 hClusterLocMax->Fill(icol1,irow1,ecell1);
960 sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId2, fCalorimeter, icol2, irow2, iRCU) ;
961 hClusterLocMax->Fill(icol2,irow2,ecell2);
962 }
29ca9cad 963
964 // Very rough way to share the cluster energy
965 Float_t eRemain = (eCluster-ecell1-ecell2)/2;
966 Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
967 Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
968
969 for(Int_t iDigit = 0; iDigit < ncells; iDigit++){
970 Int_t absId = absIdList[iDigit];
971
972 if(absId==absId1 || absId==absId2 || absId < 0) continue;
973
974 Float_t ecell = cells->GetCellAmplitude(absId);
dbba06ca 975 GetCaloUtils()->RecalibrateCellAmplitude(ecell, fCalorimeter, absId);
29ca9cad 976
5c46c992 977 if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId1,absId ))
978 {
979 absIdList1[ncells1]= absId;
29ca9cad 980
5c46c992 981 if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId2,absId ))
982 {
983 fracList1[ncells1] = shareFraction1;
984 e1 += ecell*shareFraction1;
985 }
986 else
987 {
988 fracList1[ncells1] = 1.;
989 e1 += ecell;
990 }
991
992 ncells1++;
993
29ca9cad 994 } // neigbour to cell1
243c2909 995
5c46c992 996 if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId2,absId ))
997 {
998 absIdList2[ncells2]= absId;
29ca9cad 999
5c46c992 1000 if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId1,absId ))
1001 {
1002 fracList2[ncells2] = shareFraction2;
29ca9cad 1003 e2 += ecell*shareFraction2;
5c46c992 1004 }
1005 else
1006 {
1007 fracList2[ncells2] = 1.;
29ca9cad 1008 e2 += ecell;
5c46c992 1009 }
1010
1011 ncells2++;
1012
29ca9cad 1013 } // neigbour to cell2
5c46c992 1014
29ca9cad 1015 }
1016
1017 if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::SplitEnergy() - n Local Max %d, Cluster energy = %f, Ecell1 = %f, Ecell2 = %f, Enew1 = %f, Enew2 = %f, Remain %f, \n ncells %d, ncells1 %d, ncells2 %d, f1 %f, f2 %f, sum f12 = %f \n",
1018 nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
5c46c992 1019
1020 cluster1->SetE(e1);
1021 cluster2->SetE(e2);
1022
1023 cluster1->SetNCells(ncells1);
1024 cluster2->SetNCells(ncells2);
1025
1026 cluster1->SetCellsAbsId(absIdList1);
1027 cluster2->SetCellsAbsId(absIdList2);
1028
1029 cluster1->SetCellsAmplitudeFraction(fracList1);
1030 cluster2->SetCellsAmplitudeFraction(fracList2);
1031
1032 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1033 GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1034
29ca9cad 1035
71e3889f 1036 if(fPlotCluster)
243c2909 1037 {
71e3889f 1038 //printf("Cells of cluster1: ");
1039 for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
1040 {
1041 //printf(" %d ",absIdList1[iDigit]);
1042
1043 sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList1[iDigit], fCalorimeter, icol, irow, iRCU) ;
1044
1045 if( GetCaloUtils()->AreNeighbours(fCalorimeter, absId2,absIdList1[iDigit]) )
1046 hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit])*shareFraction1);
1047 else
1048 hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
1049 }
1050
1051 //printf(" \n ");
1052 //printf("Cells of cluster2: ");
1053
1054 for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
1055 {
1056 //printf(" %d ",absIdList2[iDigit]);
1057
1058 sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList2[iDigit], fCalorimeter, icol, irow, iRCU) ;
1059 if( GetCaloUtils()->AreNeighbours(fCalorimeter, absId1,absIdList2[iDigit]) )
1060 hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit])*shareFraction2);
1061 else
1062 hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
1063
1064 }
1065 //printf(" \n ");
243c2909 1066
71e3889f 1067 gStyle->SetPadRightMargin(0.1);
1068 gStyle->SetPadLeftMargin(0.1);
1069 gStyle->SetOptStat(0);
1070 gStyle->SetOptFit(000000);
243c2909 1071
71e3889f 1072 if(maxCol-minCol > maxRow-minRow)
1073 {
1074 maxRow+= maxCol-minCol;
1075 }
29ca9cad 1076 else
71e3889f 1077 {
1078 maxCol+= maxRow-minRow;
1079 }
29ca9cad 1080
71e3889f 1081 TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1082 c->Divide(2,2);
1083 c->cd(1);
1084 gPad->SetGridy();
1085 gPad->SetGridx();
1086 hClusterMap ->SetAxisRange(minCol, maxCol,"X");
1087 hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
1088 hClusterMap ->Draw("colz");
1089 c->cd(2);
1090 gPad->SetGridy();
1091 gPad->SetGridx();
1092 hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1093 hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
1094 hClusterLocMax ->Draw("colz");
1095 c->cd(3);
1096 gPad->SetGridy();
1097 gPad->SetGridx();
1098 hCluster1 ->SetAxisRange(minCol, maxCol,"X");
1099 hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
1100 hCluster1 ->Draw("colz");
1101 c->cd(4);
1102 gPad->SetGridy();
1103 gPad->SetGridx();
1104 hCluster2 ->SetAxisRange(minCol, maxCol,"X");
1105 hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
1106 hCluster2 ->Draw("colz");
1107
1108 if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",GetEventNumber(),cluster->E(),nMax,ncells1,ncells2));
1109
1110 delete c;
1111 delete hClusterMap;
1112 delete hClusterLocMax;
1113 delete hCluster1;
1114 delete hCluster2;
29ca9cad 1115 }
992b14a7 1116}
1117