]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/PartCorrDep/AliAnaInsideClusterInvariantMass.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix for coverity 17980
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / 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>
992b14a7 33
34// --- Analysis system ---
35#include "AliAnaInsideClusterInvariantMass.h"
36#include "AliCaloTrackReader.h"
37#include "AliMCAnalysisUtils.h"
38#include "AliStack.h"
39#include "AliFiducialCut.h"
40#include "TParticle.h"
41#include "AliVCluster.h"
42#include "AliAODEvent.h"
43#include "AliAODMCParticle.h"
44#include "AliEMCALGeoParams.h"
45
c5693f62 46// --- Detectors ---
47//#include "AliPHOSGeoUtils.h"
48#include "AliEMCALGeometry.h"
49
992b14a7 50ClassImp(AliAnaInsideClusterInvariantMass)
51
52//__________________________________________________________________
53AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
745913ae 54 AliAnaCaloTrackCorrBaseClass(),
992b14a7 55 fCalorimeter(""),
243c2909 56 fM02Cut(0), fMinNCells(0),
57 fMassEtaMin(0), fMassEtaMax(0),
58 fMassPi0Min(0), fMassPi0Max(0),
59 fMassConMin(0), fMassConMax(0)
992b14a7 60{
61 //default ctor
62
63 // Init array of histograms
64 for(Int_t i = 0; i < 7; i++){
65 fhMassNLocMax1[i] = 0;
66 fhMassNLocMax2[i] = 0;
67 fhMassNLocMaxN[i] = 0;
68 fhNLocMax[i] = 0;
69 fhNLocMaxM02Cut[i] = 0;
70 fhM02NLocMax1[i] = 0;
71 fhM02NLocMax2[i] = 0;
72 fhM02NLocMaxN[i] = 0;
73 fhNCellNLocMax1[i] = 0;
74 fhNCellNLocMax2[i] = 0;
75 fhNCellNLocMaxN[i] = 0;
243c2909 76 fhM02Pi0LocMax1[i] = 0;
77 fhM02EtaLocMax1[i] = 0;
78 fhM02ConLocMax1[i] = 0;
79 fhM02Pi0LocMax2[i] = 0;
80 fhM02EtaLocMax2[i] = 0;
81 fhM02ConLocMax2[i] = 0;
82 fhM02Pi0LocMaxN[i] = 0;
83 fhM02EtaLocMaxN[i] = 0;
84 fhM02ConLocMaxN[i] = 0;
85
992b14a7 86 }
87
88 InitParameters();
89
90}
91
92//_______________________________________________________________
93TObjString * AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
94{
95 //Save parameters used for analysis
96 TString parList ; //this will be list of parameters used for this analysis.
97 const Int_t buffersize = 255;
98 char onePar[buffersize] ;
99
100 snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
101 parList+=onePar ;
102
243c2909 103 snprintf(onePar,buffersize,"Calorimeter: %s\n", fCalorimeter.Data()) ;
992b14a7 104 parList+=onePar ;
243c2909 105 snprintf(onePar,buffersize,"fM02Cut =%f \n" , fM02Cut) ;
992b14a7 106 parList+=onePar ;
243c2909 107 snprintf(onePar,buffersize,"fMinNCells =%f \n", fMinNCells) ;
992b14a7 108 parList+=onePar ;
243c2909 109 snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f\n", fMassPi0Min,fMassPi0Max);
110 parList+=onePar ;
111 snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f\n", fMassEtaMin,fMassEtaMax);
112 parList+=onePar ;
113 snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f\n", fMassConMin,fMassConMax);
114 parList+=onePar ;
115
992b14a7 116 return new TObjString(parList) ;
117
118}
119
120//____________________________________________________________________________________________________
121Bool_t AliAnaInsideClusterInvariantMass::AreNeighbours( const Int_t absId1, const Int_t absId2 ) const
122{
123 // Tells if (true) or not (false) two digits are neighbours
124 // A neighbour is defined as being two digits which share a corner
125
126 Bool_t areNeighbours = kFALSE ;
127 Int_t nSupMod =0, nModule =0, nIphi =0, nIeta =0;
128 Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0;
129 Int_t relid1[2], relid2[2] ;
130 Int_t rowdiff=0, coldiff=0;
131
132 areNeighbours = kFALSE ;
133
134 GetEMCALGeometry()->GetCellIndex(absId1, nSupMod,nModule,nIphi,nIeta);
135 GetEMCALGeometry()->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, relid1[0],relid1[1]);
136
137 GetEMCALGeometry()->GetCellIndex(absId2, nSupMod1,nModule1,nIphi1,nIeta1);
138 GetEMCALGeometry()->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, relid2[0],relid2[1]);
139
140 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
141 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
142 if(nSupMod1!=nSupMod){
143 if(nSupMod1%2) relid1[1]+=AliEMCALGeoParams::fgkEMCALCols;
144 else relid2[1]+=AliEMCALGeoParams::fgkEMCALCols;
145 }
146
147 rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;
148 coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;
149
150 if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
151 areNeighbours = kTRUE ;
152
153 return areNeighbours;
154}
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
176 if(en<=0)
177 {
178 en = cells->GetCellAmplitude(absId);
179 RecalibrateCellAmplitude(en,absId);
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
202 TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
203 TString pname[] ={"","Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
204
205 Int_t n = 1;
206
207 if(IsDataMC()) n = 7;
208
243c2909 209 Int_t nMaxBins = 10;
210
211 for(Int_t i = 0; i < n; i++){
992b14a7 212
213 fhMassNLocMax1[i] = new TH2F(Form("hMassNLocMax1%s",pname[i].Data()),
214 Form("Invariant mass of 2 highest energy cells %s",ptype[i].Data()),
215 nptbins,ptmin,ptmax,mbins,mmin,mmax);
216 fhMassNLocMax1[i]->SetYTitle("M (MeV/c^2)");
217 fhMassNLocMax1[i]->SetXTitle("E (GeV)");
218 outputContainer->Add(fhMassNLocMax1[i]) ;
219
220 fhMassNLocMax2[i] = new TH2F(Form("hMassNLocMax2%s",pname[i].Data()),
221 Form("Invariant mass of 2 local maxima cells %s",ptype[i].Data()),
222 nptbins,ptmin,ptmax,mbins,mmin,mmax);
223 fhMassNLocMax2[i]->SetYTitle("M (MeV/c^2)");
224 fhMassNLocMax2[i]->SetXTitle("E (GeV)");
225 outputContainer->Add(fhMassNLocMax2[i]) ;
226
227 fhMassNLocMaxN[i] = new TH2F(Form("hMassNLocMaxN%s",pname[i].Data()),
243c2909 228 Form("Invariant mass of N>2 local maxima cells %s",ptype[i].Data()),
992b14a7 229 nptbins,ptmin,ptmax,mbins,mmin,mmax);
230 fhMassNLocMaxN[i]->SetYTitle("M (MeV/c^2)");
243c2909 231 fhMassNLocMaxN[i]->SetXTitle("E (GeV)");
992b14a7 232 outputContainer->Add(fhMassNLocMaxN[i]) ;
233
234
235 fhNLocMax[i] = new TH2F(Form("hNLocMax%s",pname[i].Data()),
236 Form("Number of local maxima in cluster %s",ptype[i].Data()),
243c2909 237 nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
992b14a7 238 fhNLocMax[i] ->SetYTitle("N maxima");
239 fhNLocMax[i] ->SetXTitle("E (GeV)");
240 outputContainer->Add(fhNLocMax[i]) ;
241
242 fhNLocMaxM02Cut[i] = new TH2F(Form("hNLocMaxM02Cut%s",pname[i].Data()),
243 Form("Number of local maxima in cluster %s for M02 > %2.2f",ptype[i].Data(),fM02Cut),
243c2909 244 nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
992b14a7 245 fhNLocMaxM02Cut[i]->SetYTitle("N maxima");
246 fhNLocMaxM02Cut[i]->SetXTitle("E (GeV)");
247 outputContainer->Add(fhNLocMaxM02Cut[i]) ;
248
249
250 fhM02NLocMax1[i] = new TH2F(Form("hM02NLocMax1%s",pname[i].Data()),
251 Form("#lambda_{0}^{2} vs E for N max = 1 %s",ptype[i].Data()),
252 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
253 fhM02NLocMax1[i] ->SetYTitle("#lambda_{0}^{2}");
254 fhM02NLocMax1[i] ->SetXTitle("E (GeV)");
255 outputContainer->Add(fhM02NLocMax1[i]) ;
256
257 fhM02NLocMax2[i] = new TH2F(Form("hM02NLocMax2%s",pname[i].Data()),
258 Form("#lambda_{0}^{2} vs E for N max = 2 %s",ptype[i].Data()),
259 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
260 fhM02NLocMax2[i] ->SetYTitle("#lambda_{0}^{2}");
261 fhM02NLocMax2[i] ->SetXTitle("E (GeV)");
262 outputContainer->Add(fhM02NLocMax2[i]) ;
263
264
265 fhM02NLocMaxN[i] = new TH2F(Form("hM02NLocMaxN%s",pname[i].Data()),
266 Form("#lambda_{0}^{2} vs E for N max > 2 %s",ptype[i].Data()),
267 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
268 fhM02NLocMaxN[i] ->SetYTitle("#lambda_{0}^{2}");
269 fhM02NLocMaxN[i] ->SetXTitle("E (GeV)");
270 outputContainer->Add(fhM02NLocMaxN[i]) ;
271
272
273 fhNCellNLocMax1[i] = new TH2F(Form("hNCellNLocMax1%s",pname[i].Data()),
274 Form("#lambda_{0}^{2} vs E for N max = 1 %s",ptype[i].Data()),
275 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
276 fhNCellNLocMax1[i] ->SetYTitle("N cells");
277 fhNCellNLocMax1[i] ->SetXTitle("E (GeV)");
278 outputContainer->Add(fhNCellNLocMax1[i]) ;
279
280 fhNCellNLocMax2[i] = new TH2F(Form("hNCellNLocMax2%s",pname[i].Data()),
281 Form("#lambda_{0}^{2} vs E for N max = 2 %s",ptype[i].Data()),
282 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
283 fhNCellNLocMax2[i] ->SetYTitle("N cells");
284 fhNCellNLocMax2[i] ->SetXTitle("E (GeV)");
285 outputContainer->Add(fhNCellNLocMax2[i]) ;
286
287
288 fhNCellNLocMaxN[i] = new TH2F(Form("hNCellNLocMaxN%s",pname[i].Data()),
289 Form("#lambda_{0}^{2} vs E for N max > 2 %s",ptype[i].Data()),
290 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
291 fhNCellNLocMaxN[i] ->SetYTitle("N cells");
292 fhNCellNLocMaxN[i] ->SetXTitle("E (GeV)");
293 outputContainer->Add(fhNCellNLocMaxN[i]) ;
294
295
243c2909 296 fhM02Pi0LocMax1[i] = new TH2F(Form("hM02Pi0LocMax1%s",pname[i].Data()),
297 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()),
992b14a7 298 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
243c2909 299 fhM02Pi0LocMax1[i] ->SetYTitle("#lambda_{0}^{2}");
300 fhM02Pi0LocMax1[i] ->SetXTitle("E (GeV)");
301 outputContainer->Add(fhM02Pi0LocMax1[i]) ;
992b14a7 302
243c2909 303 fhM02EtaLocMax1[i] = new TH2F(Form("hM02EtaLocMax1%s",pname[i].Data()),
304 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()),
992b14a7 305 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
243c2909 306 fhM02EtaLocMax1[i] ->SetYTitle("#lambda_{0}^{2}");
307 fhM02EtaLocMax1[i] ->SetXTitle("E (GeV)");
308 outputContainer->Add(fhM02EtaLocMax1[i]) ;
992b14a7 309
243c2909 310 fhM02ConLocMax1[i] = new TH2F(Form("hM02ConLocMax1%s",pname[i].Data()),
311 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()),
992b14a7 312 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
243c2909 313 fhM02ConLocMax1[i] ->SetYTitle("#lambda_{0}^{2}");
314 fhM02ConLocMax1[i] ->SetXTitle("E (GeV)");
315 outputContainer->Add(fhM02ConLocMax1[i]) ;
316
317 fhM02Pi0LocMax2[i] = new TH2F(Form("hM02Pi0LocMax2%s",pname[i].Data()),
318 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()),
319 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
320 fhM02Pi0LocMax2[i] ->SetYTitle("#lambda_{0}^{2}");
321 fhM02Pi0LocMax2[i] ->SetXTitle("E (GeV)");
322 outputContainer->Add(fhM02Pi0LocMax2[i]) ;
323
324 fhM02EtaLocMax2[i] = new TH2F(Form("hM02EtaLocMax2%s",pname[i].Data()),
325 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()),
326 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
327 fhM02EtaLocMax2[i] ->SetYTitle("#lambda_{0}^{2}");
328 fhM02EtaLocMax2[i] ->SetXTitle("E (GeV)");
329 outputContainer->Add(fhM02EtaLocMax2[i]) ;
330
331 fhM02ConLocMax2[i] = new TH2F(Form("hM02ConLocMax2%s",pname[i].Data()),
332 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()),
333 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
334 fhM02ConLocMax2[i] ->SetYTitle("#lambda_{0}^{2}");
335 fhM02ConLocMax2[i] ->SetXTitle("E (GeV)");
336 outputContainer->Add(fhM02ConLocMax2[i]) ;
337
338 fhM02Pi0LocMaxN[i] = new TH2F(Form("hM02Pi0LocMaxN%s",pname[i].Data()),
339 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()),
340 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
341 fhM02Pi0LocMaxN[i] ->SetYTitle("#lambda_{0}^{2}");
342 fhM02Pi0LocMaxN[i] ->SetXTitle("E (GeV)");
343 outputContainer->Add(fhM02Pi0LocMaxN[i]) ;
344
345 fhM02EtaLocMaxN[i] = new TH2F(Form("hM02EtaLocMaxN%s",pname[i].Data()),
346 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()),
347 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
348 fhM02EtaLocMaxN[i] ->SetYTitle("#lambda_{0}^{2}");
349 fhM02EtaLocMaxN[i] ->SetXTitle("E (GeV)");
350 outputContainer->Add(fhM02EtaLocMaxN[i]) ;
351
352 fhM02ConLocMaxN[i] = new TH2F(Form("hM02ConLocMaxN%s",pname[i].Data()),
353 Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",fMassConMin,fMassConMax,ptype[i].Data()),
354 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
355 fhM02ConLocMaxN[i] ->SetYTitle("#lambda_{0}^{2}");
356 fhM02ConLocMaxN[i] ->SetXTitle("E (GeV)");
357 outputContainer->Add(fhM02ConLocMaxN[i]) ;
992b14a7 358
359 }
360
361 return outputContainer ;
362
363}
364
365//________________________________________________________________________________________________________
366Int_t AliAnaInsideClusterInvariantMass::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
367 Int_t *absIdList, Float_t *maxEList)
368{
369 // Find local maxima in cluster
370
371 Float_t locMaxCut = 0; // not used for the moment
372
373 Int_t iDigitN = 0 ;
374 Int_t iDigit = 0 ;
375 Int_t absId1 = -1 ;
376 Int_t absId2 = -1 ;
377 const Int_t nCells = cluster->GetNCells();
243c2909 378
379 //printf("cluster : ncells %d \n",nCells);
992b14a7 380 for(iDigit = 0; iDigit < nCells ; iDigit++){
381 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
243c2909 382 /*
383 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
384 RecalibrateCellAmplitude(en,absIdList[iDigit]);
385 Int_t icol = -1, irow = -1, iRCU = -1;
386 Int_t sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList[iDigit], fCalorimeter, icol, irow, iRCU) ;
387
388 printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
389 */
992b14a7 390 }
391
243c2909 392
992b14a7 393 for(iDigit = 0 ; iDigit < nCells; iDigit++) {
394 if(absIdList[iDigit]>=0) {
395
396 absId1 = absIdList[iDigit] ;
243c2909 397 //printf("%d : absID111 %d, %s\n",iDigit, absId1,fCalorimeter.Data());
398
992b14a7 399 Float_t en1 = cells->GetCellAmplitude(absId1);
400 RecalibrateCellAmplitude(en1,absId1);
401
402 for(iDigitN = 0; iDigitN < nCells; iDigitN++) {
243c2909 403
404 absId2 = absIdList[iDigitN] ;
405
406 if(absId2==-1) continue;
407
408 //printf("\t %d : absID222 %d, %s\n",iDigitN, absId2,fCalorimeter.Data());
992b14a7 409
410 Float_t en2 = cells->GetCellAmplitude(absId2);
411 RecalibrateCellAmplitude(en2,absId2);
412
413 if ( AreNeighbours(absId1, absId2) ) {
414
415 if (en1 > en2 ) {
416 absIdList[iDigitN] = -1 ;
243c2909 417 //printf("\t \t indexN %d not local max\n",iDigitN);
992b14a7 418 // but may be digit too is not local max ?
243c2909 419 if(en1 < en2 + locMaxCut) {
420 //printf("\t \t index %d not local max cause locMaxCut\n",iDigit);
992b14a7 421 absIdList[iDigit] = -1 ;
243c2909 422 }
992b14a7 423 }
424 else {
425 absIdList[iDigit] = -1 ;
243c2909 426 //printf("\t \t index %d not local max\n",iDigitN);
992b14a7 427 // but may be digitN too is not local max ?
428 if(en1 > en2 - locMaxCut)
243c2909 429 {
992b14a7 430 absIdList[iDigitN] = -1 ;
243c2909 431 //printf("\t \t indexN %d not local max cause locMaxCut\n",iDigit);
432 }
992b14a7 433 }
434 } // if Areneighbours
435 } // while digitN
436 } // slot not empty
437 } // while digit
438
439 iDigitN = 0 ;
440 for(iDigit = 0; iDigit < nCells; iDigit++) {
441 if(absIdList[iDigit]>=0 ){
442 absIdList[iDigitN] = absIdList[iDigit] ;
443 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
444 RecalibrateCellAmplitude(en,absIdList[iDigit]);
445 if(en < 0.1) continue; // Maxima only with seed energy at least
446 maxEList[iDigitN] = en ;
447 //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
448 iDigitN++ ;
449 }
450 }
451
243c2909 452 //printf("N maxima %d \n",iDigitN);
453 //for(Int_t imax = 0; imax < iDigitN; imax++) printf("imax %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
454
992b14a7 455 return iDigitN ;
456
457}
458
459
460//___________________________________________
461void AliAnaInsideClusterInvariantMass::Init()
462{
463 //Init
464 //Do some checks
465 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
466 printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
467 abort();
468 }
469 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
470 printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
471 abort();
472 }
473
474 if( GetReader()->GetDataType() == AliCaloTrackReader::kMC ){
475 printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
476 abort();
477
478 }
479
480}
481
482//_____________________________________________________
483void AliAnaInsideClusterInvariantMass::InitParameters()
484{
485 //Initialize the parameters of the analysis.
486 AddToHistogramsName("AnaPi0InsideClusterInvariantMass_");
487
488 fCalorimeter = "EMCAL" ;
489 fM02Cut = 0.26 ;
490 fMinNCells = 4 ;
243c2909 491
492 fMassEtaMin = 0.4;
493 fMassEtaMax = 0.6;
494
495 fMassPi0Min = 0.08;
496 fMassPi0Max = 0.20;
497
498 fMassConMin = 0.0;
499 fMassConMax = 0.05;
500
992b14a7 501}
502
503
504//__________________________________________________________________
505void AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
506{
507 //Search for pi0 in fCalorimeter with shower shape analysis
508
509 TObjArray * pl = 0x0;
510 AliVCaloCells* cells = 0x0;
511
512 //Select the Calorimeter of the photon
513 if(fCalorimeter == "PHOS"){
514 pl = GetPHOSClusters();
515 cells = GetPHOSCells();
516 }
517 else if (fCalorimeter == "EMCAL"){
518 pl = GetEMCALClusters();
519 cells = GetEMCALCells();
520 }
521
522 if(!pl || !cells) {
523 Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
524 return;
525 }
526
527 if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
528
529 for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++){
530 AliVCluster * cluster = (AliVCluster*) (pl->At(icluster));
531
532 // Study clusters with large shape parameter
533 Float_t en = cluster->E();
534 Float_t l0 = cluster->GetM02();
535 Int_t nc = cluster->GetNCells();
243c2909 536
992b14a7 537 //If too small or big E or low number of cells, skip it
538 if( ( en < GetMinEnergy() || en > GetMaxEnergy() ) && nc < fMinNCells) continue ;
539
540 Int_t absId1 = -1; Int_t absId2 = -1;
541 Int_t *absIdList = new Int_t [nc];
542 Float_t *maxEList = new Float_t[nc];
543 Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
544
545 if (nMax <= 0) {
546 printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found!");
c5693f62 547 delete [] absIdList ;
548 delete [] maxEList ;
992b14a7 549 return;
550 }
551
552 fhNLocMax[0]->Fill(en,nMax);
553
554 if ( nMax == 1 ) { fhM02NLocMax1[0]->Fill(en,l0) ; fhNCellNLocMax1[0]->Fill(en,nc) ; }
555 else if( nMax == 2 ) { fhM02NLocMax2[0]->Fill(en,l0) ; fhNCellNLocMax2[0]->Fill(en,nc) ; }
556 else if( nMax >= 3 ) { fhM02NLocMaxN[0]->Fill(en,l0) ; fhNCellNLocMaxN[0]->Fill(en,nc) ; }
557 else printf("N max smaller than 1 -> %d \n",nMax);
558
559 // Play with the MC stack if available
560 // Check origin of the candidates
561 Int_t mcindex = -1;
562 if(IsDataMC()){
563
564 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader(), 0);
565
c5693f62 566 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ) mcindex = kmcPi0;
567 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mcindex = kmcEta;
992b14a7 568 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
c5693f62 569 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPhoton;
992b14a7 570 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
c5693f62 571 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
572 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)) mcindex = kmcElectron;
573 else mcindex = kmcHadron;
992b14a7 574
575/* printf("AliAnaInsideClusterInvariantMass::FillAnalysisMakeHistograms() - tag %d, photon %d, prompt %d, frag %d, isr %d, pi0 decay %d, eta decay %d, other decay %d \n conv %d, pi0 %d, hadron %d, electron %d, unk %d, muon %d,pion %d, proton %d, neutron %d, kaon %d, antiproton %d, antineutron %d, bad %d\n",tag,
576 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton),
577 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt),
578 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation),
579 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR),
580 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay),
581 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay),
582 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay),
583
584 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
585 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
586 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
587 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
588 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
589 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon),
590 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
591 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton),
592 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
593 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon),
594 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton),
595 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
596 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCBadLabel)
597);
598*/
599
600 fhNLocMax[mcindex]->Fill(en,nMax);
601
602 if (nMax == 1 ) { fhM02NLocMax1[mcindex]->Fill(en,l0) ; fhNCellNLocMax1[mcindex]->Fill(en,nc) ; }
603 else if(nMax == 2 ) { fhM02NLocMax2[mcindex]->Fill(en,l0) ; fhNCellNLocMax2[mcindex]->Fill(en,nc) ; }
604 else if(nMax >= 3 ) { fhM02NLocMaxN[mcindex]->Fill(en,l0) ; fhNCellNLocMaxN[mcindex]->Fill(en,nc) ; }
605
606 }
607
37f53c3e 608 if( l0 < fM02Cut) {
609 delete [] absIdList ;
610 delete [] maxEList ;
611 continue;
612 }
992b14a7 613
243c2909 614 fhNLocMaxM02Cut[0]->Fill(en,nMax);
615 if(IsDataMC()) fhNLocMaxM02Cut[mcindex]->Fill(en,nMax);
616
992b14a7 617 // Get the 2 max indeces and do inv mass
618
243c2909 619 if ( nMax == 2 ) {
620 absId1 = absIdList[0];
621 absId2 = absIdList[1];
622 }
623 else if( nMax == 1 )
624 {
625
626 absId1 = absIdList[0];
992b14a7 627
992b14a7 628 //Find second highest energy cell
243c2909 629
04c5e581 630 Float_t enmax = 0 ;
992b14a7 631 for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++){
632 Int_t absId = cluster->GetCellsAbsId()[iDigit];
633 if( absId == absId1 ) continue ;
634 Float_t endig = cells->GetCellAmplitude(absId);
635 RecalibrateCellAmplitude(endig,absId);
636 if(endig > enmax) {
637 enmax = endig ;
638 absId2 = absId ;
639 }
992b14a7 640 }// cell loop
243c2909 641 }// 1 maxima
642 else { // n max > 2
643 // loop on maxima, find 2 highest
992b14a7 644
243c2909 645 // First max
646 Float_t enmax = 0 ;
992b14a7 647 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
648 Float_t endig = maxEList[iDigit];
649 if(endig > enmax) {
243c2909 650 enmax = endig ;
651 absId1 = absIdList[iDigit];
652 }
653 }// first maxima loop
654
655 // Second max
656 Float_t enmax2 = 0;
657 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
658 if(absIdList[iDigit]==absId1) continue;
659 Float_t endig = maxEList[iDigit];
660 if(endig > enmax2) {
661 enmax2 = endig ;
992b14a7 662 absId2 = absIdList[iDigit];
663 }
243c2909 664 }// second maxima loop
665
666 } // n local maxima > 2
992b14a7 667
243c2909 668 Float_t en1 = -1, en2 = -1;
669 SplitEnergy(absId1,absId2,cluster, cells,en1,en2);
992b14a7 670
243c2909 671 TLorentzVector cellMom1 = GetCellMomentum(absId1, en1, cells);
672 TLorentzVector cellMom2 = GetCellMomentum(absId2, en2, cells);
992b14a7 673
674 Float_t mass = (cellMom1+cellMom2).M();
675
243c2909 676 if (nMax==1)
677 {
678 fhMassNLocMax1[0]->Fill(en,mass);
679 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax1[0]->Fill(en,l0);
680 else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax1[0]->Fill(en,l0);
681 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax1[0]->Fill(en,l0);
682 }
683 else if(nMax==2) {
684 fhMassNLocMax2[0]->Fill(en,mass);
685 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax2[0]->Fill(en,l0);
686 else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax2[0]->Fill(en,l0);
687 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax2[0]->Fill(en,l0);
688 }
689 else if(nMax >2) {
690 fhMassNLocMaxN[0]->Fill(en,mass);
691 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMaxN[0]->Fill(en,l0);
692 else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMaxN[0]->Fill(en,l0);
693 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMaxN[0]->Fill(en,l0);
694 }
695
992b14a7 696 if(IsDataMC()){
697
243c2909 698 if (nMax==1)
699 {
700 fhMassNLocMax1[mcindex]->Fill(en,mass);
701 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax1[mcindex]->Fill(en,l0);
702 else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax1[mcindex]->Fill(en,l0);
703 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax1[mcindex]->Fill(en,l0);
704 }
705 else if(nMax==2) {
706 fhMassNLocMax2[mcindex]->Fill(en,mass);
707 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax2[mcindex]->Fill(en,l0);
708 else if(mass < fMassPi0Max && mass > fMassPi0Min)fhM02Pi0LocMax2[mcindex]->Fill(en,l0);
709 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax2[mcindex]->Fill(en,l0);
710 }
711 else if(nMax >2) {
712 fhMassNLocMaxN[mcindex]->Fill(en,mass);
713 if (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMaxN[mcindex]->Fill(en,l0);
714 else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMaxN[mcindex]->Fill(en,l0);
715 else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMaxN[mcindex]->Fill(en,l0);
716 }
992b14a7 717
718 }//Work with MC truth first
719
c5693f62 720 delete [] absIdList ;
721 delete [] maxEList ;
992b14a7 722
723 }//loop
724
725 if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");
726
727}
728
729//______________________________________________________________________
730void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
731{
732 //Print some relevant parameters set for the analysis
733 if(! opt)
734 return;
735
736 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 737 AliAnaCaloTrackCorrBaseClass::Print("");
243c2909 738 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
739 printf("lambda 0 squared > %2.1f\n", fM02Cut);
740 printf("pi0 : %2.2f<m<%2.2f \n", fMassPi0Min,fMassPi0Max);
741 printf("eta : %2.2f<m<%2.2f \n", fMassEtaMin,fMassEtaMax);
742 printf("conv: %2.2f<m<%2.2f \n", fMassConMin,fMassConMax);
743
992b14a7 744 printf(" \n") ;
745
746}
747
748//____________________________________________________________________________________________
749void AliAnaInsideClusterInvariantMass::RecalibrateCellAmplitude(Float_t & amp, const Int_t id)
750{
751 //Recaculate cell energy if recalibration factor
752
753 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
754 Int_t nModule = GetModuleNumberCellIndexes(id,fCalorimeter, icol, irow, iRCU);
755
756 if (GetCaloUtils()->IsRecalibrationOn()) {
757 if(fCalorimeter == "PHOS") {
758 amp *= GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
759 }
760 else {
761 amp *= GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
762 }
763 }
764}
765
c5693f62 766//________________________________________________________________________________________
767void AliAnaInsideClusterInvariantMass::SplitEnergy(const Int_t absId1, const Int_t absId2,
243c2909 768 AliVCluster* cluster,
769 AliVCaloCells* cells,
992b14a7 770 Float_t & e1, Float_t & e2 )
771{
772
773 // Split energy of cluster between the 2 local maxima.
243c2909 774 const Int_t ncells = cluster->GetNCells();
775 Int_t absIdList[ncells];
776 //printf("Split Local Max: 1) %d - 2) %d\n",absId1,absId2);
777 for(Int_t iDigit = 0; iDigit < ncells; iDigit++ ) {
778 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
779// printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit],cells->GetCellAmplitude(absIdList[iDigit]));
992b14a7 780 }
243c2909 781
782 // SubCluster 1
783
784 // Init counters and variables
785 Int_t ncells1 = 1;
786 Int_t absIdList1[ncells];
787 absIdList1[0] = absId1;
788 //printf("First local max : absId1 %d %d \n",absId1, absIdList1[0]);
789 for(Int_t iDigit1 = 1; iDigit1 < ncells; iDigit1++) absIdList1[iDigit1] = -1;
790
791 Float_t ecell1 = cells->GetCellAmplitude(absId1);
792 RecalibrateCellAmplitude(ecell1,absId1);
793 e1 = ecell1;
794
795 Bool_t added = kTRUE;
796 while (added)
797 {
798 added = kFALSE;
799 Int_t absId1New = absIdList1[ncells1-1];
800 //printf("\t absId %d added \n",absId1New);
801
802 Float_t e1New = cells->GetCellAmplitude(absId1New);
803 RecalibrateCellAmplitude(e1New,absId1New);
804
805 for(Int_t iDigit = 0; iDigit < ncells ; iDigit++)
806 {
807 Int_t absId = absIdList[iDigit] ;
808 if(absId!=absId1New && absId!=absId2 && absId>=0)
809 {
810 //printf("\t \t iDig %d, absId %d, absIdNew %d\n",iDigit,absId, absId1New);
811 if(AreNeighbours( absId1New,absId )){
812 //printf("\t neighbours\n");
813
814 Float_t en = cells->GetCellAmplitude(absId);
815 RecalibrateCellAmplitude(en,absId);
816 //printf("\t \t e1New %f, en %f \n",e1New,en);
817 if(e1New > en){
818 absIdList1[ncells1++] = absId;
819 absIdList [iDigit] = -1;
820 e1+=en;
821 added = kTRUE;
822 } // Decreasing energy with respect reference
823 } // Neighbours
824 } //Not local maxima or already removed
825 } // cell loop
826
827 }// while cells added to list of cells for cl1
828
829 // SubCluster 2
830
831 // Init counters and variables
832 Int_t ncells2 = 1;
833 Int_t absIdList2[ncells];
834 absIdList2[0] = absId2;
835 //printf("Second local max : absId2 %d %d \n",absId2, absIdList2[0]);
836 for(Int_t iDigit2 = 1; iDigit2 < ncells; iDigit2++) absIdList2[iDigit2] = -1;
837
838 Float_t ecell2 = cells->GetCellAmplitude(absId2);
839 RecalibrateCellAmplitude(ecell2,absId2);
840 e2 = ecell2;
841
842 added = kTRUE;
843 while (added)
844 {
845 added = kFALSE;
846 Int_t absId2New = absIdList2[ncells2-1];
847 //printf("\t absId %d added \n",absId2New);
848
849 Float_t e2New = cells->GetCellAmplitude(absId2New);
850 RecalibrateCellAmplitude(e2New,absId2New);
851
852 for(Int_t iDigit = 0; iDigit < ncells ; iDigit++)
853 {
854 Int_t absId = absIdList[iDigit] ;
855 if(absId!=absId2New && absId>=0)
856 {
857 // printf("\t \t iDig %d, absId %d, absIdNew %d\n",iDigit,absId, absId2New);
858 if(AreNeighbours( absId2New,absId )){
859 // printf("\t neighbours\n");
860
861 Float_t en = cells->GetCellAmplitude(absId);
862 RecalibrateCellAmplitude(en,absId);
863 //printf("\t \t e2New %f, en %f \n",e2New,en);
864 if(e2New > en){
865 absIdList2[ncells2++] = absId;
866 absIdList [iDigit] = -1;
867 e2+=en;
868 added = kTRUE;
869 } // Decreasing energy with respect reference
870 } // Neighbours
871 } //Not local maxima or already removed
872 } // cell loop
873
874 }// while cells added to list of cells for cl2
875
876 //printf("Cluster energy = %f, Ecell1 = %f, Ecell2 = %f, Enew1 = %f, Enew2 = %f, ncells %d, ncells1 %d, ncells2 %d \n",
877 // cluster->E(),ecell1,ecell2,e1,e2,ncells,ncells1,ncells2);
878 //if(ncells!=(ncells1+ncells2)) printf("\t Not all cells!\n");
879
992b14a7 880}
881
882