PWGGA/CaloTrackCorrelations/AliAnaInsideClusterInvariantMass.cxx

   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>
  33 #include <TCanvas.h>
  34 #include <TStyle.h>
  35 #include <TPad.h>
  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
  49 // --- Detectors ---
  50 //#include "AliPHOSGeoUtils.h"
  51 #include "AliEMCALGeometry.h"
  52
  53 ClassImp(AliAnaInsideClusterInvariantMass)
  54
  55 //__________________________________________________________________
  56 AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
  57   AliAnaCaloTrackCorrBaseClass(),
  58   fCalorimeter(""),
  59   fM02Cut(0),       fMinNCells(0),
  60   fMassEtaMin(0),   fMassEtaMax(0),
  61   fMassPi0Min(0),   fMassPi0Max(0),
  62   fMassConMin(0),   fMassConMax(0)
  63 {
  64   //default ctor
  65
  66   // Init array of histograms
  67   for(Int_t i = 0; i < 7; i++)
  68   {
  69     for(Int_t j = 0; j < 2; j++)
  70     {
  71
  72       fhMassNLocMax1[i][j]  = 0;
  73       fhMassNLocMax2[i][j]  = 0;
  74       fhMassNLocMaxN[i][j]  = 0;
  75       fhNLocMax[i][j]       = 0;
  76       fhNLocMaxEMax[i][j]   = 0;
  77       fhNLocMaxEFrac[i][j]  = 0;
  78       fhNLocMaxM02Cut[i][j] = 0;
  79       fhM02NLocMax1[i][j]   = 0;
  80       fhM02NLocMax2[i][j]   = 0;
  81       fhM02NLocMaxN[i][j]   = 0;
  82       fhNCellNLocMax1[i][j] = 0;
  83       fhNCellNLocMax2[i][j] = 0;
  84       fhNCellNLocMaxN[i][j] = 0;
  85       fhM02Pi0LocMax1[i][j] = 0;
  86       fhM02EtaLocMax1[i][j] = 0;
  87       fhM02ConLocMax1[i][j] = 0;
  88       fhM02Pi0LocMax2[i][j] = 0;
  89       fhM02EtaLocMax2[i][j] = 0;
  90       fhM02ConLocMax2[i][j] = 0;
  91       fhM02Pi0LocMaxN[i][j] = 0;
  92       fhM02EtaLocMaxN[i][j] = 0;
  93       fhM02ConLocMaxN[i][j] = 0;
  94
  95     }
  96
  97     fhTrackMatchedDEtaLocMax1[i] = 0;
  98     fhTrackMatchedDPhiLocMax1[i] = 0;
  99     fhTrackMatchedDEtaLocMax2[i] = 0;
 100     fhTrackMatchedDPhiLocMax2[i] = 0;
 101     fhTrackMatchedDEtaLocMaxN[i] = 0;
 102     fhTrackMatchedDPhiLocMaxN[i] = 0;
 103
 104   }
 105
 106   for(Int_t i = 0; i < 2; i++)
 107   {
 108     fhAnglePairLocMax1    [i] = 0;
 109     fhAnglePairLocMax2    [i] = 0;
 110     fhAnglePairLocMaxN    [i] = 0;
 111     fhAnglePairMassLocMax1[i] = 0;
 112     fhAnglePairMassLocMax2[i] = 0;
 113     fhAnglePairMassLocMaxN[i] = 0;
 114   }
 115
 116   InitParameters();
 117
 118 }
 119
 120 //_______________________________________________________________
 121 TObjString *  AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
 122 {
 123         //Save parameters used for analysis
 124   TString parList ; //this will be list of parameters used for this analysis.
 125   const Int_t buffersize = 255;
 126   char onePar[buffersize] ;
 127
 128   snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
 129   parList+=onePar ;
 130
 131   snprintf(onePar,buffersize,"Calorimeter: %s\n",        fCalorimeter.Data()) ;
 132   parList+=onePar ;
 133   snprintf(onePar,buffersize,"fLocMaxCutE =%2.2f \n",    GetCaloUtils()->GetLocalMaximaCutE()) ;
 134   parList+=onePar ;
 135   snprintf(onePar,buffersize,"fLocMaxCutEDiff =%2.2f \n",GetCaloUtils()->GetLocalMaximaCutEDiff()) ;
 136   parList+=onePar ;
 137   snprintf(onePar,buffersize,"fM02Cut =%2.2f \n",        fM02Cut) ;
 138   parList+=onePar ;
 139   snprintf(onePar,buffersize,"fMinNCells =%d \n",        fMinNCells) ;
 140   parList+=onePar ;
 141   snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f\n", fMassPi0Min,fMassPi0Max);
 142   parList+=onePar ;
 143   snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f\n", fMassEtaMin,fMassEtaMax);
 144   parList+=onePar ;
 145   snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f\n", fMassConMin,fMassConMax);
 146   parList+=onePar ;
 147
 148   return new TObjString(parList) ;
 149
 150 }
 151
 152
 153 //_____________________________________________________________________________________
 154 TLorentzVector AliAnaInsideClusterInvariantMass::GetCellMomentum(const Int_t absId,
 155                                                                  Float_t en,
 156                                                                  AliVCaloCells * cells)
 157 {
 158
 159   // Cell momentum calculation for invariant mass
 160
 161   Double_t cellpos[] = {0, 0, 0};
 162   GetEMCALGeometry()->GetGlobal(absId, cellpos);
 163
 164   if(GetVertex(0)){//calculate direction from vertex
 165     cellpos[0]-=GetVertex(0)[0];
 166     cellpos[1]-=GetVertex(0)[1];
 167     cellpos[2]-=GetVertex(0)[2];
 168   }
 169
 170   Double_t r = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]+cellpos[2]*cellpos[2] ) ;
 171
 172   //If not calculated before, get the energy
 173   if(en <=0 )
 174   {
 175     en = cells->GetCellAmplitude(absId);
 176     GetCaloUtils()->RecalibrateCellAmplitude(en,fCalorimeter,absId);
 177   }
 178
 179   TLorentzVector cellMom ;
 180   cellMom.SetPxPyPzE( en*cellpos[0]/r,  en*cellpos[1]/r, en*cellpos[2]/r,  en) ;
 181
 182   return cellMom;
 183
 184 }
 185
 186 //________________________________________________________________
 187 TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
 188 {
 189   // Create histograms to be saved in output file and
 190   // store them in outputContainer
 191   TList * outputContainer = new TList() ;
 192   outputContainer->SetName("InsideClusterHistos") ;
 193
 194   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
 195   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
 196   Int_t mbins    = GetHistogramRanges()->GetHistoMassBins();         Float_t mmax   = GetHistogramRanges()->GetHistoMassMax();         Float_t mmin   = GetHistogramRanges()->GetHistoMassMin();
 197   Int_t ncbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   ncmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   ncmin  = GetHistogramRanges()->GetHistoNClusterCellMin();
 198
 199   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
 200   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
 201   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
 202   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
 203   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
 204   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
 205
 206   TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
 207   TString pname[] ={"","Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};
 208
 209   Int_t n = 1;
 210
 211   if(IsDataMC()) n = 7;
 212
 213   Int_t nMaxBins = 10;
 214
 215   TString sMatched[] = {"","Matched"};
 216
 217   for(Int_t i = 0; i < n; i++)
 218   {
 219
 220     for(Int_t j = 0; j < 2; j++)
 221     {
 222
 223       fhMassNLocMax1[i][j]  = new TH2F(Form("hMassNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 224                                        Form("Invariant mass of 2 highest energy cells %s %s",ptype[i].Data(),sMatched[j].Data()),
 225                                        nptbins,ptmin,ptmax,mbins,mmin,mmax);
 226       fhMassNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
 227       fhMassNLocMax1[i][j]->SetXTitle("E (GeV)");
 228       outputContainer->Add(fhMassNLocMax1[i][j]) ;
 229
 230       fhMassNLocMax2[i][j]  = new TH2F(Form("hMassNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 231                                        Form("Invariant mass of 2 local maxima cells %s %s",ptype[i].Data(),sMatched[j].Data()),
 232                                        nptbins,ptmin,ptmax,mbins,mmin,mmax);
 233       fhMassNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
 234       fhMassNLocMax2[i][j]->SetXTitle("E (GeV)");
 235       outputContainer->Add(fhMassNLocMax2[i][j]) ;
 236
 237       fhMassNLocMaxN[i][j]  = new TH2F(Form("hMassNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 238                                        Form("Invariant mass of N>2 local maxima cells %s %s",ptype[i].Data(),sMatched[j].Data()),
 239                                        nptbins,ptmin,ptmax,mbins,mmin,mmax);
 240       fhMassNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
 241       fhMassNLocMaxN[i][j]->SetXTitle("E (GeV)");
 242       outputContainer->Add(fhMassNLocMaxN[i][j]) ;
 243
 244
 245       fhNLocMax[i][j]     = new TH2F(Form("hNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
 246                                      Form("Number of local maxima in cluster %s %s",ptype[i].Data(),sMatched[j].Data()),
 247                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
 248       fhNLocMax[i][j]   ->SetYTitle("N maxima");
 249       fhNLocMax[i][j]   ->SetXTitle("E (GeV)");
 250       outputContainer->Add(fhNLocMax[i][j]) ;
 251
 252       fhNLocMaxEMax[i][j]     = new TH2F(Form("hNLocMaxEMax%s%s",pname[i].Data(),sMatched[j].Data()),
 253                                          Form("Number of local maxima in cluster vs energy of maxima %s %s",ptype[i].Data(),sMatched[j].Data()),
 254                                          nptbins*10,ptmin,ptmax,nMaxBins,0,nMaxBins);
 255       fhNLocMaxEMax[i][j]   ->SetYTitle("N maxima");
 256       fhNLocMaxEMax[i][j]   ->SetXTitle("E of maxima (GeV)");
 257       outputContainer->Add(fhNLocMaxEMax[i][j]) ;
 258
 259       fhNLocMaxEFrac[i][j]     = new TH2F(Form("hNLocMaxEFrac%s%s",pname[i].Data(),sMatched[j].Data()),
 260                                           Form("Number of local maxima in cluster vs fraction of cluster energy of maxima %s %s",ptype[i].Data(),sMatched[j].Data()),
 261                                           100,0,1,nMaxBins,0,nMaxBins);
 262       fhNLocMaxEFrac[i][j]   ->SetYTitle("N maxima");
 263       fhNLocMaxEFrac[i][j]   ->SetXTitle("E maxima / E cluster");
 264       outputContainer->Add(fhNLocMaxEFrac[i][j]) ;
 265
 266       fhNLocMaxM02Cut[i][j] = new TH2F(Form("hNLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
 267                                        Form("Number of local maxima in cluster %s for M02 > %2.2f",ptype[i].Data(),fM02Cut),
 268                                        nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
 269       fhNLocMaxM02Cut[i][j]->SetYTitle("N maxima");
 270       fhNLocMaxM02Cut[i][j]->SetXTitle("E (GeV)");
 271       outputContainer->Add(fhNLocMaxM02Cut[i][j]) ;
 272
 273
 274       fhM02NLocMax1[i][j]     = new TH2F(Form("hM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 275                                          Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
 276                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 277       fhM02NLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 278       fhM02NLocMax1[i][j]   ->SetXTitle("E (GeV)");
 279       outputContainer->Add(fhM02NLocMax1[i][j]) ;
 280
 281       fhM02NLocMax2[i][j]     = new TH2F(Form("hM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 282                                          Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 283                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 284       fhM02NLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 285       fhM02NLocMax2[i][j]   ->SetXTitle("E (GeV)");
 286       outputContainer->Add(fhM02NLocMax2[i][j]) ;
 287
 288
 289       fhM02NLocMaxN[i][j]    = new TH2F(Form("hM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 290                                         Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 291                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 292       fhM02NLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 293       fhM02NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
 294       outputContainer->Add(fhM02NLocMaxN[i][j]) ;
 295
 296
 297       fhNCellNLocMax1[i][j]  = new TH2F(Form("hNCellNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 298                                         Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
 299                                         nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
 300       fhNCellNLocMax1[i][j] ->SetYTitle("N cells");
 301       fhNCellNLocMax1[i][j] ->SetXTitle("E (GeV)");
 302       outputContainer->Add(fhNCellNLocMax1[i][j]) ;
 303
 304       fhNCellNLocMax2[i][j]     = new TH2F(Form("hNCellNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 305                                            Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 306                                            nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
 307       fhNCellNLocMax2[i][j]   ->SetYTitle("N cells");
 308       fhNCellNLocMax2[i][j]   ->SetXTitle("E (GeV)");
 309       outputContainer->Add(fhNCellNLocMax2[i][j]) ;
 310
 311
 312       fhNCellNLocMaxN[i][j]     = new TH2F(Form("hNCellNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 313                                            Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 314                                            nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
 315       fhNCellNLocMaxN[i][j]   ->SetYTitle("N cells");
 316       fhNCellNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
 317       outputContainer->Add(fhNCellNLocMaxN[i][j]) ;
 318
 319
 320       fhM02Pi0LocMax1[i][j]     = new TH2F(Form("hM02Pi0LocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 321                                            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()),
 322                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 323       fhM02Pi0LocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 324       fhM02Pi0LocMax1[i][j]   ->SetXTitle("E (GeV)");
 325       outputContainer->Add(fhM02Pi0LocMax1[i][j]) ;
 326
 327       fhM02EtaLocMax1[i][j]     = new TH2F(Form("hM02EtaLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 328                                            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()),
 329                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 330       fhM02EtaLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 331       fhM02EtaLocMax1[i][j]   ->SetXTitle("E (GeV)");
 332       outputContainer->Add(fhM02EtaLocMax1[i][j]) ;
 333
 334       fhM02ConLocMax1[i][j]    = new TH2F(Form("hM02ConLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 335                                           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()),
 336                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 337       fhM02ConLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 338       fhM02ConLocMax1[i][j]   ->SetXTitle("E (GeV)");
 339       outputContainer->Add(fhM02ConLocMax1[i][j]) ;
 340
 341       fhM02Pi0LocMax2[i][j]     = new TH2F(Form("hM02Pi0LocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 342                                            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()),
 343                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 344       fhM02Pi0LocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 345       fhM02Pi0LocMax2[i][j]   ->SetXTitle("E (GeV)");
 346       outputContainer->Add(fhM02Pi0LocMax2[i][j]) ;
 347
 348       fhM02EtaLocMax2[i][j]     = new TH2F(Form("hM02EtaLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 349                                            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()),
 350                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 351       fhM02EtaLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 352       fhM02EtaLocMax2[i][j]   ->SetXTitle("E (GeV)");
 353       outputContainer->Add(fhM02EtaLocMax2[i][j]) ;
 354
 355       fhM02ConLocMax2[i][j]    = new TH2F(Form("hM02ConLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 356                                           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()),
 357                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 358       fhM02ConLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 359       fhM02ConLocMax2[i][j]   ->SetXTitle("E (GeV)");
 360       outputContainer->Add(fhM02ConLocMax2[i][j]) ;
 361
 362       fhM02Pi0LocMaxN[i][j]     = new TH2F(Form("hM02Pi0LocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 363                                            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()),
 364                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 365       fhM02Pi0LocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 366       fhM02Pi0LocMaxN[i][j]   ->SetXTitle("E (GeV)");
 367       outputContainer->Add(fhM02Pi0LocMaxN[i][j]) ;
 368
 369       fhM02EtaLocMaxN[i][j]     = new TH2F(Form("hM02EtaLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 370                                            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()),
 371                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 372       fhM02EtaLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 373       fhM02EtaLocMaxN[i][j]   ->SetXTitle("E (GeV)");
 374       outputContainer->Add(fhM02EtaLocMaxN[i][j]) ;
 375
 376       fhM02ConLocMaxN[i][j]    = new TH2F(Form("hM02ConLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 377                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",fMassConMin,fMassConMax,ptype[i].Data()),
 378                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 379       fhM02ConLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 380       fhM02ConLocMaxN[i][j]   ->SetXTitle("E (GeV)");
 381       outputContainer->Add(fhM02ConLocMaxN[i][j]) ;
 382
 383     } // matched, not matched
 384
 385   } // MC particle list
 386
 387   for(Int_t i = 0; i < n; i++)
 388   {
 389     fhTrackMatchedDEtaLocMax1[i]  = new TH2F
 390     (Form("hTrackMatchedDEtaLocMax1%s",pname[i].Data()),
 391      Form("d#eta of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
 392      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
 393     fhTrackMatchedDEtaLocMax1[i]->SetYTitle("d#eta");
 394     fhTrackMatchedDEtaLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
 395
 396     fhTrackMatchedDPhiLocMax1[i]  = new TH2F
 397     (Form("hTrackMatchedDPhiLocMax1%s",pname[i].Data()),
 398      Form("d#phi of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
 399      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
 400     fhTrackMatchedDPhiLocMax1[i]->SetYTitle("d#phi (rad)");
 401     fhTrackMatchedDPhiLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
 402
 403     outputContainer->Add(fhTrackMatchedDEtaLocMax1[i]) ;
 404     outputContainer->Add(fhTrackMatchedDPhiLocMax1[i]) ;
 405
 406     fhTrackMatchedDEtaLocMax2[i]  = new TH2F
 407     (Form("hTrackMatchedDEtaLocMax2%s",pname[i].Data()),
 408      Form("d#eta of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
 409      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
 410     fhTrackMatchedDEtaLocMax2[i]->SetYTitle("d#eta");
 411     fhTrackMatchedDEtaLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
 412
 413     fhTrackMatchedDPhiLocMax2[i]  = new TH2F
 414     (Form("hTrackMatchedDPhiLocMax2%s",pname[i].Data()),
 415      Form("d#phi of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
 416      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
 417     fhTrackMatchedDPhiLocMax2[i]->SetYTitle("d#phi (rad)");
 418     fhTrackMatchedDPhiLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
 419
 420     outputContainer->Add(fhTrackMatchedDEtaLocMax2[i]) ;
 421     outputContainer->Add(fhTrackMatchedDPhiLocMax2[i]) ;
 422
 423     fhTrackMatchedDEtaLocMaxN[i]  = new TH2F
 424     (Form("hTrackMatchedDEtaLocMaxN%s",pname[i].Data()),
 425      Form("d#eta of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
 426      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
 427     fhTrackMatchedDEtaLocMaxN[i]->SetYTitle("d#eta");
 428     fhTrackMatchedDEtaLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
 429
 430     fhTrackMatchedDPhiLocMaxN[i]  = new TH2F
 431     (Form("hTrackMatchedDPhiLocMaxN%s",pname[i].Data()),
 432      Form("d#phi of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
 433      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
 434     fhTrackMatchedDPhiLocMaxN[i]->SetYTitle("d#phi (rad)");
 435     fhTrackMatchedDPhiLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
 436
 437     outputContainer->Add(fhTrackMatchedDEtaLocMaxN[i]) ;
 438     outputContainer->Add(fhTrackMatchedDPhiLocMaxN[i]) ;
 439
 440   }
 441
 442   for(Int_t j = 0; j < 2; j++)
 443   {
 444
 445     fhAnglePairLocMax1[j]  = new TH2F(Form("hAnglePairLocMax1%s",sMatched[j].Data()),
 446                                       Form("Opening angle of 2 highest energy cells vs pair Energy, %s",sMatched[j].Data()),
 447                                       nptbins,ptmin,ptmax,200,0,0.2);
 448     fhAnglePairLocMax1[j]->SetYTitle("#alpha (rad)");
 449     fhAnglePairLocMax1[j]->SetXTitle("E (GeV)");
 450     outputContainer->Add(fhAnglePairLocMax1[j]) ;
 451
 452     fhAnglePairLocMax2[j]  = new TH2F(Form("hAnglePairLocMax2%s",sMatched[j].Data()),
 453                                       Form("Opening angle of 2 local maxima cells vs Energy, %s",sMatched[j].Data()),
 454                                       nptbins,ptmin,ptmax,200,0,0.2);
 455     fhAnglePairLocMax2[j]->SetYTitle("#alpha (rad)");
 456     fhAnglePairLocMax2[j]->SetXTitle("E (GeV)");
 457     outputContainer->Add(fhAnglePairLocMax2[j]) ;
 458
 459     fhAnglePairLocMaxN[j]  = new TH2F(Form("hAnglePairLocMaxN%s",sMatched[j].Data()),
 460                                       Form("Opening angle of N>2 local maxima cells vs Energy, %s",sMatched[j].Data()),
 461                                       nptbins,ptmin,ptmax,200,0,0.2);
 462     fhAnglePairLocMaxN[j]->SetYTitle("#alpha (rad)");
 463     fhAnglePairLocMaxN[j]->SetXTitle("E (GeV)");
 464     outputContainer->Add(fhAnglePairLocMaxN[j]) ;
 465
 466     fhAnglePairMassLocMax1[j]  = new TH2F(Form("hAnglePairMassLocMax1%s",sMatched[j].Data()),
 467                                           Form("Opening angle of 2 highest energy cells vs Mass for E > 5 GeV, %s",sMatched[j].Data()),
 468                                           mbins,mmin,mmax,200,0,0.2);
 469     fhAnglePairMassLocMax1[j]->SetXTitle("M (GeV/c^{2})");
 470     fhAnglePairMassLocMax1[j]->SetYTitle("#alpha (rad)");
 471     outputContainer->Add(fhAnglePairMassLocMax1[j]) ;
 472
 473     fhAnglePairMassLocMax2[j]  = new TH2F(Form("hAnglePairMassLocMax2%s",sMatched[j].Data()),
 474                                           Form("Opening angle of 2 local maxima cells vs Mass for E > 5 GeV, %s",sMatched[j].Data()),
 475                                           mbins,mmin,mmax,200,0,0.2);
 476     fhAnglePairMassLocMax2[j]->SetXTitle("M (GeV/c^{2})");
 477     fhAnglePairMassLocMax2[j]->SetYTitle("#alpha (rad)");
 478     outputContainer->Add(fhAnglePairMassLocMax2[j]) ;
 479
 480     fhAnglePairMassLocMaxN[j]  = new TH2F(Form("hAnglePairMassLocMaxN%s",sMatched[j].Data()),
 481                                           Form("Opening angle of N>2 local maxima cells vs Mass for E > 5 GeV, %s",sMatched[j].Data()),
 482                                           mbins,mmin,mmax,200,0,0.2);
 483     fhAnglePairMassLocMaxN[j]->SetXTitle("M (GeV/c^{2})");
 484     fhAnglePairMassLocMaxN[j]->SetYTitle("#alpha (rad)");
 485     outputContainer->Add(fhAnglePairMassLocMaxN[j]) ;
 486
 487   }
 488
 489   return outputContainer ;
 490
 491 }
 492
 493 //___________________________________________
 494 void AliAnaInsideClusterInvariantMass::Init()
 495 {
 496   //Init
 497   //Do some checks
 498   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
 499     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 500     abort();
 501   }
 502   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
 503     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 504     abort();
 505   }
 506
 507   if( GetReader()->GetDataType() == AliCaloTrackReader::kMC ){
 508     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
 509     abort();
 510
 511   }
 512
 513 }
 514
 515 //_____________________________________________________
 516 void AliAnaInsideClusterInvariantMass::InitParameters()
 517 {
 518   //Initialize the parameters of the analysis.
 519   AddToHistogramsName("AnaPi0InsideClusterInvariantMass_");
 520
 521   fCalorimeter = "EMCAL" ;
 522
 523   fM02Cut      = 0.26 ;
 524   fMinNCells   = 4 ;
 525
 526   fMassEtaMin  = 0.4;
 527   fMassEtaMax  = 0.6;
 528
 529   fMassPi0Min  = 0.08;
 530   fMassPi0Max  = 0.20;
 531
 532   fMassConMin  = 0.0;
 533   fMassConMax  = 0.05;
 534
 535 }
 536
 537
 538 //__________________________________________________________________
 539 void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
 540 {
 541   //Search for pi0 in fCalorimeter with shower shape analysis
 542
 543   TObjArray * pl       = 0x0;
 544   AliVCaloCells* cells = 0x0;
 545
 546   //Select the Calorimeter of the photon
 547   if(fCalorimeter == "PHOS"){
 548     pl    = GetPHOSClusters();
 549     cells = GetPHOSCells();
 550   }
 551   else if (fCalorimeter == "EMCAL"){
 552     pl    = GetEMCALClusters();
 553     cells = GetEMCALCells();
 554   }
 555
 556   if(!pl || !cells) {
 557     Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
 558     return;
 559   }
 560
 561         if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
 562
 563   for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++){
 564     AliVCluster * cluster = (AliVCluster*) (pl->At(icluster));
 565
 566     // Study clusters with large shape parameter
 567     Float_t en = cluster->E();
 568     Float_t l0 = cluster->GetM02();
 569     Int_t   nc = cluster->GetNCells();
 570
 571     //If too small or big E or low number of cells, skip it
 572     if( en < GetMinEnergy() || en > GetMaxEnergy() || nc < fMinNCells) continue ;
 573
 574     Int_t    absId1    = -1; Int_t absId2 = -1;
 575     Int_t   *absIdList = new Int_t  [nc];
 576     Float_t *maxEList  = new Float_t[nc];
 577     Int_t    nMax      = GetCaloUtils()->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
 578     Bool_t   matched   = IsTrackMatched(cluster,GetReader()->GetInputEvent());
 579
 580     if (nMax <= 0)
 581     {
 582       printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found!\n");
 583
 584       /*
 585       for(Int_t iDigit  = 0; iDigit < cluster->GetNCells(); iDigit++ ) {
 586         Float_t ec = cells->GetCellAmplitude(cluster->GetCellsAbsId()[iDigit]);
 587         GetCaloUtils()->RecalibrateCellAmplitude(ec,fCalorimeter,cluster->GetCellsAbsId()[iDigit]);
 588         printf("iDigit %d, absId %d, Ecell %f\n",iDigit,cluster->GetCellsAbsId()[iDigit], ec);
 589       }
 590       */
 591
 592       delete [] absIdList ;
 593       delete [] maxEList  ;
 594       return;
 595     }
 596
 597     fhNLocMax[0][matched]->Fill(en,nMax);
 598     for(Int_t imax = 0; imax < nMax; imax++)
 599     {
 600       fhNLocMaxEMax [0][matched]->Fill(maxEList[imax]   ,nMax);
 601       fhNLocMaxEFrac[0][matched]->Fill(maxEList[imax]/en,nMax);
 602     }
 603
 604
 605     if     ( nMax == 1  ) { fhM02NLocMax1[0][matched]->Fill(en,l0) ; fhNCellNLocMax1[0][matched]->Fill(en,nc) ; }
 606     else if( nMax == 2  ) { fhM02NLocMax2[0][matched]->Fill(en,l0) ; fhNCellNLocMax2[0][matched]->Fill(en,nc) ; }
 607     else if( nMax >= 3  ) { fhM02NLocMaxN[0][matched]->Fill(en,l0) ; fhNCellNLocMaxN[0][matched]->Fill(en,nc) ; }
 608     else printf("N max smaller than 1 -> %d \n",nMax);
 609
 610     Float_t dZ  = cluster->GetTrackDz();
 611     Float_t dR  = cluster->GetTrackDx();
 612
 613     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
 614     {
 615       dR = 2000., dZ = 2000.;
 616       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
 617     }
 618     //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
 619
 620     if(TMath::Abs(dR) < 999)
 621     {
 622       if     ( nMax == 1  ) { fhTrackMatchedDEtaLocMax1[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMax1[0]->Fill(en,dR); }
 623       else if( nMax == 2  ) { fhTrackMatchedDEtaLocMax2[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMax2[0]->Fill(en,dR); }
 624       else if( nMax >= 3  ) { fhTrackMatchedDEtaLocMaxN[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMaxN[0]->Fill(en,dR); }
 625     }
 626
 627     // Play with the MC stack if available
 628     // Check origin of the candidates
 629     Int_t mcindex = -1;
 630     if(IsDataMC()){
 631
 632       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader(), 0);
 633
 634       if      ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )      mcindex = kmcPi0;
 635       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )      mcindex = kmcEta;
 636       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
 637                !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPhoton;
 638       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
 639                 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
 640       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))   mcindex = kmcElectron;
 641       else                                                                                mcindex = kmcHadron;
 642
 643       //GetMCAnalysisUtils()->PrintMCTag(tag);
 644       //printf("\t MC index Assigned %d \n",mcindex);
 645
 646       fhNLocMax[mcindex][matched]->Fill(en,nMax);
 647       for(Int_t imax = 0; imax < nMax; imax++)
 648       {
 649         fhNLocMaxEMax [mcindex][matched]->Fill(maxEList[imax]   ,nMax);
 650         fhNLocMaxEFrac[mcindex][matched]->Fill(maxEList[imax]/en,nMax);
 651       }
 652       if     (nMax == 1 ) { fhM02NLocMax1[mcindex][matched]->Fill(en,l0) ; fhNCellNLocMax1[mcindex][matched]->Fill(en,nc) ; }
 653       else if(nMax == 2 ) { fhM02NLocMax2[mcindex][matched]->Fill(en,l0) ; fhNCellNLocMax2[mcindex][matched]->Fill(en,nc) ; }
 654       else if(nMax >= 3 ) { fhM02NLocMaxN[mcindex][matched]->Fill(en,l0) ; fhNCellNLocMaxN[mcindex][matched]->Fill(en,nc) ; }
 655
 656       if(TMath::Abs(dR) < 999)
 657       {
 658         if     ( nMax == 1  ) { fhTrackMatchedDEtaLocMax1[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMax1[mcindex]->Fill(en,dR); }
 659         else if( nMax == 2  ) { fhTrackMatchedDEtaLocMax2[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMax2[mcindex]->Fill(en,dR); }
 660         else if( nMax >= 3  ) { fhTrackMatchedDEtaLocMaxN[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMaxN[mcindex]->Fill(en,dR); }
 661       }
 662
 663     }
 664
 665     //---------------------------------------------------------------------
 666     // From here only if M02 is large, fill histograms or split the cluster
 667     //---------------------------------------------------------------------
 668
 669     if( l0 < fM02Cut )
 670     {
 671       delete [] absIdList ;
 672       delete [] maxEList  ;
 673       continue;
 674     }
 675
 676     fhNLocMaxM02Cut[0][matched]->Fill(en,nMax);
 677     if(IsDataMC()) fhNLocMaxM02Cut[mcindex][matched]->Fill(en,nMax);
 678
 679     //---------------------------------------------------------------------
 680     // Get the 2 max indeces and do inv mass
 681     //---------------------------------------------------------------------
 682
 683     if     ( nMax == 2 ) {
 684       absId1 = absIdList[0];
 685       absId2 = absIdList[1];
 686     }
 687     else if( nMax == 1 )
 688     {
 689
 690       absId1 = absIdList[0];
 691
 692       //Find second highest energy cell
 693
 694       Float_t enmax = 0 ;
 695       for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++){
 696         Int_t absId = cluster->GetCellsAbsId()[iDigit];
 697         if( absId == absId1 ) continue ;
 698         Float_t endig = cells->GetCellAmplitude(absId);
 699         GetCaloUtils()->RecalibrateCellAmplitude(endig,fCalorimeter,absId);
 700         if(endig > enmax) {
 701           enmax  = endig ;
 702           absId2 = absId ;
 703         }
 704       }// cell loop
 705     }// 1 maxima
 706     else {  // n max > 2
 707       // loop on maxima, find 2 highest
 708
 709       // First max
 710       Float_t enmax = 0 ;
 711       for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
 712         Float_t endig = maxEList[iDigit];
 713         if(endig > enmax) {
 714           enmax  = endig ;
 715           absId1 = absIdList[iDigit];
 716         }
 717       }// first maxima loop
 718
 719       // Second max
 720       Float_t enmax2 = 0;
 721       for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
 722         if(absIdList[iDigit]==absId1) continue;
 723         Float_t endig = maxEList[iDigit];
 724         if(endig > enmax2) {
 725           enmax2  = endig ;
 726           absId2 = absIdList[iDigit];
 727         }
 728       }// second maxima loop
 729
 730     } // n local maxima > 2
 731
 732     //---------------------------------------------------------------------
 733     // Split the cluster energy in 2, around the highest 2 local maxima
 734     //---------------------------------------------------------------------
 735
 736     //Float_t en1 = 0, en2 = 0;
 737     //SplitEnergy(absId1,absId2,cluster, cells, en1, en2, nMax /*absIdList, maxEList,*/);
 738
 739     AliAODCaloCluster *cluster1 = new AliAODCaloCluster(0, 0,NULL,0.,NULL,NULL,1,0);
 740     AliAODCaloCluster *cluster2 = new AliAODCaloCluster(1, 0,NULL,0.,NULL,NULL,1,0);
 741
 742     SplitEnergy(absId1,absId2,cluster, cells, cluster1, cluster2, nMax /*absIdList, maxEList,*/);
 743
 744     //---------------------------------------------------------------------
 745     // Get mass of pair of clusters
 746     //---------------------------------------------------------------------
 747
 748     // First set position of cluster as local maxima position,
 749     // assign splitted energy to calculate momentum
 750
 751     //TLorentzVector cellMom1 = GetCellMomentum(absId1, en1, cells);
 752     //TLorentzVector cellMom2 = GetCellMomentum(absId2, en2, cells);
 753
 754     TLorentzVector cellMom1;
 755     TLorentzVector cellMom2;
 756
 757     cluster1->GetMomentum(cellMom1,GetVertex(0));
 758     cluster2->GetMomentum(cellMom2,GetVertex(0));
 759
 760     Float_t  mass  = (cellMom1+cellMom2).M();
 761     Double_t angle = cellMom2.Angle(cellMom1.Vect());
 762
 763     if     (nMax==1)
 764     {
 765       Int_t icol1 = -1, irow1 = -1, iRCU1 = -1;
 766       Int_t icol2 = -1, irow2 = -1, iRCU2 = -1;
 767       Int_t sm1 = GetCaloUtils()->GetModuleNumberCellIndexes(absId1, fCalorimeter, icol1, irow1, iRCU1) ;
 768       Int_t sm2 = GetCaloUtils()->GetModuleNumberCellIndexes(absId2, fCalorimeter, icol2, irow2, iRCU2) ;
 769       Float_t en1 = cluster1->E();
 770       Float_t en2 = cluster2->E();
 771
 772      printf("Angle %f, E %f, E1 %f, E2 %f, ecell1 %f, ecell2 %f, Asym %f, Fcell1 %f, Fcell 2 %f, SM1 %d, SM2 %d, icol1 %d, icol2 %d, irow1 %d, irow2 %d \n  ",
 773              angle,en,en1,en2,
 774              cells->GetCellAmplitude(absId1),cells->GetCellAmplitude(absId2),
 775              TMath::Abs(en1-en2)/(en1+en2),cells->GetCellAmplitude(absId1)/en, cells->GetCellAmplitude(absId2)/en,
 776              sm1,sm2,icol1,icol2,irow1,irow2);
 777
 778       fhAnglePairLocMax1[matched]->Fill(en,angle);
 779       if( en > 5 )
 780         fhAnglePairMassLocMax1[matched]->Fill(mass,angle);
 781
 782       fhMassNLocMax1[0][matched]->Fill(en,mass);
 783       if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax1[0][matched]->Fill(en,l0);
 784       else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax1[0][matched]->Fill(en,l0);
 785       else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax1[0][matched]->Fill(en,l0);
 786     }
 787     else if(nMax==2)
 788     {
 789       fhAnglePairLocMax2[matched]->Fill(en,angle);
 790       if( en > 5 )
 791         fhAnglePairMassLocMax2[matched]->Fill(mass,angle);
 792
 793       fhMassNLocMax2[0][matched]->Fill(en,mass);
 794       if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax2[0][matched]->Fill(en,l0);
 795       else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax2[0][matched]->Fill(en,l0);
 796       else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax2[0][matched]->Fill(en,l0);
 797     }
 798     else if(nMax >2)
 799     {
 800       fhAnglePairLocMaxN[matched]->Fill(en,angle);
 801       if( en > 5 )
 802         fhAnglePairMassLocMaxN[matched]->Fill(mass,angle);
 803
 804       fhMassNLocMaxN[0][matched]->Fill(en,mass);
 805       if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMaxN[0][matched]->Fill(en,l0);
 806       else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMaxN[0][matched]->Fill(en,l0);
 807       else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMaxN[0][matched]->Fill(en,l0);
 808     }
 809
 810     if(IsDataMC()){
 811
 812       if     (nMax==1)
 813       {
 814         fhMassNLocMax1[mcindex][matched]->Fill(en,mass);
 815         if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax1[mcindex][matched]->Fill(en,l0);
 816         else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax1[mcindex][matched]->Fill(en,l0);
 817         else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax1[mcindex][matched]->Fill(en,l0);
 818       }
 819       else if(nMax==2) {
 820         fhMassNLocMax2[mcindex][matched]->Fill(en,mass);
 821         if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax2[mcindex][matched]->Fill(en,l0);
 822         else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax2[mcindex][matched]->Fill(en,l0);
 823         else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax2[mcindex][matched]->Fill(en,l0);
 824       }
 825       else if(nMax >2) {
 826         fhMassNLocMaxN[mcindex][matched]->Fill(en,mass);
 827         if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMaxN[mcindex][matched]->Fill(en,l0);
 828         else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMaxN[mcindex][matched]->Fill(en,l0);
 829         else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMaxN[mcindex][matched]->Fill(en,l0);
 830       }
 831
 832     }//Work with MC truth first
 833
 834     delete cluster1 ;
 835     delete cluster2 ;
 836     delete [] absIdList ;
 837     delete [] maxEList  ;
 838
 839   }//loop
 840
 841   if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");
 842
 843 }
 844
 845 //______________________________________________________________________
 846 void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
 847 {
 848   //Print some relevant parameters set for the analysis
 849   if(! opt)
 850     return;
 851
 852   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 853   AliAnaCaloTrackCorrBaseClass::Print("");
 854   printf("Calorimeter     =     %s\n",  fCalorimeter.Data()) ;
 855   printf("Loc. Max. E > %2.2f\n",       GetCaloUtils()->GetLocalMaximaCutE());
 856   printf("Loc. Max. E Diff > %2.2f\n",  GetCaloUtils()->GetLocalMaximaCutEDiff());
 857   printf("lambda_0^2 >  %2.1f \n",      fM02Cut);
 858   printf("pi0 : %2.2f<m<%2.2f \n",      fMassPi0Min,fMassPi0Max);
 859   printf("eta : %2.2f<m<%2.2f \n",      fMassEtaMin,fMassEtaMax);
 860   printf("conv: %2.2f<m<%2.2f \n",      fMassConMin,fMassConMax);
 861
 862   printf("    \n") ;
 863
 864 }
 865
 866
 867
 868 //________________________________________________________________________________________
 869 void AliAnaInsideClusterInvariantMass::SplitEnergy(const Int_t absId1, const Int_t absId2,
 870                                                    AliVCluster* cluster,
 871                                                    AliVCaloCells* cells,
 872                                                    //Float_t & e1, Float_t & e2,
 873                                                    AliAODCaloCluster* cluster1,
 874                                                    AliAODCaloCluster* cluster2,
 875                                                    const Int_t nMax)
 876 {
 877
 878   // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
 879   // maxima are too close and have common cells, split the energy between the 2
 880
 881 /*
 882   TH2F* hClusterMap    = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
 883   TH2F* hClusterLocMax = new TH2F("hClusterLocMax","Cluster Local Maxima",48,0,48,24,0,24);
 884   TH2F* hCluster1      = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
 885   TH2F* hCluster2      = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
 886 */
 887
 888   const Int_t ncells  = cluster->GetNCells();
 889   Int_t absIdList[ncells];
 890
 891   Float_t e1 = 0,  e2   = 0 ;
 892   Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
 893
 894   //printf("Split Local Max: 1) %d - 2) %d\n",absId1,absId2);
 895   Float_t eCluster = 0;
 896   for(Int_t iDigit  = 0; iDigit < ncells; iDigit++ ) {
 897     absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
 898
 899     //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
 900
 901     sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList[iDigit], fCalorimeter, icol, irow, iRCU) ;
 902     Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
 903     GetCaloUtils()->RecalibrateCellAmplitude(ec,fCalorimeter, absIdList[iDigit]);
 904     eCluster+=ec;
 905
 906 /*    hClusterMap->Fill(icol,irow,ec); */
 907
 908   }
 909
 910   // Init counters and variables
 911   Int_t ncells1 = 1 ;
 912   UShort_t absIdList1[9] ;
 913   Double_t fracList1 [9] ;
 914   absIdList1[0] = absId1 ;
 915   fracList1 [0] = 1. ;
 916
 917   Float_t ecell1 = cells->GetCellAmplitude(absId1);
 918   GetCaloUtils()->RecalibrateCellAmplitude(ecell1, fCalorimeter, absId1);
 919   e1 =  ecell1;
 920
 921   Int_t ncells2 = 1 ;
 922   UShort_t absIdList2[9] ;
 923   Double_t fracList2 [9] ;
 924   absIdList2[0] = absId2 ;
 925   fracList2 [0] = 1. ;
 926
 927   Float_t ecell2 = cells->GetCellAmplitude(absId2);
 928   GetCaloUtils()->RecalibrateCellAmplitude(ecell2, fCalorimeter, absId2);
 929   e2 =  ecell2;
 930
 931   /*
 932   Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
 933   sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId1, fCalorimeter, icol1, irow1, iRCU) ;
 934   hClusterLocMax->Fill(icol1,irow1,ecell1);
 935   sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId2, fCalorimeter, icol2, irow2, iRCU) ;
 936   hClusterLocMax->Fill(icol2,irow2,ecell2);
 937   */
 938
 939   // Very rough way to share the cluster energy
 940   Float_t eRemain = (eCluster-ecell1-ecell2)/2;
 941   Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
 942   Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
 943
 944   for(Int_t iDigit = 0; iDigit < ncells; iDigit++){
 945     Int_t absId = absIdList[iDigit];
 946
 947     if(absId==absId1 || absId==absId2 || absId < 0) continue;
 948
 949     Float_t ecell = cells->GetCellAmplitude(absId);
 950     GetCaloUtils()->RecalibrateCellAmplitude(ecell, fCalorimeter, absId);
 951
 952     if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId1,absId ))
 953     {
 954        absIdList1[ncells1]= absId;
 955
 956       if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId2,absId ))
 957       {
 958         fracList1[ncells1] = shareFraction1;
 959         e1 += ecell*shareFraction1;
 960       }
 961       else
 962       {
 963         fracList1[ncells1] = 1.;
 964         e1 += ecell;
 965       }
 966
 967       ncells1++;
 968
 969      } // neigbour to cell1
 970
 971     if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId2,absId ))
 972     {
 973       absIdList2[ncells2]= absId;
 974
 975       if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId1,absId ))
 976       {
 977         fracList2[ncells2] = shareFraction2;
 978         e2 += ecell*shareFraction2;
 979       }
 980       else
 981       {
 982         fracList2[ncells2] = 1.;
 983         e2 += ecell;
 984       }
 985
 986       ncells2++;
 987
 988     } // neigbour to cell2
 989
 990   }
 991
 992    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",
 993          nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
 994
 995   cluster1->SetE(e1);
 996   cluster2->SetE(e2);
 997
 998   cluster1->SetNCells(ncells1);
 999   cluster2->SetNCells(ncells2);
1000
1001   cluster1->SetCellsAbsId(absIdList1);
1002   cluster2->SetCellsAbsId(absIdList2);
1003
1004   cluster1->SetCellsAmplitudeFraction(fracList1);
1005   cluster2->SetCellsAmplitudeFraction(fracList2);
1006
1007   GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1008   GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1009
1010
1011 /*
1012   printf("Cells of cluster1: ");
1013   for(Int_t iDigit  = 0; iDigit < ncells1; iDigit++ )
1014   {
1015     printf(" %d ",absIdList1[iDigit]);
1016
1017     sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList1[iDigit], fCalorimeter, icol, irow, iRCU) ;
1018
1019     if( GetCaloUtils()->AreNeighbours(fCalorimeter, absId2,absIdList1[iDigit]) )
1020       hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit])*shareFraction1);
1021     else
1022       hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
1023   }
1024
1025   printf(" \n ");
1026   printf("Cells of cluster2: ");
1027
1028   for(Int_t iDigit  = 0; iDigit < ncells2; iDigit++ )
1029   {
1030     printf(" %d ",absIdList2[iDigit]);
1031
1032     sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList2[iDigit], fCalorimeter, icol, irow, iRCU) ;
1033     if( GetCaloUtils()->AreNeighbours(fCalorimeter, absId1,absIdList2[iDigit]) )
1034       hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit])*shareFraction2);
1035     else
1036       hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
1037
1038   }
1039   printf(" \n ");
1040
1041    gStyle->SetPadRightMargin(0.15);
1042    gStyle->SetPadLeftMargin(0.1);
1043    gStyle->SetOptStat(0);
1044    gStyle->SetOptFit(000000);
1045
1046    TCanvas  * c= new TCanvas("canvas", "canvas", 8000, 4000) ;
1047    c->Divide(2,2);
1048    c->cd(1);
1049    gPad->SetGridy();
1050    gPad->SetGridx();
1051    hClusterMap->Draw("colz");
1052    c->cd(2);
1053    gPad->SetGridy();
1054    gPad->SetGridx();
1055    hClusterLocMax ->Draw("colz");
1056    c->cd(3);
1057    gPad->SetGridy();
1058    gPad->SetGridx();
1059    hCluster1      ->Draw("colz");
1060    c->cd(4);
1061    gPad->SetGridy();
1062    gPad->SetGridx();
1063    hCluster2      ->Draw("colz");
1064
1065    c->Print(Form("Event%d_nMax%d_NCell1_%d_NCell2_%d.eps",GetEventNumber(),nMax,ncells1,ncells2));
1066
1067    delete c;
1068    delete hClusterMap;
1069    delete hClusterLocMax;
1070    delete hCluster1;
1071    delete hCluster2;
1072 */
1073 }
1074
1075
1076 //________________________________________________________________________________________
1077 //void AliAnaInsideClusterInvariantMass::SplitEnergy(Int_t absId1, Int_t absId2,
1078 //                                                   AliVCluster* cluster,
1079 //                                                   AliVCaloCells* cells,
1080 //                                                   Float_t & e1, Float_t & e2,
1081 //                                                   const Int_t nMax, Int_t *listMax, Float_t *eListMax,)
1082 //{
1083 //
1084 //  // Split energy of cluster between the 2 local maxima.
1085 //  const Int_t ncells  = cluster->GetNCells();
1086 //  Int_t     absIdList[ncells];
1087 //  //Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
1088 ///*
1089 //  TH2F* hClusterMap    = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1090 //  TH2F* hClusterLocMax = new TH2F("hClusterLocMax","Cluster Local Maxima",48,0,48,24,0,24);
1091 //  TH2F* hClusterLocMax2= new TH2F("hClusterLocMax2","Cluster Highest Local Maxima",48,0,48,24,0,24);
1092 //  TH2F* hCluster1      = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1093 //  TH2F* hCluster2      = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1094 //  TH2F* hClusterNo     = new TH2F("hClusterNo","Cluster No",48,0,48,24,0,24);
1095 // */
1096 //  Float_t ec = 0;
1097 //  //printf("Split Local Max: 1) %d - 2) %d\n",absId1,absId2);
1098 //  for(Int_t iDigit  = 0; iDigit < ncells; iDigit++ ) {
1099 //    absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1100 //
1101 //    //sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList[iDigit], fCalorimeter, icol, irow, iRCU) ;
1102 //    //ec = cells->GetCellAmplitude(absIdList[iDigit]);
1103 //    //GetCaloUtils()->RecalibrateCellAmplitude(ec,fCalorimeter, absIdList[iDigit]);
1104 //    //hClusterMap->Fill(icol,irow,ec);
1105 //
1106 //    //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1107 //
1108 //  }
1109 //   /*
1110 //  printf("N Local Maxima %d \n",nMax);
1111 //  for(Int_t imax = 0; imax < nMax; imax++)
1112 //  {
1113 //    sm = GetCaloUtils()->GetModuleNumberCellIndexes(listMax[imax], fCalorimeter, icol, irow, iRCU) ;
1114 //    printf("LocalMaxima absId %d, Ecell %f\n",absIdList[imax], cells->GetCellAmplitude(listMax[imax]));
1115 //    hClusterLocMax->Fill(icol,irow,eListMax[imax]);
1116 //  }
1117 //  */
1118 //
1119 //  //Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1120 //  //sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId1, fCalorimeter, icol1, irow1, iRCU) ;
1121 //  Float_t ec1 = cells->GetCellAmplitude(absId1);
1122 //  GetCaloUtils()->RecalibrateCellAmplitude(ec1,fCalorimeter, absId1);
1123 //  //hClusterLocMax2->Fill(icol1,irow1,ec1);
1124 //
1125 //  //sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId2, fCalorimeter, icol2, irow2, iRCU) ;
1126 //  Float_t ec2 = cells->GetCellAmplitude(absId2);
1127 //  GetCaloUtils()->RecalibrateCellAmplitude(ec2,fCalorimeter, absId2);
1128 //  //hClusterLocMax2->Fill(icol2,irow2,ec2);
1129 //
1130 //  Int_t absIdtmp = 0;
1131 //  if(ec2>ec1){
1132 //    absIdtmp = absId2;
1133 //    absId2   = absId1;
1134 //    absId1   = absIdtmp;
1135 //  }
1136 //
1137 //  // SubCluster 1
1138 //
1139 //  // Init counters and variables
1140 //  Int_t ncells1 = 1;
1141 //  Int_t absIdList1[ncells];
1142 //  absIdList1[0] = absId1;
1143 //  //printf("First local max : absId1 %d %d \n",absId1, absIdList1[0]);
1144 //  for(Int_t iDigit1 = 1; iDigit1 < ncells; iDigit1++) absIdList1[iDigit1] = -1;
1145 //
1146 //  Float_t ecell1 = cells->GetCellAmplitude(absId1);
1147 //  GetCaloUtils()->RecalibrateCellAmplitude(ecell1,fCalorimeter, absId1);
1148 //  e1 =  ecell1;
1149 //
1150 //  //Int_t icolNew = -1, irowNew = -1, iRCUNew = -1;
1151 //  //Int_t jcol = -1, jrow = -1, jRCU = -1;
1152 //
1153 //  Bool_t added = kTRUE;
1154 //  Int_t cellj = 0;
1155 //  while (cellj < ncells1)
1156 //  {
1157 //    added = kFALSE;
1158 //    Int_t absId1New = absIdList1[cellj];
1159 //    //printf("\t absId %d added \n",absId1New);
1160 //
1161 //    Float_t e1New = cells->GetCellAmplitude(absId1New);
1162 //    GetCaloUtils()->RecalibrateCellAmplitude(e1New,fCalorimeter, absId1New);
1163 //
1164 //    //sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId1New, fCalorimeter, icolNew, irowNew, iRCUNew) ;
1165 //
1166 //    for(Int_t iDigit = 0; iDigit < ncells ; iDigit++)
1167 //    {
1168 //      Int_t absId = absIdList[iDigit] ;
1169 //      if(absId!=absId1New && absId!=absId2 && absId>=0)
1170 //      {
1171 //        Float_t en = cells->GetCellAmplitude(absId);
1172 //        GetCaloUtils()->RecalibrateCellAmplitude(en,fCalorimeter, absId);
1173 //        //printf("\t \t iDig %d, absId %d, absIdNew %d, en %f, enNew %f\n",iDigit,absId, absId1New,en, e1New);
1174 //        //sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId, fCalorimeter, jcol, jrow, jRCU) ;
1175 //        //printf("\t \t \t (col,row) New  (%d,%d), check (%d,%d) \n",icolNew, irowNew,jcol,jrow);
1176 //        if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId1New,absId )){
1177 //          //printf("\t \t \t neighbours\n");
1178 //          if(e1New > en-GetCaloUtils()->GetLocalMaximaCutEDiff()){
1179 //            absIdList1[ncells1++] = absId;
1180 //
1181 //            if((absId1New==absId1 && GetCaloUtils()->AreNeighbours(fCalorimeter, absId1,absId ) && GetCaloUtils()->AreNeighbours( absId2,absId ))) {
1182 //              e1+=en/2;
1183 //            }
1184 //            else {
1185 //              absIdList [iDigit]    = -1;
1186 //              e1+=en;
1187 //            }
1188 //          } // Decreasing energy with respect reference
1189 //        } // Neighbours
1190 //      } //Not local maxima or already removed
1191 //    } // cell loop
1192 //    cellj++;
1193 //  }// while cells added to list of cells for cl1
1194 //
1195 //  // SubCluster 2
1196 //
1197 //  // Init counters and variables
1198 //  Int_t ncells2 = 1;
1199 //  Int_t absIdList2[ncells];
1200 //  absIdList2[0] = absId2;
1201 //  //printf("Second local max : absId2 %d %d \n",absId2, absIdList2[0]);
1202 //  for(Int_t iDigit2 = 1; iDigit2 < ncells; iDigit2++) absIdList2[iDigit2] = -1;
1203 //
1204 //  Float_t ecell2 = cells->GetCellAmplitude(absId2);
1205 //  GetCaloUtils()->RecalibrateCellAmplitude(ecell2,fCalorimeter, absId2);
1206 //  e2 =  ecell2;
1207 //
1208 //  added = kTRUE;
1209 //  cellj = 0;
1210 //  while (cellj < ncells2)
1211 //  {
1212 //    added = kFALSE;
1213 //    Int_t absId2New = absIdList2[cellj];
1214 //    //printf("\t absId %d added \n",absId2New);
1215 //
1216 //    Float_t e2New = cells->GetCellAmplitude(absId2New);
1217 //    GetCaloUtils()->RecalibrateCellAmplitude(e2New,fCalorimeter,absId2New);
1218 //    //sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId2New, fCalorimeter, icolNew, irowNew, iRCU) ;
1219 //
1220 //    for(Int_t iDigit = 0; iDigit < ncells ; iDigit++)
1221 //    {
1222 //      Int_t absId = absIdList[iDigit] ;
1223 //      if(absId!=absId2New && absId>=0)
1224 //      {
1225 //        Float_t en = cells->GetCellAmplitude(absId);
1226 //        GetCaloUtils()->RecalibrateCellAmplitude(en,fCalorimeter, absId);
1227 //        //printf("\t \t iDig %d, absId %d, absIdNew %d, en %f, enNew %f\n",iDigit,absId, absId2New,en, e2New);
1228 //        //sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId, fCalorimeter, jcol, jrow, jRCU) ;
1229 //        //printf("\t \t \t (col,row) New  (%d,%d), check (%d,%d) \n",icolNew, irowNew,jcol,jrow);
1230 //        if(GetCaloUtils()->AreNeighbours( fCalorimeter, absId2New,absId )){
1231 //          //printf("\t \t \t neighbours\n");
1232 //          if(e2New > en-GetCaloUtils()->GetLocalMaximaCutEDiff()){
1233 //            absIdList2[ncells2++] = absId;
1234 //            absIdList [iDigit]    = -1;
1235 //            if(absId2New==absId2 && GetCaloUtils()->AreNeighbours(fCalorimeter, absId1,absId ) && GetCaloUtils()->AreNeighbours( fCalorimeter,absId2,absId )){
1236 //              e2+=en/2;
1237 //            }
1238 //            else {
1239 //              e2+=en;
1240 //            }
1241 //          } // Decreasing energy with respect reference
1242 //        } // Neighbours
1243 //      } //Not local maxima or already removed
1244 //    } // cell loop
1245 //    cellj++;
1246 //  }// while cells added to list of cells for cl2
1247 // /*
1248 //  for(Int_t iDigit  = 0; iDigit < ncells1; iDigit++ ) {
1249 //
1250 //    sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList1[iDigit], fCalorimeter, icol, irow, iRCU) ;
1251 //
1252 //    hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
1253 //
1254 //
1255 //  }
1256 //
1257 //  for(Int_t iDigit  = 0; iDigit < ncells2; iDigit++ ) {
1258 //
1259 //    sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList2[iDigit], fCalorimeter, icol, irow, iRCU) ;
1260 //
1261 //    hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
1262 //
1263 //
1264 //  }
1265 //
1266 //
1267 //  for(Int_t iDigit  = 0; iDigit < ncells; iDigit++ ) {
1268 //    if(absIdList[iDigit] < 0 || absIdList[iDigit]==absId1 || absIdList[iDigit]==absId2) continue;
1269 //    sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList[iDigit], fCalorimeter, icol, irow, iRCU) ;
1270 //    hClusterNo->Fill(icol,irow,cells->GetCellAmplitude(absIdList[iDigit]));
1271 //  }
1272 //
1273 //
1274 //  printf("Cluster energy  = %f, Ecell1 = %f, Ecell2 = %f, Enew1 = %f, Enew2 = %f, ncells %d, ncells1 %d, ncells2 %d \n",
1275 //         cluster->E(),ecell1,ecell2,e1,e2,ncells,ncells1,ncells2);
1276 //  if(ncells!=(ncells1+ncells2)) printf("\t Not all cells!\n");
1277 //
1278 //  gStyle->SetPadRightMargin(0.15);
1279 //  gStyle->SetPadLeftMargin(0.1);
1280 //  gStyle->SetOptStat(0);
1281 //  gStyle->SetOptFit(000000);
1282 //
1283 //  TCanvas  * c= new TCanvas("canvas", "canvas", 8000, 4000) ;
1284 //  c->Divide(3,2);
1285 //  c->cd(1);
1286 //  gPad->SetGridy();
1287 //  gPad->SetGridx();
1288 //  hClusterMap->Draw("colz");
1289 //  c->cd(2);
1290 //  gPad->SetGridy();
1291 //  gPad->SetGridx();
1292 //  hClusterLocMax ->Draw("colz");
1293 //  c->cd(3);
1294 //  gPad->SetGridy();
1295 //  gPad->SetGridx();
1296 //  hClusterLocMax2->Draw("colz");
1297 //  c->cd(4);
1298 //  gPad->SetGridy();
1299 //  gPad->SetGridx();
1300 //  hCluster1      ->Draw("colz");
1301 //  c->cd(5);
1302 //  gPad->SetGridy();
1303 //  gPad->SetGridx();
1304 //  hCluster2      ->Draw("colz");
1305 //  c->cd(6);
1306 //  gPad->SetGridy();
1307 //  gPad->SetGridx();
1308 //  hClusterNo     ->Draw("colz");
1309 //
1310 //  c->Print(Form("Event%d_nMax%d_NCell1_%d_NCell2_%d_Left%d.eps",GetEventNumber(),nMax,ncells1,ncells2,ncells-ncells1-ncells2));
1311 //
1312 //  delete c;
1313 //  delete hClusterMap;
1314 //  delete hClusterLocMax;
1315 //  delete hClusterLocMax2;
1316 //  delete hCluster1;
1317 //  delete hCluster2;
1318 //  delete hClusterNo;
1319 //*/
1320 //}
1321 //
1322