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