]> 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
Rename base classes from PartCorr to CaloTrackCorr, agreed new naming and directory...
[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(""),
56 fM02Cut(0),
57 fMinNCells(0)
58{
59 //default ctor
60
61 // Init array of histograms
62 for(Int_t i = 0; i < 7; i++){
63 fhMassNLocMax1[i] = 0;
64 fhMassNLocMax2[i] = 0;
65 fhMassNLocMaxN[i] = 0;
66 fhNLocMax[i] = 0;
67 fhNLocMaxM02Cut[i] = 0;
68 fhM02NLocMax1[i] = 0;
69 fhM02NLocMax2[i] = 0;
70 fhM02NLocMaxN[i] = 0;
71 fhNCellNLocMax1[i] = 0;
72 fhNCellNLocMax2[i] = 0;
73 fhNCellNLocMaxN[i] = 0;
74 fhM02Pi0[i] = 0;
75 fhM02Eta[i] = 0;
76 fhM02Con[i] = 0;
77 fhInvMassAllCells[i]=0;
78 }
79
80 InitParameters();
81
82}
83
84//_______________________________________________________________
85TObjString * AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
86{
87 //Save parameters used for analysis
88 TString parList ; //this will be list of parameters used for this analysis.
89 const Int_t buffersize = 255;
90 char onePar[buffersize] ;
91
92 snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
93 parList+=onePar ;
94
95 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
96 parList+=onePar ;
97 snprintf(onePar,buffersize,"fM02Cut =%f \n" ,fM02Cut) ;
98 parList+=onePar ;
99 snprintf(onePar,buffersize,"fMinNCells =%f \n",fMinNCells) ;
100 parList+=onePar ;
101
102 return new TObjString(parList) ;
103
104}
105
106//____________________________________________________________________________________________________
107Bool_t AliAnaInsideClusterInvariantMass::AreNeighbours( const Int_t absId1, const Int_t absId2 ) const
108{
109 // Tells if (true) or not (false) two digits are neighbours
110 // A neighbour is defined as being two digits which share a corner
111
112 Bool_t areNeighbours = kFALSE ;
113 Int_t nSupMod =0, nModule =0, nIphi =0, nIeta =0;
114 Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0;
115 Int_t relid1[2], relid2[2] ;
116 Int_t rowdiff=0, coldiff=0;
117
118 areNeighbours = kFALSE ;
119
120 GetEMCALGeometry()->GetCellIndex(absId1, nSupMod,nModule,nIphi,nIeta);
121 GetEMCALGeometry()->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, relid1[0],relid1[1]);
122
123 GetEMCALGeometry()->GetCellIndex(absId2, nSupMod1,nModule1,nIphi1,nIeta1);
124 GetEMCALGeometry()->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, relid2[0],relid2[1]);
125
126 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
127 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
128 if(nSupMod1!=nSupMod){
129 if(nSupMod1%2) relid1[1]+=AliEMCALGeoParams::fgkEMCALCols;
130 else relid2[1]+=AliEMCALGeoParams::fgkEMCALCols;
131 }
132
133 rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;
134 coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;
135
136 if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
137 areNeighbours = kTRUE ;
138
139 return areNeighbours;
140}
141
142//_____________________________________________________________________________________
143TLorentzVector AliAnaInsideClusterInvariantMass::GetCellMomentum(const Int_t absId,
144 AliVCaloCells * cells)
145{
146
147 // Cell momentum calculation for invariant mass
148
149 Double_t cellpos[] = {0, 0, 0};
150 GetEMCALGeometry()->GetGlobal(absId, cellpos);
151
152 if(GetVertex(0)){//calculate direction from vertex
153 cellpos[0]-=GetVertex(0)[0];
154 cellpos[1]-=GetVertex(0)[1];
155 cellpos[2]-=GetVertex(0)[2];
156 }
157
158 Double_t r = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]+cellpos[2]*cellpos[2] ) ;
159
160 Float_t en = cells->GetCellAmplitude(absId);
161 RecalibrateCellAmplitude(en,absId);
162 TLorentzVector cellMom ;
163 cellMom.SetPxPyPzE( en*cellpos[0]/r, en*cellpos[1]/r, en*cellpos[2]/r, en) ;
164
165 return cellMom;
166
167}
168
169//________________________________________________________________
170TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
171{
172 // Create histograms to be saved in output file and
173 // store them in outputContainer
174 TList * outputContainer = new TList() ;
175 outputContainer->SetName("InsideClusterHistos") ;
176
745913ae 177 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins(); Float_t ptmax = GetHistogramRanges()->GetHistoPtMax(); Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
178 Int_t ssbins = GetHistogramRanges()->GetHistoShowerShapeBins(); Float_t ssmax = GetHistogramRanges()->GetHistoShowerShapeMax(); Float_t ssmin = GetHistogramRanges()->GetHistoShowerShapeMin();
179 Int_t mbins = GetHistogramRanges()->GetHistoMassBins(); Float_t mmax = GetHistogramRanges()->GetHistoMassMax(); Float_t mmin = GetHistogramRanges()->GetHistoMassMin();
180 Int_t ncbins = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t ncmax = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t ncmin = GetHistogramRanges()->GetHistoNClusterCellMin();
992b14a7 181
182 TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
183 TString pname[] ={"","Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
184
185 Int_t n = 1;
186
187 if(IsDataMC()) n = 7;
188
189 for(Int_t i = 0; i < n; i++){
190
191 fhInvMassAllCells[i] = new TH2F(Form("hInvMassAllCells%s",pname[i].Data()),
192 Form("Invariant mass of all cells %s",ptype[i].Data()),
193 nptbins,ptmin,ptmax,mbins,mmin,mmax);
194 fhInvMassAllCells[i]->SetYTitle("M (MeV/c^2)");
195 fhInvMassAllCells[i]->SetXTitle("E (GeV)");
196 outputContainer->Add(fhInvMassAllCells[i]) ;
197
198
199 fhMassNLocMax1[i] = new TH2F(Form("hMassNLocMax1%s",pname[i].Data()),
200 Form("Invariant mass of 2 highest energy cells %s",ptype[i].Data()),
201 nptbins,ptmin,ptmax,mbins,mmin,mmax);
202 fhMassNLocMax1[i]->SetYTitle("M (MeV/c^2)");
203 fhMassNLocMax1[i]->SetXTitle("E (GeV)");
204 outputContainer->Add(fhMassNLocMax1[i]) ;
205
206 fhMassNLocMax2[i] = new TH2F(Form("hMassNLocMax2%s",pname[i].Data()),
207 Form("Invariant mass of 2 local maxima cells %s",ptype[i].Data()),
208 nptbins,ptmin,ptmax,mbins,mmin,mmax);
209 fhMassNLocMax2[i]->SetYTitle("M (MeV/c^2)");
210 fhMassNLocMax2[i]->SetXTitle("E (GeV)");
211 outputContainer->Add(fhMassNLocMax2[i]) ;
212
213 fhMassNLocMaxN[i] = new TH2F(Form("hMassNLocMaxN%s",pname[i].Data()),
214 Form("Invariant mass of 2 local maxima cells %s",ptype[i].Data()),
215 nptbins,ptmin,ptmax,mbins,mmin,mmax);
216 fhMassNLocMaxN[i]->SetYTitle("M (MeV/c^2)");
217 fhMassNLocMax2[i]->SetXTitle("E (GeV)");
218 outputContainer->Add(fhMassNLocMaxN[i]) ;
219
220
221 fhNLocMax[i] = new TH2F(Form("hNLocMax%s",pname[i].Data()),
222 Form("Number of local maxima in cluster %s",ptype[i].Data()),
223 nptbins,ptmin,ptmax,5,0,5);
224 fhNLocMax[i] ->SetYTitle("N maxima");
225 fhNLocMax[i] ->SetXTitle("E (GeV)");
226 outputContainer->Add(fhNLocMax[i]) ;
227
228 fhNLocMaxM02Cut[i] = new TH2F(Form("hNLocMaxM02Cut%s",pname[i].Data()),
229 Form("Number of local maxima in cluster %s for M02 > %2.2f",ptype[i].Data(),fM02Cut),
230 nptbins,ptmin,ptmax,5,0,5);
231 fhNLocMaxM02Cut[i]->SetYTitle("N maxima");
232 fhNLocMaxM02Cut[i]->SetXTitle("E (GeV)");
233 outputContainer->Add(fhNLocMaxM02Cut[i]) ;
234
235
236 fhM02NLocMax1[i] = new TH2F(Form("hM02NLocMax1%s",pname[i].Data()),
237 Form("#lambda_{0}^{2} vs E for N max = 1 %s",ptype[i].Data()),
238 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
239 fhM02NLocMax1[i] ->SetYTitle("#lambda_{0}^{2}");
240 fhM02NLocMax1[i] ->SetXTitle("E (GeV)");
241 outputContainer->Add(fhM02NLocMax1[i]) ;
242
243 fhM02NLocMax2[i] = new TH2F(Form("hM02NLocMax2%s",pname[i].Data()),
244 Form("#lambda_{0}^{2} vs E for N max = 2 %s",ptype[i].Data()),
245 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
246 fhM02NLocMax2[i] ->SetYTitle("#lambda_{0}^{2}");
247 fhM02NLocMax2[i] ->SetXTitle("E (GeV)");
248 outputContainer->Add(fhM02NLocMax2[i]) ;
249
250
251 fhM02NLocMaxN[i] = new TH2F(Form("hM02NLocMaxN%s",pname[i].Data()),
252 Form("#lambda_{0}^{2} vs E for N max > 2 %s",ptype[i].Data()),
253 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
254 fhM02NLocMaxN[i] ->SetYTitle("#lambda_{0}^{2}");
255 fhM02NLocMaxN[i] ->SetXTitle("E (GeV)");
256 outputContainer->Add(fhM02NLocMaxN[i]) ;
257
258
259 fhNCellNLocMax1[i] = new TH2F(Form("hNCellNLocMax1%s",pname[i].Data()),
260 Form("#lambda_{0}^{2} vs E for N max = 1 %s",ptype[i].Data()),
261 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
262 fhNCellNLocMax1[i] ->SetYTitle("N cells");
263 fhNCellNLocMax1[i] ->SetXTitle("E (GeV)");
264 outputContainer->Add(fhNCellNLocMax1[i]) ;
265
266 fhNCellNLocMax2[i] = new TH2F(Form("hNCellNLocMax2%s",pname[i].Data()),
267 Form("#lambda_{0}^{2} vs E for N max = 2 %s",ptype[i].Data()),
268 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
269 fhNCellNLocMax2[i] ->SetYTitle("N cells");
270 fhNCellNLocMax2[i] ->SetXTitle("E (GeV)");
271 outputContainer->Add(fhNCellNLocMax2[i]) ;
272
273
274 fhNCellNLocMaxN[i] = new TH2F(Form("hNCellNLocMaxN%s",pname[i].Data()),
275 Form("#lambda_{0}^{2} vs E for N max > 2 %s",ptype[i].Data()),
276 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
277 fhNCellNLocMaxN[i] ->SetYTitle("N cells");
278 fhNCellNLocMaxN[i] ->SetXTitle("E (GeV)");
279 outputContainer->Add(fhNCellNLocMaxN[i]) ;
280
281
282 fhM02Pi0[i] = new TH2F(Form("hM02Pi0%s",pname[i].Data()),
283 Form("#lambda_{0}^{2} vs E ffor mass [0.1-0.2] MeV/c2 %s",ptype[i].Data()),
284 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
285 fhM02Pi0[i] ->SetYTitle("#lambda_{0}^{2}");
286 fhM02Pi0[i] ->SetXTitle("E (GeV)");
287 outputContainer->Add(fhM02Pi0[i]) ;
288
289 fhM02Eta[i] = new TH2F(Form("hM02Eta%s",pname[i].Data()),
290 Form("#lambda_{0}^{2} vs E for mass [0.4-0.7] MeV/c2, %s",ptype[i].Data()),
291 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
292 fhM02Eta[i] ->SetYTitle("#lambda_{0}^{2}");
293 fhM02Eta[i] ->SetXTitle("E (GeV)");
294 outputContainer->Add(fhM02Eta[i]) ;
295
296
297 fhM02Con[i] = new TH2F(Form("hM02Con%s",pname[i].Data()),
298 Form("#lambda_{0}^{2} vs E for mass < 0.5 MeV/c2, %s",ptype[i].Data()),
299 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
300 fhM02Con[i] ->SetYTitle("#lambda_{0}^{2}");
301 fhM02Con[i] ->SetXTitle("E (GeV)");
302 outputContainer->Add(fhM02Con[i]) ;
303
304 }
305
306 return outputContainer ;
307
308}
309
310//________________________________________________________________________________________________________
311Int_t AliAnaInsideClusterInvariantMass::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
312 Int_t *absIdList, Float_t *maxEList)
313{
314 // Find local maxima in cluster
315
316 Float_t locMaxCut = 0; // not used for the moment
317
318 Int_t iDigitN = 0 ;
319 Int_t iDigit = 0 ;
320 Int_t absId1 = -1 ;
321 Int_t absId2 = -1 ;
322 const Int_t nCells = cluster->GetNCells();
323 //printf("cluster :");
324 for(iDigit = 0; iDigit < nCells ; iDigit++){
325 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
326
327 //Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
328 //RecalibrateCellAmplitude(en,absIdList[iDigit]);
329 //Int_t icol = -1, irow = -1, iRCU = -1;
330 //Int_t sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList[iDigit], fCalorimeter, icol, irow, iRCU) ;
331
332 //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
333 }
334
335 for(iDigit = 0 ; iDigit < nCells; iDigit++) {
336 if(absIdList[iDigit]>=0) {
337
338 absId1 = absIdList[iDigit] ;
339 Float_t en1 = cells->GetCellAmplitude(absId1);
340 RecalibrateCellAmplitude(en1,absId1);
341
342 for(iDigitN = 0; iDigitN < nCells; iDigitN++) {
343 absId2 = absIdList[iDigitN] ;
344
345 Float_t en2 = cells->GetCellAmplitude(absId2);
346 RecalibrateCellAmplitude(en2,absId2);
347
348 if ( AreNeighbours(absId1, absId2) ) {
349
350 if (en1 > en2 ) {
351 absIdList[iDigitN] = -1 ;
352 // but may be digit too is not local max ?
353 if(en1 < en2 + locMaxCut)
354 absIdList[iDigit] = -1 ;
355 }
356 else {
357 absIdList[iDigit] = -1 ;
358 // but may be digitN too is not local max ?
359 if(en1 > en2 - locMaxCut)
360 absIdList[iDigitN] = -1 ;
361 }
362 } // if Areneighbours
363 } // while digitN
364 } // slot not empty
365 } // while digit
366
367 iDigitN = 0 ;
368 for(iDigit = 0; iDigit < nCells; iDigit++) {
369 if(absIdList[iDigit]>=0 ){
370 absIdList[iDigitN] = absIdList[iDigit] ;
371 Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
372 RecalibrateCellAmplitude(en,absIdList[iDigit]);
373 if(en < 0.1) continue; // Maxima only with seed energy at least
374 maxEList[iDigitN] = en ;
375 //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
376 iDigitN++ ;
377 }
378 }
379
380 return iDigitN ;
381
382}
383
384
385//___________________________________________
386void AliAnaInsideClusterInvariantMass::Init()
387{
388 //Init
389 //Do some checks
390 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
391 printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
392 abort();
393 }
394 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
395 printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
396 abort();
397 }
398
399 if( GetReader()->GetDataType() == AliCaloTrackReader::kMC ){
400 printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
401 abort();
402
403 }
404
405}
406
407//_____________________________________________________
408void AliAnaInsideClusterInvariantMass::InitParameters()
409{
410 //Initialize the parameters of the analysis.
411 AddToHistogramsName("AnaPi0InsideClusterInvariantMass_");
412
413 fCalorimeter = "EMCAL" ;
414 fM02Cut = 0.26 ;
415 fMinNCells = 4 ;
416}
417
418
419//__________________________________________________________________
420void AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
421{
422 //Search for pi0 in fCalorimeter with shower shape analysis
423
424 TObjArray * pl = 0x0;
425 AliVCaloCells* cells = 0x0;
426
427 //Select the Calorimeter of the photon
428 if(fCalorimeter == "PHOS"){
429 pl = GetPHOSClusters();
430 cells = GetPHOSCells();
431 }
432 else if (fCalorimeter == "EMCAL"){
433 pl = GetEMCALClusters();
434 cells = GetEMCALCells();
435 }
436
437 if(!pl || !cells) {
438 Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
439 return;
440 }
441
442 if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
443
444 for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++){
445 AliVCluster * cluster = (AliVCluster*) (pl->At(icluster));
446
447 // Study clusters with large shape parameter
448 Float_t en = cluster->E();
449 Float_t l0 = cluster->GetM02();
450 Int_t nc = cluster->GetNCells();
451
452 //If too small or big E or low number of cells, skip it
453 if( ( en < GetMinEnergy() || en > GetMaxEnergy() ) && nc < fMinNCells) continue ;
454
455 Int_t absId1 = -1; Int_t absId2 = -1;
456 Int_t *absIdList = new Int_t [nc];
457 Float_t *maxEList = new Float_t[nc];
458 Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
459
460 if (nMax <= 0) {
461 printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found!");
c5693f62 462 delete [] absIdList ;
463 delete [] maxEList ;
992b14a7 464 return;
465 }
466
467 fhNLocMax[0]->Fill(en,nMax);
468
469 if ( nMax == 1 ) { fhM02NLocMax1[0]->Fill(en,l0) ; fhNCellNLocMax1[0]->Fill(en,nc) ; }
470 else if( nMax == 2 ) { fhM02NLocMax2[0]->Fill(en,l0) ; fhNCellNLocMax2[0]->Fill(en,nc) ; }
471 else if( nMax >= 3 ) { fhM02NLocMaxN[0]->Fill(en,l0) ; fhNCellNLocMaxN[0]->Fill(en,nc) ; }
472 else printf("N max smaller than 1 -> %d \n",nMax);
473
474 // Play with the MC stack if available
475 // Check origin of the candidates
476 Int_t mcindex = -1;
477 if(IsDataMC()){
478
479 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader(), 0);
480
c5693f62 481 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ) mcindex = kmcPi0;
482 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mcindex = kmcEta;
992b14a7 483 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
c5693f62 484 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPhoton;
992b14a7 485 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
c5693f62 486 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
487 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)) mcindex = kmcElectron;
488 else mcindex = kmcHadron;
992b14a7 489
490/* 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,
491 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton),
492 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt),
493 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation),
494 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR),
495 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay),
496 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay),
497 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay),
498
499 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
500 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
501 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
502 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
503 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
504 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon),
505 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
506 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton),
507 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
508 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon),
509 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton),
510 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
511 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCBadLabel)
512);
513*/
514
515 fhNLocMax[mcindex]->Fill(en,nMax);
516
517 if (nMax == 1 ) { fhM02NLocMax1[mcindex]->Fill(en,l0) ; fhNCellNLocMax1[mcindex]->Fill(en,nc) ; }
518 else if(nMax == 2 ) { fhM02NLocMax2[mcindex]->Fill(en,l0) ; fhNCellNLocMax2[mcindex]->Fill(en,nc) ; }
519 else if(nMax >= 3 ) { fhM02NLocMaxN[mcindex]->Fill(en,l0) ; fhNCellNLocMaxN[mcindex]->Fill(en,nc) ; }
520
521 }
522
37f53c3e 523 if( l0 < fM02Cut) {
524 delete [] absIdList ;
525 delete [] maxEList ;
526 continue;
527 }
992b14a7 528
529 // Get the 2 max indeces and do inv mass
530
531 absId1 = absIdList[0];
532 TLorentzVector cellMomi = GetCellMomentum(absId1, cells);
533
534 if ( nMax == 2 ) absId2 = absIdList[1];
535 else if( nMax == 1 ){
536 //Find second highest energy cell
04c5e581 537 Float_t enmax = 0 ;
992b14a7 538 for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++){
539 Int_t absId = cluster->GetCellsAbsId()[iDigit];
540 if( absId == absId1 ) continue ;
541 Float_t endig = cells->GetCellAmplitude(absId);
542 RecalibrateCellAmplitude(endig,absId);
543 if(endig > enmax) {
544 enmax = endig ;
545 absId2 = absId ;
546 }
547
548 //Get mass of all cell in cluster combinations
549
550
551 Float_t enj = cells->GetCellAmplitude(absId);
552 RecalibrateCellAmplitude(enj,absId);
553
37f53c3e 554 if(enj < 0.3) continue;
992b14a7 555
556 TLorentzVector cellMomj = GetCellMomentum(absId, cells);
557
558 fhInvMassAllCells[0]->Fill(en,(cellMomj+cellMomi).M());
559
560 if(IsDataMC()) fhInvMassAllCells[mcindex]->Fill(en,(cellMomj+cellMomi).M());
561
562 }// cell loop
563
564 }
565 else { // loop on maxima, find 2 highest
566
567 Int_t enmax = 0 ;
568 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
569 Float_t endig = maxEList[iDigit];
570 if(endig > enmax) {
571 endig = enmax ;
572 absId2 = absIdList[iDigit];
573 }
574 }// maxima loop
575
576 // First max is not highest, check if there is other higher
577 if(maxEList[0] < enmax){
578
579 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
580 if(absId2 == absIdList[iDigit]) continue;
581 Float_t endig = maxEList[iDigit];
582 if(endig > enmax) {
583 endig = enmax ;
584 absId1 = absIdList[iDigit];
585 }
586 }// maxima loop
587
588 }
589
590 }
591
592 fhNLocMaxM02Cut[0]->Fill(en,nMax);
593
594 TLorentzVector cellMom1 = GetCellMomentum(absId1, cells);
595 TLorentzVector cellMom2 = GetCellMomentum(absId2, cells);
596
597 Float_t mass = (cellMom1+cellMom2).M();
598
599 if (nMax==1) fhMassNLocMax1[0]->Fill(en,mass);
600 else if(nMax==2) fhMassNLocMax2[0]->Fill(en,mass);
601 else if(nMax >2) fhMassNLocMaxN[0]->Fill(en,mass);
602
603 if (mass < 0.1) fhM02Con[0]->Fill(en,l0);
604 else if(mass < 0.2) fhM02Pi0[0]->Fill(en,l0);
605 else if(mass < 0.6 && mass > 0.4) fhM02Eta[0]->Fill(en,l0);
606
607 if(IsDataMC()){
608
609 fhNLocMaxM02Cut[mcindex]->Fill(en,nMax);
610
611 if (nMax==1) fhMassNLocMax1[mcindex]->Fill(en,mass);
612 else if(nMax==2) fhMassNLocMax2[mcindex]->Fill(en,mass);
613 else if(nMax >2) fhMassNLocMaxN[mcindex]->Fill(en,mass);
614
615 if (mass < 0.1) fhM02Con[mcindex]->Fill(en,l0);
616 else if(mass < 0.2) fhM02Pi0[mcindex]->Fill(en,l0);
617 else if(mass < 0.6 && mass > 0.4) fhM02Eta[mcindex]->Fill(en,l0);
618
619
620 }//Work with MC truth first
621
c5693f62 622 delete [] absIdList ;
623 delete [] maxEList ;
992b14a7 624
625 }//loop
626
627 if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");
628
629}
630
631//______________________________________________________________________
632void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
633{
634 //Print some relevant parameters set for the analysis
635 if(! opt)
636 return;
637
638 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
745913ae 639 AliAnaCaloTrackCorrBaseClass::Print("");
992b14a7 640 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
641 printf("lambda 0 sqared > %2.1f\n", fM02Cut);
642 printf(" \n") ;
643
644}
645
646//____________________________________________________________________________________________
647void AliAnaInsideClusterInvariantMass::RecalibrateCellAmplitude(Float_t & amp, const Int_t id)
648{
649 //Recaculate cell energy if recalibration factor
650
651 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
652 Int_t nModule = GetModuleNumberCellIndexes(id,fCalorimeter, icol, irow, iRCU);
653
654 if (GetCaloUtils()->IsRecalibrationOn()) {
655 if(fCalorimeter == "PHOS") {
656 amp *= GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
657 }
658 else {
659 amp *= GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
660 }
661 }
662}
663
c5693f62 664//________________________________________________________________________________________
665void AliAnaInsideClusterInvariantMass::SplitEnergy(const Int_t absId1, const Int_t absId2,
666 const AliVCaloCells* cells,
992b14a7 667 Float_t & e1, Float_t & e2 )
668{
669
670 // Split energy of cluster between the 2 local maxima.
671
672 Int_t icol1 = -1, irow1 = -1, iRCU1 = -1;
673 Int_t sm1 = GetCaloUtils()->GetModuleNumberCellIndexes(absId1, fCalorimeter, icol1, irow1, iRCU1) ;
674 Int_t icol2 = -1, irow2 = -1, iRCU2 = -1;
675 Int_t sm2 = GetCaloUtils()->GetModuleNumberCellIndexes(absId2, fCalorimeter, icol2, irow2, iRCU2) ;
676
677 if(sm1!=sm2) {
678 if(sm1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
679 else icol2+=AliEMCALGeoParams::fgkEMCALCols;
680 }
681
04c5e581 682 // just to avoid compilation warning
683 Int_t nTotCells = cells->GetNumberOfCells();
684 if(GetDebug() > 2)printf("N cells %d, e1 %f, e2 %f\n", nTotCells,e1, e2);
992b14a7 685/// continue here
686
687}
688
689