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 <TH2F.h>
  33
  34 // --- Analysis system ---
  35 #include "AliAnaInsideClusterInvariantMass.h"
  36 #include "AliCaloTrackReader.h"
  37 #include "AliMCAnalysisUtils.h"
  38 #include "AliStack.h"
  39 #include "AliFiducialCut.h"
  40 #include "TParticle.h"
  41 #include "AliVCluster.h"
  42 #include "AliAODEvent.h"
  43 #include "AliAODMCParticle.h"
  44 #include "AliEMCALGeoParams.h"
  45
  46 // --- Detectors ---
  47 //#include "AliPHOSGeoUtils.h"
  48 #include "AliEMCALGeometry.h"
  49
  50 ClassImp(AliAnaInsideClusterInvariantMass)
  51
  52 //__________________________________________________________________
  53 AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
  54   AliAnaCaloTrackCorrBaseClass(),
  55   fCalorimeter(""),
  56   fM02MaxCut(0),    fM02MinCut(0),
  57   fMinNCells(0),    fMinBadDist(0),
  58   fFillAngleHisto(kFALSE),
  59   fFillTMResidualHisto(kFALSE),
  60   fFillSSExtraHisto(kFALSE)
  61 {
  62   //default ctor
  63
  64   // Init array of histograms
  65   for(Int_t i = 0; i < 7; i++)
  66   {
  67     for(Int_t j = 0; j < 2; j++)
  68     {
  69       fhMassNLocMax1[i][j]  = 0;
  70       fhMassNLocMax2[i][j]  = 0;
  71       fhMassNLocMaxN[i][j]  = 0;
  72       fhNLocMax[i][j]       = 0;
  73       fhNLocMaxM02Cut[i][j] = 0;
  74       fhM02NLocMax1[i][j]   = 0;
  75       fhM02NLocMax2[i][j]   = 0;
  76       fhM02NLocMaxN[i][j]   = 0;
  77       fhNCellNLocMax1[i][j] = 0;
  78       fhNCellNLocMax2[i][j] = 0;
  79       fhNCellNLocMaxN[i][j] = 0;
  80       fhM02Pi0LocMax1[i][j] = 0;
  81       fhM02EtaLocMax1[i][j] = 0;
  82       fhM02ConLocMax1[i][j] = 0;
  83       fhM02Pi0LocMax2[i][j] = 0;
  84       fhM02EtaLocMax2[i][j] = 0;
  85       fhM02ConLocMax2[i][j] = 0;
  86       fhM02Pi0LocMaxN[i][j] = 0;
  87       fhM02EtaLocMaxN[i][j] = 0;
  88       fhM02ConLocMaxN[i][j] = 0;
  89
  90       fhMassM02NLocMax1[i][j]= 0;
  91       fhMassM02NLocMax2[i][j]= 0;
  92       fhMassM02NLocMaxN[i][j]= 0;
  93       fhMassDispEtaNLocMax1[i][j]= 0;
  94       fhMassDispEtaNLocMax2[i][j]= 0;
  95       fhMassDispEtaNLocMaxN[i][j]= 0;
  96       fhMassDispPhiNLocMax1[i][j]= 0;
  97       fhMassDispPhiNLocMax2[i][j]= 0;
  98       fhMassDispPhiNLocMaxN[i][j]= 0;
  99       fhMassDispAsyNLocMax1[i][j]= 0;
 100       fhMassDispAsyNLocMax2[i][j]= 0;
 101       fhMassDispAsyNLocMaxN[i][j]= 0;
 102
 103       fhSplitEFractionNLocMax1[i][j]=0;
 104       fhSplitEFractionNLocMax2[i][j]=0;
 105       fhSplitEFractionNLocMaxN[i][j]=0;
 106
 107       fhMCGenFracNLocMax1[i][j]= 0;
 108       fhMCGenFracNLocMax2[i][j]= 0;
 109       fhMCGenFracNLocMaxN[i][j]= 0;
 110     }
 111
 112     for(Int_t jj = 0; jj < 4; jj++)
 113     {
 114       fhM02MCGenFracNLocMax1Ebin[i][jj] = 0;
 115       fhM02MCGenFracNLocMax2Ebin[i][jj] = 0;
 116       fhM02MCGenFracNLocMaxNEbin[i][jj] = 0;
 117
 118       fhMassMCGenFracNLocMax1Ebin[i][jj]= 0;
 119       fhMassMCGenFracNLocMax2Ebin[i][jj]= 0;
 120       fhMassMCGenFracNLocMaxNEbin[i][jj]= 0;
 121
 122       fhMCGenFracNLocMaxEbin[i][jj]       = 0;
 123       fhMCGenFracNLocMaxEbinMatched[i][jj]= 0;
 124
 125       fhMassSplitEFractionNLocMax1Ebin[i][jj] = 0;
 126       fhMassSplitEFractionNLocMax2Ebin[i][jj] = 0;
 127       fhMassSplitEFractionNLocMaxNEbin[i][jj] = 0;
 128     }
 129
 130     fhTrackMatchedDEtaLocMax1[i] = 0;
 131     fhTrackMatchedDPhiLocMax1[i] = 0;
 132     fhTrackMatchedDEtaLocMax2[i] = 0;
 133     fhTrackMatchedDPhiLocMax2[i] = 0;
 134     fhTrackMatchedDEtaLocMaxN[i] = 0;
 135     fhTrackMatchedDPhiLocMaxN[i] = 0;
 136
 137   }
 138
 139   for(Int_t i = 0; i < 2; i++)
 140   {
 141     fhAnglePairLocMax1    [i] = 0;
 142     fhAnglePairLocMax2    [i] = 0;
 143     fhAnglePairLocMaxN    [i] = 0;
 144     fhAnglePairMassLocMax1[i] = 0;
 145     fhAnglePairMassLocMax2[i] = 0;
 146     fhAnglePairMassLocMaxN[i] = 0;
 147   }
 148
 149   for(Int_t i = 0; i < 4; i++)
 150   {
 151     fhMassM02NLocMax1Ebin[i] = 0 ;
 152     fhMassM02NLocMax2Ebin[i] = 0 ;
 153     fhMassM02NLocMaxNEbin[i] = 0 ;
 154
 155     fhMassDispEtaNLocMax1Ebin[i] = 0 ;
 156     fhMassDispEtaNLocMax2Ebin[i] = 0 ;
 157     fhMassDispEtaNLocMaxNEbin[i] = 0 ;
 158
 159     fhMassDispPhiNLocMax1Ebin[i] = 0 ;
 160     fhMassDispPhiNLocMax2Ebin[i] = 0 ;
 161     fhMassDispPhiNLocMaxNEbin[i] = 0 ;
 162
 163     fhMassDispAsyNLocMax1Ebin[i] = 0 ;
 164     fhMassDispAsyNLocMax2Ebin[i] = 0 ;
 165     fhMassDispAsyNLocMaxNEbin[i] = 0 ;
 166
 167     fhMCAsymM02NLocMax1MCPi0Ebin[i] = 0 ;
 168     fhMCAsymM02NLocMax2MCPi0Ebin[i] = 0 ;
 169     fhMCAsymM02NLocMaxNMCPi0Ebin[i] = 0 ;
 170   }
 171
 172   InitParameters();
 173
 174 }
 175
 176 //_______________________________________________________________
 177 TObjString *  AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
 178 {
 179         //Save parameters used for analysis
 180   TString parList ; //this will be list of parameters used for this analysis.
 181   const Int_t buffersize = 255;
 182   char onePar[buffersize] ;
 183
 184   snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
 185   parList+=onePar ;
 186
 187   snprintf(onePar,buffersize,"Calorimeter: %s\n",        fCalorimeter.Data()) ;
 188   parList+=onePar ;
 189   snprintf(onePar,buffersize,"fLocMaxCutE =%2.2f \n",    GetCaloUtils()->GetLocalMaximaCutE()) ;
 190   parList+=onePar ;
 191   snprintf(onePar,buffersize,"fLocMaxCutEDiff =%2.2f \n",GetCaloUtils()->GetLocalMaximaCutEDiff()) ;
 192   parList+=onePar ;
 193   snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n",    fM02MinCut, fM02MaxCut) ;
 194   parList+=onePar ;
 195   snprintf(onePar,buffersize,"fMinNCells =%d \n",        fMinNCells) ;
 196   parList+=onePar ;
 197   snprintf(onePar,buffersize,"fMinBadDist =%1.1f \n",    fMinBadDist) ;
 198   parList+=onePar ;
 199
 200   return new TObjString(parList) ;
 201
 202 }
 203
 204 //________________________________________________________________
 205 TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
 206 {
 207   // Create histograms to be saved in output file and
 208   // store them in outputContainer
 209   TList * outputContainer = new TList() ;
 210   outputContainer->SetName("InsideClusterHistos") ;
 211
 212   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
 213   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
 214   Int_t mbins    = GetHistogramRanges()->GetHistoMassBins();         Float_t mmax   = GetHistogramRanges()->GetHistoMassMax();         Float_t mmin   = GetHistogramRanges()->GetHistoMassMin();
 215   Int_t ncbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   ncmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   ncmin  = GetHistogramRanges()->GetHistoNClusterCellMin();
 216
 217   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
 218   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
 219   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
 220   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
 221   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
 222   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
 223
 224   TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
 225   TString pname[] ={"","Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};
 226
 227   Int_t n = 1;
 228
 229   if(IsDataMC()) n = 7;
 230
 231   Int_t nMaxBins = 10;
 232
 233   TString sMatched[] = {"","Matched"};
 234
 235   for(Int_t i = 0; i < n; i++)
 236   {
 237     for(Int_t j = 0; j < 2; j++)
 238     {
 239
 240       fhMassNLocMax1[i][j]  = new TH2F(Form("hMassNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 241                                        Form("Invariant mass of 2 highest energy cells vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
 242                                        nptbins,ptmin,ptmax,mbins,mmin,mmax);
 243       fhMassNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
 244       fhMassNLocMax1[i][j]->SetXTitle("E (GeV)");
 245       outputContainer->Add(fhMassNLocMax1[i][j]) ;
 246
 247       fhMassNLocMax2[i][j]  = new TH2F(Form("hMassNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 248                                        Form("Invariant mass of 2 local maxima cells vs E,%s %s",ptype[i].Data(),sMatched[j].Data()),
 249                                        nptbins,ptmin,ptmax,mbins,mmin,mmax);
 250       fhMassNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
 251       fhMassNLocMax2[i][j]->SetXTitle("E (GeV)");
 252       outputContainer->Add(fhMassNLocMax2[i][j]) ;
 253
 254       fhMassNLocMaxN[i][j]  = new TH2F(Form("hMassNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 255                                        Form("Invariant mass of N>2 local maxima cells vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
 256                                        nptbins,ptmin,ptmax,mbins,mmin,mmax);
 257       fhMassNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
 258       fhMassNLocMaxN[i][j]->SetXTitle("E (GeV)");
 259       outputContainer->Add(fhMassNLocMaxN[i][j]) ;
 260
 261       fhMassM02NLocMax1[i][j]  = new TH2F(Form("hMassM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 262                                           Form("Invariant mass of 2 highest energy cells #lambda_{0}^{2}, E > 7 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
 263                                           ssbins,ssmin,ssmax,mbins,mmin,mmax);
 264       fhMassM02NLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
 265       fhMassM02NLocMax1[i][j]->SetXTitle("#lambda_{0}^{2}");
 266       outputContainer->Add(fhMassM02NLocMax1[i][j]) ;
 267
 268       fhMassM02NLocMax2[i][j]  = new TH2F(Form("hMassM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 269                                           Form("Invariant mass of 2 local maxima cells #lambda_{0}^{2}, E > 7 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
 270                                           ssbins,ssmin,ssmax,mbins,mmin,mmax);
 271       fhMassM02NLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
 272       fhMassM02NLocMax2[i][j]->SetXTitle("#lambda_{0}^{2}");
 273       outputContainer->Add(fhMassM02NLocMax2[i][j]) ;
 274
 275       fhMassM02NLocMaxN[i][j]  = new TH2F(Form("hMassM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 276                                           Form("Invariant mass of N>2 local maxima cells vs #lambda_{0}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
 277                                           ssbins,ssmin,ssmax,mbins,mmin,mmax);
 278       fhMassM02NLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
 279       fhMassM02NLocMaxN[i][j]->SetXTitle("#lambda_{0}^{2}");
 280       outputContainer->Add(fhMassM02NLocMaxN[i][j]) ;
 281
 282       if(fFillSSExtraHisto)
 283       {
 284         fhMassDispEtaNLocMax1[i][j]  = new TH2F(Form("hMassDispEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 285                                                 Form("Invariant mass of 2 highest energy cells #sigma_{#eta #eta}^{2}, E > 7 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
 286                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
 287         fhMassDispEtaNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
 288         fhMassDispEtaNLocMax1[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
 289         outputContainer->Add(fhMassDispEtaNLocMax1[i][j]) ;
 290
 291         fhMassDispEtaNLocMax2[i][j]  = new TH2F(Form("hMassDispEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 292                                                 Form("Invariant mass of 2 local maxima cells #sigma_{#eta #eta}^{2}, E > 7 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
 293                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
 294         fhMassDispEtaNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
 295         fhMassDispEtaNLocMax2[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
 296         outputContainer->Add(fhMassDispEtaNLocMax2[i][j]) ;
 297
 298         fhMassDispEtaNLocMaxN[i][j]  = new TH2F(Form("hMassDispEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 299                                                 Form("Invariant mass of N>2 local maxima cells vs #sigma_{#eta #eta}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
 300                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
 301         fhMassDispEtaNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
 302         fhMassDispEtaNLocMaxN[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
 303         outputContainer->Add(fhMassDispEtaNLocMaxN[i][j]) ;
 304
 305         fhMassDispPhiNLocMax1[i][j]  = new TH2F(Form("hMassDispPhiNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 306                                                 Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E > 7 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
 307                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
 308         fhMassDispPhiNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
 309         fhMassDispPhiNLocMax1[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
 310         outputContainer->Add(fhMassDispPhiNLocMax1[i][j]) ;
 311
 312         fhMassDispPhiNLocMax2[i][j]  = new TH2F(Form("hMassDispPhiNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 313                                                 Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E > 7 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
 314                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
 315         fhMassDispPhiNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
 316         fhMassDispPhiNLocMax2[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
 317         outputContainer->Add(fhMassDispPhiNLocMax2[i][j]) ;
 318
 319         fhMassDispPhiNLocMaxN[i][j]  = new TH2F(Form("hMassDispPhiNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 320                                                 Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
 321                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
 322         fhMassDispPhiNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
 323         fhMassDispPhiNLocMaxN[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
 324         outputContainer->Add(fhMassDispPhiNLocMaxN[i][j]) ;
 325
 326         fhMassDispAsyNLocMax1[i][j]  = new TH2F(Form("hMassDispAsyNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 327                                                 Form("Invariant mass of 2 highest energy cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E > 7 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
 328                                                 200,-1,1,mbins,mmin,mmax);
 329         fhMassDispAsyNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
 330         fhMassDispAsyNLocMax1[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
 331         outputContainer->Add(fhMassDispAsyNLocMax1[i][j]) ;
 332
 333         fhMassDispAsyNLocMax2[i][j]  = new TH2F(Form("hMassDispAsyNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 334                                                 Form("Invariant mass of 2 local maxima cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E > 7 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
 335                                                 200,-1,1,mbins,mmin,mmax);
 336         fhMassDispAsyNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
 337         fhMassDispAsyNLocMax2[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
 338         outputContainer->Add(fhMassDispAsyNLocMax2[i][j]) ;
 339
 340         fhMassDispAsyNLocMaxN[i][j]  = new TH2F(Form("hMassDispAsyNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 341                                                 Form("Invariant mass of N>2 local maxima cells vsA = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), %s %s",ptype[i].Data(),sMatched[j].Data()),
 342                                                 200,-1,1,mbins,mmin,mmax);
 343         fhMassDispAsyNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
 344         fhMassDispAsyNLocMaxN[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
 345         outputContainer->Add(fhMassDispAsyNLocMaxN[i][j]) ;
 346       }
 347
 348       fhNLocMax[i][j]     = new TH2F(Form("hNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
 349                                      Form("Number of local maxima in cluster %s %s",ptype[i].Data(),sMatched[j].Data()),
 350                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
 351       fhNLocMax[i][j]   ->SetYTitle("N maxima");
 352       fhNLocMax[i][j]   ->SetXTitle("E (GeV)");
 353       outputContainer->Add(fhNLocMax[i][j]) ;
 354
 355       fhNLocMaxM02Cut[i][j] = new TH2F(Form("hNLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
 356                                        Form("Number of local maxima in cluster %s for %2.2f < M02 < %2.2f",ptype[i].Data(),fM02MinCut,fM02MaxCut),
 357                                        nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
 358       fhNLocMaxM02Cut[i][j]->SetYTitle("N maxima");
 359       fhNLocMaxM02Cut[i][j]->SetXTitle("E (GeV)");
 360       outputContainer->Add(fhNLocMaxM02Cut[i][j]) ;
 361
 362
 363       fhM02NLocMax1[i][j]     = new TH2F(Form("hM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 364                                          Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
 365                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 366       fhM02NLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 367       fhM02NLocMax1[i][j]   ->SetXTitle("E (GeV)");
 368       outputContainer->Add(fhM02NLocMax1[i][j]) ;
 369
 370       fhM02NLocMax2[i][j]     = new TH2F(Form("hM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 371                                          Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 372                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 373       fhM02NLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 374       fhM02NLocMax2[i][j]   ->SetXTitle("E (GeV)");
 375       outputContainer->Add(fhM02NLocMax2[i][j]) ;
 376
 377       fhM02NLocMaxN[i][j]    = new TH2F(Form("hM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 378                                         Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 379                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 380       fhM02NLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 381       fhM02NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
 382       outputContainer->Add(fhM02NLocMaxN[i][j]) ;
 383
 384
 385       fhSplitEFractionNLocMax1[i][j]     = new TH2F(Form("hSplitEFractionNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 386                                          Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
 387                                          nptbins,ptmin,ptmax,120,0,1.2);
 388       fhSplitEFractionNLocMax1[i][j]   ->SetXTitle("E_{cluster} (GeV)");
 389       fhSplitEFractionNLocMax1[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
 390       outputContainer->Add(fhSplitEFractionNLocMax1[i][j]) ;
 391
 392       fhSplitEFractionNLocMax2[i][j]     = new TH2F(Form("hSplitEFractionNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 393                                          Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 394                                          nptbins,ptmin,ptmax,120,0,1.2);
 395       fhSplitEFractionNLocMax2[i][j]   ->SetXTitle("E_{cluster} (GeV)");
 396       fhSplitEFractionNLocMax2[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
 397       outputContainer->Add(fhSplitEFractionNLocMax2[i][j]) ;
 398
 399       fhSplitEFractionNLocMaxN[i][j]    = new TH2F(Form("hSplitEFractionNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 400                                         Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 401                                         nptbins,ptmin,ptmax,120,0,1.2);
 402       fhSplitEFractionNLocMaxN[i][j]   ->SetXTitle("E_{cluster} (GeV)");
 403       fhSplitEFractionNLocMaxN[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
 404       outputContainer->Add(fhSplitEFractionNLocMaxN[i][j]) ;
 405
 406
 407       if(i > 0) // skip first entry in array, general case not filled
 408       {
 409         fhMCGenFracNLocMax1[i][j]     = new TH2F(Form("hMCGenFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 410                                                  Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
 411                                                  nptbins,ptmin,ptmax,200,0,2);
 412         fhMCGenFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
 413         fhMCGenFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
 414         outputContainer->Add(fhMCGenFracNLocMax1[i][j]) ;
 415
 416         fhMCGenFracNLocMax2[i][j]     = new TH2F(Form("hMCGenFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 417                                                  Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 418                                                  nptbins,ptmin,ptmax,200,0,2);
 419         fhMCGenFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
 420         fhMCGenFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
 421         outputContainer->Add(fhMCGenFracNLocMax2[i][j]) ;
 422
 423
 424         fhMCGenFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 425                                                 Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 426                                                 nptbins,ptmin,ptmax,200,0,2);
 427         fhMCGenFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
 428         fhMCGenFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
 429         outputContainer->Add(fhMCGenFracNLocMaxN[i][j]) ;
 430
 431       }
 432
 433       if(fFillSSExtraHisto)
 434       {
 435         fhNCellNLocMax1[i][j]  = new TH2F(Form("hNCellNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 436                                           Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
 437                                           nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
 438         fhNCellNLocMax1[i][j] ->SetYTitle("N cells");
 439         fhNCellNLocMax1[i][j] ->SetXTitle("E (GeV)");
 440         outputContainer->Add(fhNCellNLocMax1[i][j]) ;
 441
 442         fhNCellNLocMax2[i][j]     = new TH2F(Form("hNCellNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 443                                              Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 444                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
 445         fhNCellNLocMax2[i][j]   ->SetYTitle("N cells");
 446         fhNCellNLocMax2[i][j]   ->SetXTitle("E (GeV)");
 447         outputContainer->Add(fhNCellNLocMax2[i][j]) ;
 448
 449
 450         fhNCellNLocMaxN[i][j]     = new TH2F(Form("hNCellNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 451                                              Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
 452                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
 453         fhNCellNLocMaxN[i][j]   ->SetYTitle("N cells");
 454         fhNCellNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
 455         outputContainer->Add(fhNCellNLocMaxN[i][j]) ;
 456       }
 457
 458       fhM02Pi0LocMax1[i][j]     = new TH2F(Form("hM02Pi0LocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 459                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 1",
 460                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
 461                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 462       fhM02Pi0LocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 463       fhM02Pi0LocMax1[i][j]   ->SetXTitle("E (GeV)");
 464       outputContainer->Add(fhM02Pi0LocMax1[i][j]) ;
 465
 466       fhM02EtaLocMax1[i][j]     = new TH2F(Form("hM02EtaLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 467                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
 468                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
 469                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 470       fhM02EtaLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 471       fhM02EtaLocMax1[i][j]   ->SetXTitle("E (GeV)");
 472       outputContainer->Add(fhM02EtaLocMax1[i][j]) ;
 473
 474       fhM02ConLocMax1[i][j]    = new TH2F(Form("hM02ConLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
 475                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
 476                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
 477                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 478       fhM02ConLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 479       fhM02ConLocMax1[i][j]   ->SetXTitle("E (GeV)");
 480       outputContainer->Add(fhM02ConLocMax1[i][j]) ;
 481
 482       fhM02Pi0LocMax2[i][j]     = new TH2F(Form("hM02Pi0LocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 483                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 2",
 484                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
 485                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 486       fhM02Pi0LocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 487       fhM02Pi0LocMax2[i][j]   ->SetXTitle("E (GeV)");
 488       outputContainer->Add(fhM02Pi0LocMax2[i][j]) ;
 489
 490       fhM02EtaLocMax2[i][j]     = new TH2F(Form("hM02EtaLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 491                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
 492                                                GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
 493                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 494       fhM02EtaLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 495       fhM02EtaLocMax2[i][j]   ->SetXTitle("E (GeV)");
 496       outputContainer->Add(fhM02EtaLocMax2[i][j]) ;
 497
 498       fhM02ConLocMax2[i][j]    = new TH2F(Form("hM02ConLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
 499                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
 500                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
 501                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 502       fhM02ConLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 503       fhM02ConLocMax2[i][j]   ->SetXTitle("E (GeV)");
 504       outputContainer->Add(fhM02ConLocMax2[i][j]) ;
 505
 506       fhM02Pi0LocMaxN[i][j]     = new TH2F(Form("hM02Pi0LocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 507                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max > 2",
 508                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
 509                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 510       fhM02Pi0LocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 511       fhM02Pi0LocMaxN[i][j]   ->SetXTitle("E (GeV)");
 512       outputContainer->Add(fhM02Pi0LocMaxN[i][j]) ;
 513
 514       fhM02EtaLocMaxN[i][j]     = new TH2F(Form("hM02EtaLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 515                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max > 2",
 516                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
 517                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 518       fhM02EtaLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 519       fhM02EtaLocMaxN[i][j]   ->SetXTitle("E (GeV)");
 520       outputContainer->Add(fhM02EtaLocMaxN[i][j]) ;
 521
 522       fhM02ConLocMaxN[i][j]    = new TH2F(Form("hM02ConLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
 523                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
 524                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
 525                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 526       fhM02ConLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 527       fhM02ConLocMaxN[i][j]   ->SetXTitle("E (GeV)");
 528       outputContainer->Add(fhM02ConLocMaxN[i][j]) ;
 529
 530     } // matched, not matched
 531
 532       for(Int_t j = 0; j < 4; j++)
 533       {
 534
 535         fhMassSplitEFractionNLocMax1Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax1%sEbin%d",pname[i].Data(),j),
 536                                                            Form("Invariant mass of 2 highest energy cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
 537                                                            120,0,1.2,mbins,mmin,mmax);
 538         fhMassSplitEFractionNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
 539         fhMassSplitEFractionNLocMax1Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
 540         outputContainer->Add(fhMassSplitEFractionNLocMax1Ebin[i][j]) ;
 541
 542         fhMassSplitEFractionNLocMax2Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax2%sEbin%d",pname[i].Data(),j),
 543                                                            Form("Invariant mass of 2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
 544                                                            120,0,1.2,mbins,mmin,mmax);
 545         fhMassSplitEFractionNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
 546         fhMassSplitEFractionNLocMax2Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
 547         outputContainer->Add(fhMassSplitEFractionNLocMax2Ebin[i][j]) ;
 548
 549         fhMassSplitEFractionNLocMaxNEbin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMaxN%sEbin%d",pname[i].Data(),j),
 550                                                            Form("Invariant mass of N>2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
 551                                                            120,0,1.2,mbins,mmin,mmax);
 552         fhMassSplitEFractionNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
 553         fhMassSplitEFractionNLocMaxNEbin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
 554         outputContainer->Add(fhMassSplitEFractionNLocMaxNEbin[i][j]) ;
 555
 556         if(i>0) // skip first entry in array, general case not filled
 557         {
 558           fhMCGenFracNLocMaxEbin[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%d",pname[i].Data(),j),
 559                                                    Form("NLM vs E, %s, E bin %d",ptype[i].Data(),j),
 560                                                    200,0,2,nMaxBins,0,nMaxBins);
 561           fhMCGenFracNLocMaxEbin[i][j]->SetYTitle("NLM");
 562           fhMCGenFracNLocMaxEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
 563           outputContainer->Add(fhMCGenFracNLocMaxEbin[i][j]) ;
 564
 565           fhMCGenFracNLocMaxEbinMatched[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%dMatched",pname[i].Data(),j),
 566                                                           Form("NLM vs E, %s, E bin %d, matched to a track",ptype[i].Data(),j),
 567                                                           200,0,2,nMaxBins,0,nMaxBins);
 568           fhMCGenFracNLocMaxEbinMatched[i][j]->SetYTitle("NLM");
 569           fhMCGenFracNLocMaxEbinMatched[i][j]->SetXTitle("E_{gen} / E_{reco}");
 570           outputContainer->Add(fhMCGenFracNLocMaxEbinMatched[i][j]) ;
 571
 572           fhMassMCGenFracNLocMax1Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
 573                                                         Form("Invariant mass of 2 highest energy cells vs E, %s, E bin %d",ptype[i].Data(),j),
 574                                                         200,0,2,mbins,mmin,mmax);
 575           fhMassMCGenFracNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
 576           fhMassMCGenFracNLocMax1Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
 577           outputContainer->Add(fhMassMCGenFracNLocMax1Ebin[i][j]) ;
 578
 579           fhMassMCGenFracNLocMax2Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
 580                                                         Form("Invariant mass of 2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
 581                                                         200,0,2,mbins,mmin,mmax);
 582           fhMassMCGenFracNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
 583           fhMassMCGenFracNLocMax2Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
 584           outputContainer->Add(fhMassMCGenFracNLocMax2Ebin[i][j]) ;
 585
 586           fhMassMCGenFracNLocMaxNEbin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
 587                                                         Form("Invariant mass of N>2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
 588                                                         200,0,2,mbins,mmin,mmax);
 589           fhMassMCGenFracNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
 590           fhMassMCGenFracNLocMaxNEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
 591           outputContainer->Add(fhMassMCGenFracNLocMaxNEbin[i][j]) ;
 592
 593           fhM02MCGenFracNLocMax1Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
 594                                                           Form("#lambda_{0}^{2} vs E for N max  = 1 %s, E bin %d",ptype[i].Data(), j),
 595                                                           200,0,2,ssbins,ssmin,ssmax);
 596           fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 597           fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
 598           outputContainer->Add(fhM02MCGenFracNLocMax1Ebin[i][j]) ;
 599
 600           fhM02MCGenFracNLocMax2Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
 601                                                           Form("#lambda_{0}^{2} vs E for N max  = 2 %s, E bin %d",ptype[i].Data(),j),
 602                                                           200,0,2,ssbins,ssmin,ssmax);
 603           fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 604           fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
 605           outputContainer->Add(fhM02MCGenFracNLocMax2Ebin[i][j]) ;
 606
 607           fhM02MCGenFracNLocMaxNEbin[i][j]    = new TH2F(Form("hM02MCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
 608                                                          Form("#lambda_{0}^{2} vs E for N max  > 2 %s, E bin %d",ptype[i].Data(),j),
 609                                                          200,0,2,ssbins,ssmin,ssmax);
 610           fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
 611           fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
 612           outputContainer->Add(fhM02MCGenFracNLocMaxNEbin[i][j]) ;
 613         }
 614       }
 615   } // MC particle list
 616
 617   for(Int_t i = 0; i < 4; i++)
 618   {
 619
 620     if(IsDataMC())
 621     {
 622       fhMCAsymM02NLocMax1MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax1MCPi0Ebin%d",i),
 623                                                  Form("Asymmetry of MC #pi^{0} of 2 highest energy cells #lambda_{0}^{2}, E bin %d",i),
 624                                                  ssbins,ssmin,ssmax,100,0,1);
 625       fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
 626       fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
 627       outputContainer->Add(fhMCAsymM02NLocMax1MCPi0Ebin[i]) ;
 628
 629       fhMCAsymM02NLocMax2MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax2MCPi0Ebin%d",i),
 630                                                  Form("Asymmetry of MC #pi^{0} of 2 local maxima cells #lambda_{0}^{2}, E bin %d",i),
 631                                                  ssbins,ssmin,ssmax,100,0,1);
 632       fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
 633       fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
 634       outputContainer->Add(fhMCAsymM02NLocMax2MCPi0Ebin[i]) ;
 635
 636       fhMCAsymM02NLocMaxNMCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMaxNMCPi0Ebin%d",i),
 637                                                  Form("Asymmetry of MC #pi^{0} of N>2 local maxima cells vs #lambda_{0}^{2}, E bin %d",i),
 638                                                  ssbins,ssmin,ssmax,100,0,1);
 639       fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetYTitle("Decay asymmetry");
 640       fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
 641       outputContainer->Add(fhMCAsymM02NLocMaxNMCPi0Ebin[i]) ;
 642     }
 643
 644     fhMassM02NLocMax1Ebin[i]  = new TH2F(Form("hMassM02NLocMax1Ebin%d",i),
 645                                         Form("Invariant mass of 2 highest energy cells #lambda_{0}^{2}, E bin %d",i),
 646                                         ssbins,ssmin,ssmax,mbins,mmin,mmax);
 647     fhMassM02NLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
 648     fhMassM02NLocMax1Ebin[i]->SetXTitle("#lambda_{0}^{2}");
 649     outputContainer->Add(fhMassM02NLocMax1Ebin[i]) ;
 650
 651     fhMassM02NLocMax2Ebin[i]  = new TH2F(Form("hMassM02NLocMax2Ebin%d",i),
 652                                         Form("Invariant mass of 2 local maxima cells #lambda_{0}^{2}, E bin %d",i),
 653                                         ssbins,ssmin,ssmax,mbins,mmin,mmax);
 654     fhMassM02NLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
 655     fhMassM02NLocMax2Ebin[i]->SetXTitle("#lambda_{0}^{2}");
 656     outputContainer->Add(fhMassM02NLocMax2Ebin[i]) ;
 657
 658     fhMassM02NLocMaxNEbin[i]  = new TH2F(Form("hMassM02NLocMaxNEbin%d",i),
 659                                         Form("Invariant mass of N>2 local maxima cells vs #lambda_{0}^{2}, E bin %d",i),
 660                                         ssbins,ssmin,ssmax,mbins,mmin,mmax);
 661     fhMassM02NLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
 662     fhMassM02NLocMaxNEbin[i]->SetXTitle("#lambda_{0}^{2}");
 663     outputContainer->Add(fhMassM02NLocMaxNEbin[i]) ;
 664
 665     if(fFillSSExtraHisto)
 666     {
 667       fhMassDispEtaNLocMax1Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax1Ebin%d",i),
 668                                                Form("Invariant mass of 2 highest energy cells #sigma_{#eta #eta}^{2}, E bin %d",i),
 669                                                ssbins,ssmin,ssmax,mbins,mmin,mmax);
 670       fhMassDispEtaNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
 671       fhMassDispEtaNLocMax1Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
 672       outputContainer->Add(fhMassDispEtaNLocMax1Ebin[i]) ;
 673
 674       fhMassDispEtaNLocMax2Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax2Ebin%d",i),
 675                                                Form("Invariant mass of 2 local maxima cells #sigma_{#eta #eta}^{2}, E bin %d",i),
 676                                                ssbins,ssmin,ssmax,mbins,mmin,mmax);
 677       fhMassDispEtaNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
 678       fhMassDispEtaNLocMax2Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
 679       outputContainer->Add(fhMassDispEtaNLocMax2Ebin[i]) ;
 680
 681       fhMassDispEtaNLocMaxNEbin[i]  = new TH2F(Form("hMassDispEtaNLocMaxNEbin%d",i),
 682                                                Form("Invariant mass of N>2 local maxima cells vs #sigma_{#eta #eta}^{2}, E bin %d",i),
 683                                                ssbins,ssmin,ssmax,mbins,mmin,mmax);
 684       fhMassDispEtaNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
 685       fhMassDispEtaNLocMaxNEbin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
 686       outputContainer->Add(fhMassDispEtaNLocMaxNEbin[i]) ;
 687
 688       fhMassDispPhiNLocMax1Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax1Ebin%d",i),
 689                                                Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E bin %d",i),
 690                                                ssbins,ssmin,ssmax,mbins,mmin,mmax);
 691       fhMassDispPhiNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
 692       fhMassDispPhiNLocMax1Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
 693       outputContainer->Add(fhMassDispPhiNLocMax1Ebin[i]) ;
 694
 695       fhMassDispPhiNLocMax2Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax2Ebin%d",i),
 696                                                Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E bin %d",i),
 697                                                ssbins,ssmin,ssmax,mbins,mmin,mmax);
 698       fhMassDispPhiNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
 699       fhMassDispPhiNLocMax2Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
 700       outputContainer->Add(fhMassDispPhiNLocMax2Ebin[i]) ;
 701
 702       fhMassDispPhiNLocMaxNEbin[i]  = new TH2F(Form("hMassDispPhiNLocMaxNEbin%d",i),
 703                                                Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, E bin %d",i),
 704                                                ssbins,ssmin,ssmax,mbins,mmin,mmax);
 705       fhMassDispPhiNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
 706       fhMassDispPhiNLocMaxNEbin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
 707       outputContainer->Add(fhMassDispPhiNLocMaxNEbin[i]) ;
 708
 709       fhMassDispAsyNLocMax1Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax1Ebin%d",i),
 710                                                Form("Invariant mass of 2 highest energy cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E bin %d",i),
 711                                                200,-1,1,mbins,mmin,mmax);
 712       fhMassDispAsyNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
 713       fhMassDispAsyNLocMax1Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
 714       outputContainer->Add(fhMassDispAsyNLocMax1Ebin[i]) ;
 715
 716       fhMassDispAsyNLocMax2Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax2Ebin%d",i),
 717                                                Form("Invariant mass of 2 local maxima cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E bin %d",i),
 718                                                200,-1,1,mbins,mmin,mmax);
 719       fhMassDispAsyNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
 720       fhMassDispAsyNLocMax2Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
 721       outputContainer->Add(fhMassDispAsyNLocMax2Ebin[i]) ;
 722
 723       fhMassDispAsyNLocMaxNEbin[i]  = new TH2F(Form("hMassDispAsyNLocMaxNEbin%d",i),
 724                                                Form("Invariant mass of N>2 local maxima cells vs A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E bin %d",i),
 725                                                200,-1,1,mbins,mmin,mmax);
 726       fhMassDispAsyNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
 727       fhMassDispAsyNLocMaxNEbin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
 728       outputContainer->Add(fhMassDispAsyNLocMaxNEbin[i]) ;
 729     }
 730   }
 731
 732   if(fFillTMResidualHisto)
 733   {
 734     for(Int_t i = 0; i < n; i++)
 735     {
 736
 737       fhTrackMatchedDEtaLocMax1[i]  = new TH2F
 738       (Form("hTrackMatchedDEtaLocMax1%s",pname[i].Data()),
 739        Form("d#eta of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
 740        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
 741       fhTrackMatchedDEtaLocMax1[i]->SetYTitle("d#eta");
 742       fhTrackMatchedDEtaLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
 743
 744       fhTrackMatchedDPhiLocMax1[i]  = new TH2F
 745       (Form("hTrackMatchedDPhiLocMax1%s",pname[i].Data()),
 746        Form("d#phi of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
 747        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
 748       fhTrackMatchedDPhiLocMax1[i]->SetYTitle("d#phi (rad)");
 749       fhTrackMatchedDPhiLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
 750
 751       outputContainer->Add(fhTrackMatchedDEtaLocMax1[i]) ;
 752       outputContainer->Add(fhTrackMatchedDPhiLocMax1[i]) ;
 753
 754       fhTrackMatchedDEtaLocMax2[i]  = new TH2F
 755       (Form("hTrackMatchedDEtaLocMax2%s",pname[i].Data()),
 756        Form("d#eta of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
 757        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
 758       fhTrackMatchedDEtaLocMax2[i]->SetYTitle("d#eta");
 759       fhTrackMatchedDEtaLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
 760
 761       fhTrackMatchedDPhiLocMax2[i]  = new TH2F
 762       (Form("hTrackMatchedDPhiLocMax2%s",pname[i].Data()),
 763        Form("d#phi of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
 764        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
 765       fhTrackMatchedDPhiLocMax2[i]->SetYTitle("d#phi (rad)");
 766       fhTrackMatchedDPhiLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
 767
 768       outputContainer->Add(fhTrackMatchedDEtaLocMax2[i]) ;
 769       outputContainer->Add(fhTrackMatchedDPhiLocMax2[i]) ;
 770
 771       fhTrackMatchedDEtaLocMaxN[i]  = new TH2F
 772       (Form("hTrackMatchedDEtaLocMaxN%s",pname[i].Data()),
 773        Form("d#eta of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
 774        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
 775       fhTrackMatchedDEtaLocMaxN[i]->SetYTitle("d#eta");
 776       fhTrackMatchedDEtaLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
 777
 778       fhTrackMatchedDPhiLocMaxN[i]  = new TH2F
 779       (Form("hTrackMatchedDPhiLocMaxN%s",pname[i].Data()),
 780        Form("d#phi of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
 781        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
 782       fhTrackMatchedDPhiLocMaxN[i]->SetYTitle("d#phi (rad)");
 783       fhTrackMatchedDPhiLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
 784
 785       outputContainer->Add(fhTrackMatchedDEtaLocMaxN[i]) ;
 786       outputContainer->Add(fhTrackMatchedDPhiLocMaxN[i]) ;
 787     }
 788   }
 789
 790   if(fFillAngleHisto)
 791   {
 792     for(Int_t j = 0; j < 2; j++)
 793     {
 794
 795       fhAnglePairLocMax1[j]  = new TH2F(Form("hAnglePairLocMax1%s",sMatched[j].Data()),
 796                                         Form("Opening angle of 2 highest energy cells vs pair Energy, %s",sMatched[j].Data()),
 797                                         nptbins,ptmin,ptmax,200,0,0.2);
 798       fhAnglePairLocMax1[j]->SetYTitle("#alpha (rad)");
 799       fhAnglePairLocMax1[j]->SetXTitle("E (GeV)");
 800       outputContainer->Add(fhAnglePairLocMax1[j]) ;
 801
 802       fhAnglePairLocMax2[j]  = new TH2F(Form("hAnglePairLocMax2%s",sMatched[j].Data()),
 803                                         Form("Opening angle of 2 local maxima cells vs Energy, %s",sMatched[j].Data()),
 804                                         nptbins,ptmin,ptmax,200,0,0.2);
 805       fhAnglePairLocMax2[j]->SetYTitle("#alpha (rad)");
 806       fhAnglePairLocMax2[j]->SetXTitle("E (GeV)");
 807       outputContainer->Add(fhAnglePairLocMax2[j]) ;
 808
 809       fhAnglePairLocMaxN[j]  = new TH2F(Form("hAnglePairLocMaxN%s",sMatched[j].Data()),
 810                                         Form("Opening angle of N>2 local maxima cells vs Energy, %s",sMatched[j].Data()),
 811                                         nptbins,ptmin,ptmax,200,0,0.2);
 812       fhAnglePairLocMaxN[j]->SetYTitle("#alpha (rad)");
 813       fhAnglePairLocMaxN[j]->SetXTitle("E (GeV)");
 814       outputContainer->Add(fhAnglePairLocMaxN[j]) ;
 815
 816       fhAnglePairMassLocMax1[j]  = new TH2F(Form("hAnglePairMassLocMax1%s",sMatched[j].Data()),
 817                                             Form("Opening angle of 2 highest energy cells vs Mass for E > 7 GeV, %s",sMatched[j].Data()),
 818                                             mbins,mmin,mmax,200,0,0.2);
 819       fhAnglePairMassLocMax1[j]->SetXTitle("M (GeV/c^{2})");
 820       fhAnglePairMassLocMax1[j]->SetYTitle("#alpha (rad)");
 821       outputContainer->Add(fhAnglePairMassLocMax1[j]) ;
 822
 823       fhAnglePairMassLocMax2[j]  = new TH2F(Form("hAnglePairMassLocMax2%s",sMatched[j].Data()),
 824                                             Form("Opening angle of 2 local maxima cells vs Mass for E > 7 GeV, %s",sMatched[j].Data()),
 825                                             mbins,mmin,mmax,200,0,0.2);
 826       fhAnglePairMassLocMax2[j]->SetXTitle("M (GeV/c^{2})");
 827       fhAnglePairMassLocMax2[j]->SetYTitle("#alpha (rad)");
 828       outputContainer->Add(fhAnglePairMassLocMax2[j]) ;
 829
 830       fhAnglePairMassLocMaxN[j]  = new TH2F(Form("hAnglePairMassLocMaxN%s",sMatched[j].Data()),
 831                                             Form("Opening angle of N>2 local maxima cells vs Mass for E > 7 GeV, %s",sMatched[j].Data()),
 832                                             mbins,mmin,mmax,200,0,0.2);
 833       fhAnglePairMassLocMaxN[j]->SetXTitle("M (GeV/c^{2})");
 834       fhAnglePairMassLocMaxN[j]->SetYTitle("#alpha (rad)");
 835       outputContainer->Add(fhAnglePairMassLocMaxN[j]) ;
 836
 837     }
 838   }
 839
 840   return outputContainer ;
 841
 842 }
 843
 844 //___________________________________________
 845 void AliAnaInsideClusterInvariantMass::Init()
 846 {
 847   //Init
 848   //Do some checks
 849   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD())
 850   {
 851     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 852     abort();
 853   }
 854   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD())
 855   {
 856     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 857     abort();
 858   }
 859
 860   if( GetReader()->GetDataType() == AliCaloTrackReader::kMC )
 861   {
 862     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
 863     abort();
 864
 865   }
 866
 867 }
 868
 869 //_____________________________________________________
 870 void AliAnaInsideClusterInvariantMass::InitParameters()
 871 {
 872   //Initialize the parameters of the analysis.
 873   AddToHistogramsName("AnaPi0InsideClusterInvariantMass_");
 874
 875   fCalorimeter = "EMCAL" ;
 876
 877   fM02MinCut   = 0.26 ;
 878   fM02MaxCut   = 10 ;
 879
 880   fMinNCells   = 4 ;
 881   fMinBadDist  = 2 ;
 882
 883 }
 884
 885
 886 //__________________________________________________________________
 887 void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
 888 {
 889   //Search for pi0 in fCalorimeter with shower shape analysis
 890
 891   TObjArray * pl       = 0x0;
 892   AliVCaloCells* cells = 0x0;
 893
 894   //Select the Calorimeter of the photon
 895   if(fCalorimeter == "PHOS")
 896   {
 897     pl    = GetPHOSClusters();
 898     cells = GetPHOSCells();
 899   }
 900   else if (fCalorimeter == "EMCAL")
 901   {
 902     pl    = GetEMCALClusters();
 903     cells = GetEMCALCells();
 904   }
 905
 906   if(!pl || !cells)
 907   {
 908     Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
 909     return;
 910   }
 911
 912         if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
 913
 914   for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++)
 915   {
 916     AliVCluster * cluster = (AliVCluster*) (pl->At(icluster));
 917
 918     // Study clusters with large shape parameter
 919     Float_t en = cluster->E();
 920     Float_t l0 = cluster->GetM02();
 921     Int_t   nc = cluster->GetNCells();
 922     Float_t bd = cluster->GetDistanceToBadChannel() ;
 923
 924
 925     //If too small or big E or low number of cells, or close to a bad channel skip it
 926     if( en < GetMinEnergy() || en > GetMaxEnergy() || nc < fMinNCells || bd < fMinBadDist) continue ;
 927
 928     //printf("en %2.2f, GetMinEnergy() %2.2f, GetMaxEnergy() %2.2f, nc %d, fMinNCells %d,  bd %2.2f, fMinBadDist %2.2f\n",
 929     //       en,GetMinEnergy(), GetMaxEnergy(), nc, fMinNCells, bd, fMinBadDist);
 930
 931     // Get more Shower Shape parameters
 932     Float_t ll0  = 0., ll1  = 0.;
 933     Float_t disp= 0., dispEta = 0., dispPhi    = 0.;
 934     Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
 935
 936     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
 937                                                                                  ll0, ll1, disp, dispEta, dispPhi, sEta, sPhi, sEtaPhi);
 938
 939     Float_t dispAsy = -1;
 940     if(dispEta+dispPhi >0 ) dispAsy = (dispPhi-dispEta) / (dispPhi+dispEta);
 941
 942     Int_t    nMax = 0;
 943     Double_t mass = 0., angle = 0.;
 944     Double_t e1   = 0., e2    = 0.;
 945     Int_t pidTag = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(cluster,cells,GetCaloUtils(),
 946                                                                                GetVertex(0), nMax, mass, angle,e1,e2);
 947     if (nMax <= 0)
 948     {
 949       if(GetDebug() > 0 )
 950         printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found! It did not pass CaloPID selection criteria \n");
 951
 952       return;
 953     }
 954
 955     Float_t splitFrac = (e1+e2)/en;
 956     //printf("e1 %f, e2 %f, sum %f, cluster %f, fract %f \n",e1,e2,e1+e2,en,splitFrac);
 957
 958     Bool_t  matched   = IsTrackMatched(cluster,GetReader()->GetInputEvent());
 959
 960     fhNLocMax[0][matched]->Fill(en,nMax);
 961
 962     if     ( nMax == 1  ) { fhM02NLocMax1[0][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax1[0][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax1[0][matched]->Fill(en,nc) ; }
 963     else if( nMax == 2  ) { fhM02NLocMax2[0][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax2[0][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax2[0][matched]->Fill(en,nc) ; }
 964     else if( nMax >= 3  ) { fhM02NLocMaxN[0][matched]->Fill(en,l0) ; fhSplitEFractionNLocMaxN[0][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMaxN[0][matched]->Fill(en,nc) ; }
 965     else printf("N max smaller than 1 -> %d \n",nMax);
 966
 967     Float_t dZ  = cluster->GetTrackDz();
 968     Float_t dR  = cluster->GetTrackDx();
 969
 970     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
 971     {
 972       dR = 2000., dZ = 2000.;
 973       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
 974     }
 975     //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
 976
 977     if(TMath::Abs(dR) < 999 && fFillTMResidualHisto)
 978     {
 979       if     ( nMax == 1  ) { fhTrackMatchedDEtaLocMax1[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMax1[0]->Fill(en,dR); }
 980       else if( nMax == 2  ) { fhTrackMatchedDEtaLocMax2[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMax2[0]->Fill(en,dR); }
 981       else if( nMax >= 3  ) { fhTrackMatchedDEtaLocMaxN[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMaxN[0]->Fill(en,dR); }
 982     }
 983
 984     // Play with the MC stack if available
 985     // Check origin of the candidates
 986     Int_t mcindex   = -1;
 987     Float_t eprim   =  0;
 988     Float_t asymGen = -2;
 989     Int_t mcLabel   = cluster->GetLabel();
 990     if(IsDataMC())
 991     {
 992       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader(), 0);
 993
 994       if      ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )      mcindex = kmcPi0;
 995       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )      mcindex = kmcEta;
 996       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
 997                !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPhoton;
 998       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
 999                 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
1000       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))   mcindex = kmcElectron;
1001       else                                                                                mcindex = kmcHadron;
1002
1003       fhNLocMax[mcindex][matched]->Fill(en,nMax);
1004
1005       if     (nMax == 1 ) { fhM02NLocMax1[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax1[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax1[mcindex][matched]->Fill(en,nc) ; }
1006       else if(nMax == 2 ) { fhM02NLocMax2[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax2[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax2[mcindex][matched]->Fill(en,nc) ; }
1007       else if(nMax >= 3 ) { fhM02NLocMaxN[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMaxN[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMaxN[mcindex][matched]->Fill(en,nc) ; }
1008
1009       if(TMath::Abs(dR) < 999 && fFillTMResidualHisto)
1010       {
1011         if     ( nMax == 1  ) { fhTrackMatchedDEtaLocMax1[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMax1[mcindex]->Fill(en,dR); }
1012         else if( nMax == 2  ) { fhTrackMatchedDEtaLocMax2[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMax2[mcindex]->Fill(en,dR); }
1013         else if( nMax >= 3  ) { fhTrackMatchedDEtaLocMaxN[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMaxN[mcindex]->Fill(en,dR); }
1014       }
1015
1016       Bool_t ok = kFALSE;
1017       TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
1018       eprim = primary.E();
1019
1020       if(mcindex == kmcPi0 || mcindex == kmcEta)
1021       {
1022         if(mcindex == kmcPi0)
1023         {
1024           asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,111,GetReader(),ok));
1025           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok);
1026           if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
1027         }
1028         else
1029         {
1030           asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,221,GetReader(),ok));
1031           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok);
1032           if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
1033         }
1034       }
1035     }
1036
1037     Float_t efrac = eprim/en;
1038
1039     Int_t ebin = -1;
1040     if(en > 8  && en <= 12) ebin = 0;
1041     if(en > 12 && en <= 16) ebin = 1;
1042     if(en > 16 && en <= 20) ebin = 2;
1043     if(en > 20)             ebin = 3;
1044
1045     if(ebin >= 0 && IsDataMC())
1046     {
1047       if( !matched ) fhMCGenFracNLocMaxEbin       [mcindex][ebin]->Fill(efrac,nMax);
1048       else           fhMCGenFracNLocMaxEbinMatched[mcindex][ebin]->Fill(efrac,nMax);
1049     }
1050
1051     if     (nMax==1)
1052     {
1053       if( en > 7 )
1054       {
1055         fhMassM02NLocMax1    [0][matched]->Fill(l0     ,  mass );
1056         if(fFillSSExtraHisto)
1057         {
1058           fhMassDispEtaNLocMax1[0][matched]->Fill(dispEta,  mass );
1059           fhMassDispPhiNLocMax1[0][matched]->Fill(dispPhi,  mass );
1060           fhMassDispAsyNLocMax1[0][matched]->Fill(dispAsy,  mass );
1061         }
1062
1063         if(IsDataMC())
1064         {
1065           fhMassM02NLocMax1          [mcindex][matched]->Fill(l0     ,  mass  );
1066           fhMCGenFracNLocMax1        [mcindex][matched]->Fill(en     ,  efrac );
1067           if(!matched && ebin >= 0)
1068           {
1069             fhM02MCGenFracNLocMax1Ebin [mcindex][ebin]->Fill(efrac  ,  l0    );
1070             fhMassMCGenFracNLocMax1Ebin[mcindex][ebin]->Fill(efrac  ,  mass  );
1071             fhMCAsymM02NLocMax1MCPi0Ebin        [ebin]->Fill(l0     ,  asymGen );
1072           }
1073
1074           if(fFillSSExtraHisto)
1075           {
1076             fhMassDispEtaNLocMax1[mcindex][matched]->Fill(dispEta,  mass );
1077             fhMassDispPhiNLocMax1[mcindex][matched]->Fill(dispPhi,  mass );
1078             fhMassDispAsyNLocMax1[mcindex][matched]->Fill(dispAsy,  mass );
1079           }
1080         }
1081       }
1082
1083       if(!matched && ebin >= 0)
1084       {
1085         fhMassSplitEFractionNLocMax1Ebin[0][ebin]->Fill(splitFrac,  mass);
1086         if(IsDataMC())fhMassSplitEFractionNLocMax1Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
1087
1088         fhMassM02NLocMax1Ebin    [ebin]->Fill(l0     ,  mass    );
1089         if(fFillSSExtraHisto)
1090         {
1091           fhMassDispEtaNLocMax1Ebin[ebin]->Fill(dispEta,  mass );
1092           fhMassDispPhiNLocMax1Ebin[ebin]->Fill(dispPhi,  mass );
1093           fhMassDispAsyNLocMax1Ebin[ebin]->Fill(dispAsy,  mass );
1094         }
1095       }
1096     }
1097     else if(nMax==2)
1098     {
1099       if( en > 7 )
1100       {
1101         fhMassM02NLocMax2    [0][matched]->Fill(l0     ,  mass );
1102         if(fFillSSExtraHisto)
1103         {
1104           fhMassDispEtaNLocMax2[0][matched]->Fill(dispEta,  mass );
1105           fhMassDispPhiNLocMax2[0][matched]->Fill(dispPhi,  mass );
1106           fhMassDispAsyNLocMax2[0][matched]->Fill(dispAsy,  mass );
1107         }
1108
1109         if(IsDataMC())
1110         {
1111           fhMassM02NLocMax2    [mcindex][matched]->Fill(l0     ,  mass );
1112           fhMCGenFracNLocMax2  [mcindex][matched]->Fill(en     ,  efrac );
1113           if(!matched && ebin >= 0)
1114           {
1115             fhM02MCGenFracNLocMax2Ebin [mcindex][ebin]->Fill(efrac  ,  l0    );
1116             fhMassMCGenFracNLocMax2Ebin[mcindex][ebin]->Fill(efrac  ,  mass  );
1117             fhMCAsymM02NLocMax2MCPi0Ebin        [ebin]->Fill(l0     ,  asymGen );
1118
1119           }
1120           if(fFillSSExtraHisto)
1121           {
1122             fhMassDispEtaNLocMax2[mcindex][matched]->Fill(dispEta,  mass );
1123             fhMassDispPhiNLocMax2[mcindex][matched]->Fill(dispPhi,  mass );
1124             fhMassDispAsyNLocMax2[mcindex][matched]->Fill(dispAsy,  mass );
1125           }
1126         }
1127       }
1128
1129       if(!matched && ebin >= 0)
1130       {
1131         fhMassSplitEFractionNLocMax2Ebin[0][ebin]->Fill(splitFrac,  mass);
1132         if(IsDataMC())fhMassSplitEFractionNLocMax2Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
1133
1134         fhMassM02NLocMax2Ebin    [ebin]->Fill(l0     ,  mass );
1135         if(fFillSSExtraHisto)
1136         {
1137           fhMassDispEtaNLocMax2Ebin[ebin]->Fill(dispEta,  mass );
1138           fhMassDispPhiNLocMax2Ebin[ebin]->Fill(dispPhi,  mass );
1139           fhMassDispAsyNLocMax2Ebin[ebin]->Fill(dispAsy,  mass );
1140         }
1141       }
1142     }
1143     else if(nMax > 2 )
1144     {
1145       if( en > 7 )
1146       {
1147         fhMassM02NLocMaxN    [0][matched]->Fill(l0     ,  mass );
1148         if(fFillSSExtraHisto)
1149         {
1150           fhMassDispEtaNLocMaxN[0][matched]->Fill(dispEta,  mass );
1151           fhMassDispPhiNLocMaxN[0][matched]->Fill(dispPhi,  mass );
1152           fhMassDispAsyNLocMaxN[0][matched]->Fill(dispAsy,  mass );
1153         }
1154
1155         if(IsDataMC())
1156         {
1157           fhMassM02NLocMaxN    [mcindex][matched]->Fill(l0     ,  mass );
1158           fhMCGenFracNLocMaxN  [mcindex][matched]->Fill(en     ,  efrac );
1159           if(!matched && ebin >= 0)
1160           {
1161             fhM02MCGenFracNLocMaxNEbin [mcindex][ebin]->Fill(efrac  ,  l0     );
1162             fhMassMCGenFracNLocMaxNEbin[mcindex][ebin]->Fill(efrac  ,  mass   );
1163             fhMCAsymM02NLocMaxNMCPi0Ebin        [ebin]->Fill(l0     ,  asymGen);
1164           }
1165           if(fFillSSExtraHisto)
1166           {
1167             fhMassDispEtaNLocMaxN[mcindex][matched]->Fill(dispEta,  mass );
1168             fhMassDispPhiNLocMaxN[mcindex][matched]->Fill(dispPhi,  mass );
1169             fhMassDispAsyNLocMaxN[mcindex][matched]->Fill(dispAsy,  mass );
1170           }
1171         }
1172       }
1173
1174       if(!matched && ebin >= 0)
1175       {
1176         fhMassSplitEFractionNLocMaxNEbin[0][ebin]->Fill(splitFrac,  mass);
1177         if(IsDataMC())fhMassSplitEFractionNLocMaxNEbin[mcindex][ebin]->Fill(splitFrac,  mass);
1178
1179         fhMassM02NLocMaxNEbin    [ebin]->Fill(l0     ,  mass );
1180         if(fFillSSExtraHisto)
1181         {
1182           fhMassDispEtaNLocMaxNEbin[ebin]->Fill(dispEta,  mass );
1183           fhMassDispPhiNLocMaxNEbin[ebin]->Fill(dispPhi,  mass );
1184           fhMassDispAsyNLocMaxNEbin[ebin]->Fill(dispAsy,  mass );
1185         }
1186       }
1187     }
1188
1189     //---------------------------------------------------------------------
1190     // From here only if M02 is large but not too large, fill histograms
1191     //---------------------------------------------------------------------
1192
1193     if( l0 < fM02MinCut || l0 > fM02MaxCut ) continue ;
1194
1195     fhNLocMaxM02Cut[0][matched]->Fill(en,nMax);
1196     if(IsDataMC()) fhNLocMaxM02Cut[mcindex][matched]->Fill(en,nMax);
1197
1198     if     (nMax==1)
1199     {
1200       fhMassNLocMax1[0][matched] ->Fill(en,mass );
1201
1202       if(fFillAngleHisto)
1203       {
1204         fhAnglePairLocMax1[matched]->Fill(en,angle);
1205       if( en > 7 )
1206         fhAnglePairMassLocMax1[matched]->Fill(mass,angle);
1207       }
1208
1209       if     (pidTag==AliCaloPID::kPhoton) fhM02ConLocMax1[0][matched]->Fill(en,l0);
1210       else if(pidTag==AliCaloPID::kPi0   ) fhM02Pi0LocMax1[0][matched]->Fill(en,l0);
1211       else if(pidTag==AliCaloPID::kEta)    fhM02EtaLocMax1[0][matched]->Fill(en,l0);
1212     }
1213     else if(nMax==2)
1214     {
1215       fhMassNLocMax2[0] [matched]->Fill(en,mass );
1216
1217       if(fFillAngleHisto)
1218       {
1219         fhAnglePairLocMax2[matched]->Fill(en,angle);
1220         if( en > 7 )
1221           fhAnglePairMassLocMax2[matched]->Fill(mass,angle);
1222       }
1223
1224       if     (pidTag==AliCaloPID::kPhoton) fhM02ConLocMax2[0][matched]->Fill(en,l0);
1225       else if(pidTag==AliCaloPID::kPi0   ) fhM02Pi0LocMax2[0][matched]->Fill(en,l0);
1226       else if(pidTag==AliCaloPID::kEta)    fhM02EtaLocMax2[0][matched]->Fill(en,l0);
1227     }
1228     else if(nMax >2)
1229     {
1230       fhMassNLocMaxN[0] [matched]->Fill(en,mass );
1231
1232       if(fFillAngleHisto)
1233       {
1234         fhAnglePairLocMaxN[matched]->Fill(en,angle);
1235         if( en > 7 )
1236           fhAnglePairMassLocMaxN[matched]->Fill(mass,angle);
1237       }
1238
1239       if     (pidTag==AliCaloPID::kPhoton) fhM02ConLocMaxN[0][matched]->Fill(en,l0);
1240       else if(pidTag==AliCaloPID::kPi0   ) fhM02Pi0LocMaxN[0][matched]->Fill(en,l0);
1241       else if(pidTag==AliCaloPID::kEta   ) fhM02EtaLocMaxN[0][matched]->Fill(en,l0);
1242     }
1243
1244
1245     if(IsDataMC())
1246     {
1247       if     (nMax==1)
1248       {
1249         fhMassNLocMax1[mcindex][matched]->Fill(en,mass);
1250         if     (pidTag==AliCaloPID::kPhoton) fhM02ConLocMax1[mcindex][matched]->Fill(en,l0);
1251         else if(pidTag==AliCaloPID::kPi0   ) fhM02Pi0LocMax1[mcindex][matched]->Fill(en,l0);
1252         else if(pidTag==AliCaloPID::kEta   ) fhM02EtaLocMax1[mcindex][matched]->Fill(en,l0);
1253       }
1254       else if(nMax==2)
1255       {
1256         fhMassNLocMax2[mcindex][matched]->Fill(en,mass);
1257         if     (pidTag==AliCaloPID::kPhoton) fhM02ConLocMax2[mcindex][matched]->Fill(en,l0);
1258         else if(pidTag==AliCaloPID::kPi0   ) fhM02Pi0LocMax2[mcindex][matched]->Fill(en,l0);
1259         else if(pidTag==AliCaloPID::kEta   ) fhM02EtaLocMax2[mcindex][matched]->Fill(en,l0);
1260       }
1261       else if(nMax >2)
1262       {
1263         fhMassNLocMaxN[mcindex][matched]->Fill(en,mass);
1264         if     (pidTag==AliCaloPID::kPhoton) fhM02ConLocMaxN[mcindex][matched]->Fill(en,l0);
1265         else if(pidTag==AliCaloPID::kPi0   ) fhM02Pi0LocMaxN[mcindex][matched]->Fill(en,l0);
1266         else if(pidTag==AliCaloPID::kEta   ) fhM02EtaLocMaxN[mcindex][matched]->Fill(en,l0);
1267       }
1268
1269     }//Work with MC truth first
1270
1271   }//loop
1272
1273   if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");
1274
1275 }
1276
1277 //______________________________________________________________________
1278 void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
1279 {
1280   //Print some relevant parameters set for the analysis
1281   if(! opt)
1282     return;
1283
1284   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1285   AliAnaCaloTrackCorrBaseClass::Print("");
1286   printf("Calorimeter     =     %s\n",  fCalorimeter.Data()) ;
1287   printf("Loc. Max. E > %2.2f\n",       GetCaloUtils()->GetLocalMaximaCutE());
1288   printf("Loc. Max. E Diff > %2.2f\n",  GetCaloUtils()->GetLocalMaximaCutEDiff());
1289   printf("Min. N Cells =%d \n",         fMinNCells) ;
1290   printf("Min. Dist. to Bad =%1.1f \n", fMinBadDist) ;
1291   printf("%2.2f < lambda_0^2 <%2.2f \n",fM02MinCut,fM02MaxCut);
1292
1293   printf("    \n") ;
1294
1295 }
1296
1297
1298