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