1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //_________________________________________________________________________
18 // Split clusters with some criteria and calculate invariant mass
19 // to identify them as pi0 or conversion
22 //-- Author: Gustavo Conesa (LPSC-Grenoble)
23 //_________________________________________________________________________
25 //////////////////////////////////////////////////////////////////////////////
28 // --- ROOT system ---
30 #include <TClonesArray.h>
31 #include <TObjString.h>
34 // --- Analysis system ---
35 #include "AliAnaInsideClusterInvariantMass.h"
36 #include "AliCaloTrackReader.h"
37 #include "AliMCAnalysisUtils.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"
47 //#include "AliPHOSGeoUtils.h"
48 #include "AliEMCALGeometry.h"
50 ClassImp(AliAnaInsideClusterInvariantMass)
52 //__________________________________________________________________
53 AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
54 AliAnaPartCorrBaseClass(),
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;
67 fhNLocMaxM02Cut[i] = 0;
71 fhNCellNLocMax1[i] = 0;
72 fhNCellNLocMax2[i] = 0;
73 fhNCellNLocMaxN[i] = 0;
77 fhInvMassAllCells[i]=0;
84 //_______________________________________________________________
85 TObjString * AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
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] ;
92 snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
95 snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
97 snprintf(onePar,buffersize,"fM02Cut =%f \n" ,fM02Cut) ;
99 snprintf(onePar,buffersize,"fMinNCells =%f \n",fMinNCells) ;
102 return new TObjString(parList) ;
106 //____________________________________________________________________________________________________
107 Bool_t AliAnaInsideClusterInvariantMass::AreNeighbours( const Int_t absId1, const Int_t absId2 ) const
109 // Tells if (true) or not (false) two digits are neighbours
110 // A neighbour is defined as being two digits which share a corner
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;
118 areNeighbours = kFALSE ;
120 GetEMCALGeometry()->GetCellIndex(absId1, nSupMod,nModule,nIphi,nIeta);
121 GetEMCALGeometry()->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, relid1[0],relid1[1]);
123 GetEMCALGeometry()->GetCellIndex(absId2, nSupMod1,nModule1,nIphi1,nIeta1);
124 GetEMCALGeometry()->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, relid2[0],relid2[1]);
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;
133 rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;
134 coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;
136 if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
137 areNeighbours = kTRUE ;
139 return areNeighbours;
142 //_____________________________________________________________________________________
143 TLorentzVector AliAnaInsideClusterInvariantMass::GetCellMomentum(const Int_t absId,
144 AliVCaloCells * cells)
147 // Cell momentum calculation for invariant mass
149 Double_t cellpos[] = {0, 0, 0};
150 GetEMCALGeometry()->GetGlobal(absId, cellpos);
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];
158 Double_t r = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]+cellpos[2]*cellpos[2] ) ;
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) ;
169 //________________________________________________________________
170 TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
172 // Create histograms to be saved in output file and
173 // store them in outputContainer
174 TList * outputContainer = new TList() ;
175 outputContainer->SetName("InsideClusterHistos") ;
177 Int_t nptbins = GetHistoPtBins(); Float_t ptmax = GetHistoPtMax(); Float_t ptmin = GetHistoPtMin();
178 Int_t ssbins = GetHistoShowerShapeBins(); Float_t ssmax = GetHistoShowerShapeMax(); Float_t ssmin = GetHistoShowerShapeMin();
179 Int_t mbins = GetHistoMassBins(); Float_t mmax = GetHistoMassMax(); Float_t mmin = GetHistoMassMin();
180 Int_t ncbins = GetHistoNClusterCellBins(); Int_t ncmax = GetHistoNClusterCellMax(); Int_t ncmin = GetHistoNClusterCellMin();
182 TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
183 TString pname[] ={"","Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
187 if(IsDataMC()) n = 7;
189 for(Int_t i = 0; i < n; i++){
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
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]) ;
306 return outputContainer ;
310 //________________________________________________________________________________________________________
311 Int_t AliAnaInsideClusterInvariantMass::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
312 Int_t *absIdList, Float_t *maxEList)
314 // Find local maxima in cluster
316 Float_t locMaxCut = 0; // not used for the moment
322 const Int_t nCells = cluster->GetNCells();
323 //printf("cluster :");
324 for(iDigit = 0; iDigit < nCells ; iDigit++){
325 absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
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) ;
332 //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
335 for(iDigit = 0 ; iDigit < nCells; iDigit++) {
336 if(absIdList[iDigit]>=0) {
338 absId1 = absIdList[iDigit] ;
339 Float_t en1 = cells->GetCellAmplitude(absId1);
340 RecalibrateCellAmplitude(en1,absId1);
342 for(iDigitN = 0; iDigitN < nCells; iDigitN++) {
343 absId2 = absIdList[iDigitN] ;
345 Float_t en2 = cells->GetCellAmplitude(absId2);
346 RecalibrateCellAmplitude(en2,absId2);
348 if ( AreNeighbours(absId1, absId2) ) {
351 absIdList[iDigitN] = -1 ;
352 // but may be digit too is not local max ?
353 if(en1 < en2 + locMaxCut)
354 absIdList[iDigit] = -1 ;
357 absIdList[iDigit] = -1 ;
358 // but may be digitN too is not local max ?
359 if(en1 > en2 - locMaxCut)
360 absIdList[iDigitN] = -1 ;
362 } // if Areneighbours
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);
385 //___________________________________________
386 void AliAnaInsideClusterInvariantMass::Init()
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");
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");
399 if( GetReader()->GetDataType() == AliCaloTrackReader::kMC ){
400 printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
407 //_____________________________________________________
408 void AliAnaInsideClusterInvariantMass::InitParameters()
410 //Initialize the parameters of the analysis.
411 AddToHistogramsName("AnaPi0InsideClusterInvariantMass_");
413 fCalorimeter = "EMCAL" ;
419 //__________________________________________________________________
420 void AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
422 //Search for pi0 in fCalorimeter with shower shape analysis
424 TObjArray * pl = 0x0;
425 AliVCaloCells* cells = 0x0;
427 //Select the Calorimeter of the photon
428 if(fCalorimeter == "PHOS"){
429 pl = GetPHOSClusters();
430 cells = GetPHOSCells();
432 else if (fCalorimeter == "EMCAL"){
433 pl = GetEMCALClusters();
434 cells = GetEMCALCells();
438 Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
442 if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
444 for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++){
445 AliVCluster * cluster = (AliVCluster*) (pl->At(icluster));
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();
452 //If too small or big E or low number of cells, skip it
453 if( ( en < GetMinEnergy() || en > GetMaxEnergy() ) && nc < fMinNCells) continue ;
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) ;
461 printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found!");
462 delete [] absIdList ;
467 fhNLocMax[0]->Fill(en,nMax);
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);
474 // Play with the MC stack if available
475 // Check origin of the candidates
479 Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader(), 0);
481 if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ) mcindex = kmcPi0;
482 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ) mcindex = kmcEta;
483 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
484 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPhoton;
485 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
486 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
487 else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)) mcindex = kmcElectron;
488 else mcindex = kmcHadron;
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),
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)
515 fhNLocMax[mcindex]->Fill(en,nMax);
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) ; }
524 delete [] absIdList ;
529 // Get the 2 max indeces and do inv mass
531 absId1 = absIdList[0];
532 TLorentzVector cellMomi = GetCellMomentum(absId1, cells);
534 if ( nMax == 2 ) absId2 = absIdList[1];
535 else if( nMax == 1 ){
536 //Find second highest energy cell
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);
548 //Get mass of all cell in cluster combinations
551 Float_t enj = cells->GetCellAmplitude(absId);
552 RecalibrateCellAmplitude(enj,absId);
554 if(enj < 0.3) continue;
556 TLorentzVector cellMomj = GetCellMomentum(absId, cells);
558 fhInvMassAllCells[0]->Fill(en,(cellMomj+cellMomi).M());
560 if(IsDataMC()) fhInvMassAllCells[mcindex]->Fill(en,(cellMomj+cellMomi).M());
565 else { // loop on maxima, find 2 highest
568 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
569 Float_t endig = maxEList[iDigit];
572 absId2 = absIdList[iDigit];
576 // First max is not highest, check if there is other higher
577 if(maxEList[0] < enmax){
579 for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
580 if(absId2 == absIdList[iDigit]) continue;
581 Float_t endig = maxEList[iDigit];
584 absId1 = absIdList[iDigit];
592 fhNLocMaxM02Cut[0]->Fill(en,nMax);
594 TLorentzVector cellMom1 = GetCellMomentum(absId1, cells);
595 TLorentzVector cellMom2 = GetCellMomentum(absId2, cells);
597 Float_t mass = (cellMom1+cellMom2).M();
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);
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);
609 fhNLocMaxM02Cut[mcindex]->Fill(en,nMax);
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);
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);
620 }//Work with MC truth first
622 delete [] absIdList ;
627 if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");
631 //______________________________________________________________________
632 void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
634 //Print some relevant parameters set for the analysis
638 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
639 AliAnaPartCorrBaseClass::Print("");
640 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
641 printf("lambda 0 sqared > %2.1f\n", fM02Cut);
646 //____________________________________________________________________________________________
647 void AliAnaInsideClusterInvariantMass::RecalibrateCellAmplitude(Float_t & amp, const Int_t id)
649 //Recaculate cell energy if recalibration factor
651 Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
652 Int_t nModule = GetModuleNumberCellIndexes(id,fCalorimeter, icol, irow, iRCU);
654 if (GetCaloUtils()->IsRecalibrationOn()) {
655 if(fCalorimeter == "PHOS") {
656 amp *= GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
659 amp *= GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
664 //________________________________________________________________________________________
665 void AliAnaInsideClusterInvariantMass::SplitEnergy(const Int_t absId1, const Int_t absId2,
666 const AliVCaloCells* cells,
667 Float_t & e1, Float_t & e2 )
670 // Split energy of cluster between the 2 local maxima.
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) ;
678 if(sm1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
679 else icol2+=AliEMCALGeoParams::fgkEMCALCols;
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);