PWGGA/CaloTrackCorrelations/AliAnaInsideClusterInvariantMass.cxx

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