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 #include <TDatabasePDG.h>
  34
  35 // --- Analysis system ---
  36 #include "AliAnaInsideClusterInvariantMass.h"
  37 #include "AliCaloTrackReader.h"
  38 #include "AliMCAnalysisUtils.h"
  39 #include "AliStack.h"
  40 #include "AliFiducialCut.h"
  41 #include "TParticle.h"
  42 #include "AliVCluster.h"
  43 #include "AliAODEvent.h"
  44 #include "AliAODMCParticle.h"
  45 #include "AliEMCALGeoParams.h"
  46
  47 // --- Detectors ---
  48 //#include "AliPHOSGeoUtils.h"
  49 #include "AliEMCALGeometry.h"
  50
  51 ClassImp(AliAnaInsideClusterInvariantMass)
  52
  53 //__________________________________________________________________
  54 AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
  55   AliAnaCaloTrackCorrBaseClass(),
  56   fCalorimeter(""),
  57   fM02MaxCut(0),    fM02MinCut(0),
  58   fMinNCells(0),    fMinBadDist(0),
  59   fHistoECut(0),
  60   fFillAngleHisto(kFALSE),
  61   fFillTMHisto(kFALSE),
  62   fFillTMResidualHisto(kFALSE),
  63   fFillSSExtraHisto(kFALSE),
  64   fFillMCHisto(kFALSE),
  65   fFillSSWeightHisto(kFALSE),
  66   fFillEbinHisto(0),
  67   fFillMCOverlapHisto(0),
  68   fSSWeightN(0),              fSSECellCutN(0),
  69   fWSimu(0),
  70   fhMassM02CutNLocMax1(0),    fhMassM02CutNLocMax2(0),    fhMassM02CutNLocMaxN(0),
  71   fhAsymM02CutNLocMax1(0),    fhAsymM02CutNLocMax2(0),    fhAsymM02CutNLocMaxN(0),
  72   fhMassSplitECutNLocMax1(0), fhMassSplitECutNLocMax2(0), fhMassSplitECutNLocMaxN(0),
  73   fhMCGenFracAfterCutsNLocMax1MCPi0(0),
  74   fhMCGenFracAfterCutsNLocMax2MCPi0(0),
  75   fhMCGenFracAfterCutsNLocMaxNMCPi0(0),
  76   fhMCGenSplitEFracAfterCutsNLocMax1MCPi0(0),
  77   fhMCGenSplitEFracAfterCutsNLocMax2MCPi0(0),
  78   fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0(0),
  79   fhEventPlanePi0NLocMax1(0), fhEventPlaneEtaNLocMax1(0),
  80   fhEventPlanePi0NLocMax2(0), fhEventPlaneEtaNLocMax2(0),
  81   fhEventPlanePi0NLocMaxN(0), fhEventPlaneEtaNLocMaxN(0),
  82   fhClusterEtaPhiNLocMax1(0), fhClusterEtaPhiNLocMax2(0),  fhClusterEtaPhiNLocMaxN(0),
  83   fhPi0EtaPhiNLocMax1(0),     fhPi0EtaPhiNLocMax2(0),      fhPi0EtaPhiNLocMaxN(0),
  84   fhEtaEtaPhiNLocMax1(0),     fhEtaEtaPhiNLocMax2(0),      fhEtaEtaPhiNLocMaxN(0),
  85   fhPi0EPairDiffTimeNLM1(0),  fhPi0EPairDiffTimeNLM2(0),   fhPi0EPairDiffTimeNLMN(0),
  86   fhEtaEPairDiffTimeNLM1(0),  fhEtaEPairDiffTimeNLM2(0),   fhEtaEPairDiffTimeNLMN(0),
  87   fhMCPi0HighNLMPair(0),      fhMCPi0LowNLMPair(0),
  88   fhMCPi0AnyNLMPair(0),       fhMCPi0NoneNLMPair(0),
  89   fhMCPi0HighNLMPairNoMCMatch(0), fhMCPi0LowNLMPairNoMCMatch(0),
  90   fhMCPi0AnyNLMPairNoMCMatch(0),  fhMCPi0NoneNLMPairNoMCMatch(0),
  91   fhMCEOverlapType(0),            fhMCEOverlapTypeMatch(0)
  92 {
  93   //default ctor
  94
  95   // Init array of histograms
  96   for(Int_t i = 0; i < 8; i++)
  97   {
  98     for(Int_t j = 0; j < 2; j++)
  99     {
 100       fhMassNLocMax1[i][j]  = 0;
 101       fhMassNLocMax2[i][j]  = 0;
 102       fhMassNLocMaxN[i][j]  = 0;
 103       fhNLocMax[i][j]       = 0;
 104       fhNLocMaxM02Cut[i][j] = 0;
 105       fhM02NLocMax1[i][j]   = 0;
 106       fhM02NLocMax2[i][j]   = 0;
 107       fhM02NLocMaxN[i][j]   = 0;
 108       fhNCellNLocMax1[i][j] = 0;
 109       fhNCellNLocMax2[i][j] = 0;
 110       fhNCellNLocMaxN[i][j] = 0;
 111       fhM02Pi0NLocMax1[i][j] = 0;
 112       fhM02EtaNLocMax1[i][j] = 0;
 113       fhM02ConNLocMax1[i][j] = 0;
 114       fhM02Pi0NLocMax2[i][j] = 0;
 115       fhM02EtaNLocMax2[i][j] = 0;
 116       fhM02ConNLocMax2[i][j] = 0;
 117       fhM02Pi0NLocMaxN[i][j] = 0;
 118       fhM02EtaNLocMaxN[i][j] = 0;
 119       fhM02ConNLocMaxN[i][j] = 0;
 120
 121       fhMassPi0NLocMax1[i][j] = 0;
 122       fhMassEtaNLocMax1[i][j] = 0;
 123       fhMassConNLocMax1[i][j] = 0;
 124       fhMassPi0NLocMax2[i][j] = 0;
 125       fhMassEtaNLocMax2[i][j] = 0;
 126       fhMassConNLocMax2[i][j] = 0;
 127       fhMassPi0NLocMaxN[i][j] = 0;
 128       fhMassEtaNLocMaxN[i][j] = 0;
 129       fhMassConNLocMaxN[i][j] = 0;
 130
 131
 132       fhAsyPi0NLocMax1[i][j] = 0;
 133       fhAsyEtaNLocMax1[i][j] = 0;
 134       fhAsyConNLocMax1[i][j] = 0;
 135       fhAsyPi0NLocMax2[i][j] = 0;
 136       fhAsyEtaNLocMax2[i][j] = 0;
 137       fhAsyConNLocMax2[i][j] = 0;
 138       fhAsyPi0NLocMaxN[i][j] = 0;
 139       fhAsyEtaNLocMaxN[i][j] = 0;
 140       fhAsyConNLocMaxN[i][j] = 0;
 141
 142       fhMassM02NLocMax1[i][j]= 0;
 143       fhMassM02NLocMax2[i][j]= 0;
 144       fhMassM02NLocMaxN[i][j]= 0;
 145       fhMassDispEtaNLocMax1[i][j]= 0;
 146       fhMassDispEtaNLocMax2[i][j]= 0;
 147       fhMassDispEtaNLocMaxN[i][j]= 0;
 148       fhMassDispPhiNLocMax1[i][j]= 0;
 149       fhMassDispPhiNLocMax2[i][j]= 0;
 150       fhMassDispPhiNLocMaxN[i][j]= 0;
 151       fhMassDispAsyNLocMax1[i][j]= 0;
 152       fhMassDispAsyNLocMax2[i][j]= 0;
 153       fhMassDispAsyNLocMaxN[i][j]= 0;
 154
 155       fhSplitEFractionNLocMax1[i][j]=0;
 156       fhSplitEFractionNLocMax2[i][j]=0;
 157       fhSplitEFractionNLocMaxN[i][j]=0;
 158
 159       fhMCGenFracNLocMax1[i][j]= 0;
 160       fhMCGenFracNLocMax2[i][j]= 0;
 161       fhMCGenFracNLocMaxN[i][j]= 0;
 162
 163       fhMCGenSplitEFracNLocMax1[i][j]= 0;
 164       fhMCGenSplitEFracNLocMax2[i][j]= 0;
 165       fhMCGenSplitEFracNLocMaxN[i][j]= 0;
 166
 167       fhMCGenEFracvsSplitEFracNLocMax1[i][j]= 0;
 168       fhMCGenEFracvsSplitEFracNLocMax2[i][j]= 0;
 169       fhMCGenEFracvsSplitEFracNLocMaxN[i][j]= 0;
 170
 171       fhMCGenEvsSplitENLocMax1[i][j]= 0;
 172       fhMCGenEvsSplitENLocMax2[i][j]= 0;
 173       fhMCGenEvsSplitENLocMaxN[i][j]= 0;
 174
 175       fhAsymNLocMax1 [i][j] = 0;
 176       fhAsymNLocMax2 [i][j] = 0;
 177       fhAsymNLocMaxN [i][j] = 0;
 178
 179       fhMassAfterCutsNLocMax1[i][j] = 0;
 180       fhMassAfterCutsNLocMax2[i][j] = 0;
 181       fhMassAfterCutsNLocMaxN[i][j] = 0;
 182
 183       fhSplitEFractionAfterCutsNLocMax1[i][j] = 0 ;
 184       fhSplitEFractionAfterCutsNLocMax2[i][j] = 0 ;
 185       fhSplitEFractionAfterCutsNLocMaxN[i][j] = 0 ;
 186
 187       fhCentralityPi0NLocMax1[i][j] = 0 ;
 188       fhCentralityEtaNLocMax1[i][j] = 0 ;
 189
 190       fhCentralityPi0NLocMax2[i][j] = 0 ;
 191       fhCentralityEtaNLocMax2[i][j] = 0 ;
 192
 193       fhCentralityPi0NLocMaxN[i][j] = 0 ;
 194       fhCentralityEtaNLocMaxN[i][j] = 0 ;
 195     }
 196
 197     for(Int_t jj = 0; jj < 4; jj++)
 198     {
 199       fhM02MCGenFracNLocMax1Ebin[i][jj] = 0;
 200       fhM02MCGenFracNLocMax2Ebin[i][jj] = 0;
 201       fhM02MCGenFracNLocMaxNEbin[i][jj] = 0;
 202
 203       fhMassMCGenFracNLocMax1Ebin[i][jj]= 0;
 204       fhMassMCGenFracNLocMax2Ebin[i][jj]= 0;
 205       fhMassMCGenFracNLocMaxNEbin[i][jj]= 0;
 206
 207       fhMCGenFracNLocMaxEbin[i][jj]       = 0;
 208       fhMCGenFracNLocMaxEbinMatched[i][jj]= 0;
 209
 210       fhMassSplitEFractionNLocMax1Ebin[i][jj] = 0;
 211       fhMassSplitEFractionNLocMax2Ebin[i][jj] = 0;
 212       fhMassSplitEFractionNLocMaxNEbin[i][jj] = 0;
 213     }
 214
 215     fhTrackMatchedDEtaNLocMax1[i] = 0;
 216     fhTrackMatchedDPhiNLocMax1[i] = 0;
 217     fhTrackMatchedDEtaNLocMax2[i] = 0;
 218     fhTrackMatchedDPhiNLocMax2[i] = 0;
 219     fhTrackMatchedDEtaNLocMaxN[i] = 0;
 220     fhTrackMatchedDPhiNLocMaxN[i] = 0;
 221
 222     fhTrackMatchedDEtaNLocMax1Pos[i] = 0;
 223     fhTrackMatchedDPhiNLocMax1Pos[i] = 0;
 224     fhTrackMatchedDEtaNLocMax2Pos[i] = 0;
 225     fhTrackMatchedDPhiNLocMax2Pos[i] = 0;
 226     fhTrackMatchedDEtaNLocMaxNPos[i] = 0;
 227     fhTrackMatchedDPhiNLocMaxNPos[i] = 0;
 228
 229     fhTrackMatchedDEtaNLocMax1Neg[i] = 0;
 230     fhTrackMatchedDPhiNLocMax1Neg[i] = 0;
 231     fhTrackMatchedDEtaNLocMax2Neg[i] = 0;
 232     fhTrackMatchedDPhiNLocMax2Neg[i] = 0;
 233     fhTrackMatchedDEtaNLocMaxNNeg[i] = 0;
 234     fhTrackMatchedDPhiNLocMaxNNeg[i] = 0;
 235
 236     for(Int_t nlm = 0; nlm < 3; nlm++)
 237     {
 238       fhMCEM02Overlap0     [nlm][i] = 0;
 239       fhMCEM02Overlap1     [nlm][i] = 0;
 240       fhMCEM02OverlapN     [nlm][i] = 0;
 241       fhMCEM02Overlap0Match[nlm][i] = 0;
 242       fhMCEM02Overlap1Match[nlm][i] = 0;
 243       fhMCEM02OverlapNMatch[nlm][i] = 0;
 244
 245       fhMCEMassOverlap0     [nlm][i] = 0;
 246       fhMCEMassOverlap1     [nlm][i] = 0;
 247       fhMCEMassOverlapN     [nlm][i] = 0;
 248       fhMCEMassOverlap0Match[nlm][i] = 0;
 249       fhMCEMassOverlap1Match[nlm][i] = 0;
 250       fhMCEMassOverlapNMatch[nlm][i] = 0;
 251
 252       fhMCENOverlaps       [nlm][i] = 0;
 253       fhMCENOverlapsMatch  [nlm][i] = 0;
 254
 255       if(i > 3) continue ;
 256
 257       fhMCPi0MassM02Overlap0     [nlm][i] = 0;
 258       fhMCPi0MassM02Overlap1     [nlm][i] = 0;
 259       fhMCPi0MassM02OverlapN     [nlm][i] = 0;
 260       fhMCPi0MassM02Overlap0Match[nlm][i] = 0;
 261       fhMCPi0MassM02Overlap1Match[nlm][i] = 0;
 262       fhMCPi0MassM02OverlapNMatch[nlm][i] = 0;
 263     }
 264   }
 265
 266   for(Int_t i = 0; i < 2; i++)
 267   {
 268     fhAnglePairNLocMax1    [i] = 0;
 269     fhAnglePairNLocMax2    [i] = 0;
 270     fhAnglePairNLocMaxN    [i] = 0;
 271     fhAnglePairMassNLocMax1[i] = 0;
 272     fhAnglePairMassNLocMax2[i] = 0;
 273     fhAnglePairMassNLocMaxN[i] = 0;
 274     fhSplitEFractionvsAsyNLocMax1[i] = 0;
 275     fhSplitEFractionvsAsyNLocMax2[i] = 0;
 276     fhSplitEFractionvsAsyNLocMaxN[i] = 0;
 277   }
 278
 279   for(Int_t i = 0; i < 4; i++)
 280   {
 281     fhMassM02NLocMax1Ebin[i] = 0 ;
 282     fhMassM02NLocMax2Ebin[i] = 0 ;
 283     fhMassM02NLocMaxNEbin[i] = 0 ;
 284
 285     fhMassAsyNLocMax1Ebin[i] = 0 ;
 286     fhMassAsyNLocMax2Ebin[i] = 0 ;
 287     fhMassAsyNLocMaxNEbin[i] = 0 ;
 288
 289     fhAsyMCGenRecoNLocMax1EbinPi0[i] = 0 ;
 290     fhAsyMCGenRecoNLocMax2EbinPi0[i] = 0 ;
 291     fhAsyMCGenRecoNLocMaxNEbinPi0[i] = 0 ;
 292
 293     fhMassDispEtaNLocMax1Ebin[i] = 0 ;
 294     fhMassDispEtaNLocMax2Ebin[i] = 0 ;
 295     fhMassDispEtaNLocMaxNEbin[i] = 0 ;
 296
 297     fhMassDispPhiNLocMax1Ebin[i] = 0 ;
 298     fhMassDispPhiNLocMax2Ebin[i] = 0 ;
 299     fhMassDispPhiNLocMaxNEbin[i] = 0 ;
 300
 301     fhMassDispAsyNLocMax1Ebin[i] = 0 ;
 302     fhMassDispAsyNLocMax2Ebin[i] = 0 ;
 303     fhMassDispAsyNLocMaxNEbin[i] = 0 ;
 304
 305     fhMCAsymM02NLocMax1MCPi0Ebin[i] = 0 ;
 306     fhMCAsymM02NLocMax2MCPi0Ebin[i] = 0 ;
 307     fhMCAsymM02NLocMaxNMCPi0Ebin[i] = 0 ;
 308   }
 309
 310   for(Int_t nlm = 0; nlm < 3; nlm++)
 311   {
 312     fhPi0CellE       [nlm] = 0 ;
 313     fhPi0CellEFrac   [nlm] = 0 ;
 314     fhPi0CellLogEFrac[nlm] = 0 ;
 315
 316     fhPi0CellEMaxEMax2Frac   [nlm] = 0 ;
 317     fhPi0CellEMaxClusterFrac [nlm] = 0 ;
 318     fhPi0CellEMax2ClusterFrac[nlm] = 0 ;
 319
 320     fhPi0CellEMaxFrac [nlm] = 0 ;
 321     fhPi0CellEMax2Frac[nlm] = 0 ;
 322
 323     for(Int_t i = 0; i < 10; i++)
 324     {
 325       fhM02WeightPi0  [nlm][i] = 0;
 326       fhM02ECellCutPi0[nlm][i] = 0;
 327     }
 328   }
 329
 330   InitParameters();
 331
 332 }
 333
 334 //_______________________________________________________________________________________________________
 335 void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* cluster, const Int_t mcindex)
 336 {
 337   // Check origin NLM tower of the cluster, when MC gives merged pi0
 338
 339   if(!fFillMCOverlapHisto) return;
 340
 341   if(!IsDataMC()) return;
 342
 343   if(mcindex != kmcPi0 && mcindex != kmcPi0Conv) return;
 344
 345   const UInt_t nc = cluster->GetNCells();
 346   Int_t   list[nc];
 347   Float_t elist[nc];
 348   Int_t nMax = GetCaloUtils()->GetNumberOfLocalMaxima(cluster, GetEMCALCells(),list, elist);
 349
 350 //  printf("AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin() - Cluster E %2.2f; NLM = %d, cluster MC labels:\n",cluster->E(),nMax);
 351 //
 352 //  for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
 353 //  {
 354 //    Bool_t ok  =kFALSE,gok = kFALSE;
 355 //    Int_t pdg    = -22222, status   = -1;
 356 //    Int_t gpdg   = -22222, gstatus  = -1;
 357 //    Int_t ggpdg  = -22222, ggstatus = -1;
 358 //    Int_t gLabel = -1, ggLabel = -1;
 359 //
 360 //    Int_t label = cluster->GetLabels()[ilab];
 361 //    TLorentzVector primary   =GetMCAnalysisUtils()->GetMother     (label,GetReader(),  pdg,  status, ok);
 362 //    TLorentzVector gprimary  =GetMCAnalysisUtils()->GetGrandMother(label,GetReader(), gpdg, gstatus,gok, gLabel,ggLabel);
 363 //    TLorentzVector ggprimary =GetMCAnalysisUtils()->GetMother(ggLabel  ,GetReader(),ggpdg,ggstatus,gok);
 364 //    printf("\t %d; mother: Label %d; PDG %d; E %2.2f - grand mother label %d; PDG %d; E %2.2f- great grand mother label %d; PDG %d; E %2.2f\n",
 365 //           ilab,label,pdg,primary.E(), gLabel,gpdg,gprimary.E(), ggLabel,ggpdg,ggprimary.E());
 366 //
 367 //  }
 368 //
 369 //  printf("AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin() - Cluster Cells MC labels:\n");
 370 //
 371 //  for (UInt_t icell = 0; icell < nc; icell++ )
 372 //  {
 373 //    Bool_t ok  =kFALSE,gok = kFALSE;
 374 //    Int_t pdg    = -22222, status   = -1;
 375 //    Int_t gpdg   = -22222, gstatus  = -1;
 376 //    Int_t ggpdg  = -22222, ggstatus = -1;
 377 //    Int_t gLabel = -1, ggLabel = -1;
 378 //
 379 //    Int_t absId = cluster->GetCellAbsId(icell);
 380 //
 381 //    printf("cell abs Id %d, amplitude %f\n",absId,GetEMCALCells()->GetCellAmplitude(absId));
 382 //
 383 //    Int_t label = GetEMCALCells()->GetCellMCLabel(absId);
 384 //    TLorentzVector primary   =GetMCAnalysisUtils()->GetMother     (label,GetReader(),  pdg,  status, ok);
 385 //    TLorentzVector gprimary  =GetMCAnalysisUtils()->GetGrandMother(label,GetReader(), gpdg, gstatus,gok, gLabel,ggLabel);
 386 //    TLorentzVector ggprimary =GetMCAnalysisUtils()->GetMother(ggLabel  ,GetReader(),ggpdg,ggstatus,gok);
 387 //    printf(" %d; mother: Label %d; PDG %d; E %2.2f - grand mother label %d; PDG %d; E %2.2f- great grand mother label %d; PDG %d; E %2.2f\n",
 388 //           icell,label,pdg,primary.E(), gLabel,gpdg,gprimary.E(), ggLabel,ggpdg,ggprimary.E());
 389
 390 //  }
 391
 392   //If only one maxima, consider all the towers in the cluster
 393   if(nMax==1)
 394   {
 395       for (UInt_t icell = 0; icell < nc; icell++ )
 396       {
 397         list [icell] = cluster->GetCellAbsId(icell);
 398         elist[icell] = GetEMCALCells()->GetCellAmplitude(list[icell]);
 399       }
 400   }
 401
 402   //Find highest energy Local Maxima Towers
 403   Int_t   imax  = -1;
 404   Int_t   imax2 = -1;
 405   Float_t emax  = -1;
 406   Float_t emax2 = -1;
 407   for(Int_t i = 0; i < nMax; i++)
 408   {
 409     //printf("i %d: AbsId %d; E %2.3f\n",i,list[i],elist[i]);
 410     if(elist[i] > emax)
 411     {
 412       imax = i;
 413       emax = elist[i];
 414     }
 415   }
 416   //Find second highest
 417   for(Int_t i = 0; i < nMax; i++)
 418   {
 419     if(i==imax) continue;
 420
 421     if(elist[i] > emax2)
 422     {
 423       imax2 = i;
 424       emax2 = elist[i];
 425     }
 426   }
 427
 428   //printf("Highest : %d and %d\n",imax,imax2);
 429
 430   // Check that the highest mc label and the max cluster label are the same
 431   Int_t mcLabelMax = GetEMCALCells()->GetCellMCLabel(list[imax]);
 432   GetReader()->RemapMCLabelForAODs(mcLabelMax);
 433   Int_t mcLabelMax2 = GetEMCALCells()->GetCellMCLabel(list[imax2]);
 434   GetReader()->RemapMCLabelForAODs(mcLabelMax2);
 435
 436   Int_t mcLabelclusterMax = cluster->GetLabels()[0];
 437   Bool_t matchHighLMAndHighMC = kFALSE;
 438
 439   if(mcLabelclusterMax==mcLabelMax)
 440   {
 441     matchHighLMAndHighMC = kTRUE;
 442     //printf("*** MATCH cluster and LM maximum ***\n");
 443   }
 444   else
 445   {
 446      //printf("*** NO MATCH cluster and LM maximum, check second ***\n");
 447     if(mcLabelclusterMax==mcLabelMax2)
 448     {
 449       //printf("\t *** MATCH cluster and 2nd LM maximum ***\n");
 450       matchHighLMAndHighMC = kTRUE;
 451     }
 452     else
 453     {
 454       //printf("\t *** NO MATCH***\n");
 455       matchHighLMAndHighMC = kFALSE;
 456     }
 457   }
 458
 459   // Compare the common ancestors of the 2 highest energy local maxima
 460   Int_t ancPDG = 0, ancStatus = -1;
 461   TLorentzVector momentum; TVector3 prodVertex;
 462   Int_t ancLabel = 0;
 463   Bool_t high = kFALSE;
 464   Bool_t low  = kFALSE;
 465
 466 //  // print maxima origin
 467 //  for(Int_t i = 0; i < nMax; i++)
 468 //  {
 469 //    Int_t mcLabel1 = GetEMCALCells()->GetCellMCLabel(list[i]);
 470 //    GetReader()->RemapMCLabelForAODs(mcLabel1);
 471 //
 472 //    Bool_t ok  =kFALSE,gok = kFALSE;
 473 //    Int_t pdg    = -22222, status   = -1;
 474 //    Int_t gpdg   = -22222, gstatus  = -1;
 475 //    Int_t ggpdg  = -22222, ggstatus = -1;
 476 //    Int_t gLabel = -1, ggLabel = -1;
 477 //    TLorentzVector primary   =GetMCAnalysisUtils()->GetMother     (mcLabel1,GetReader(),  pdg,  status, ok);
 478 //    TLorentzVector gprimary  =GetMCAnalysisUtils()->GetGrandMother(mcLabel1,GetReader(), gpdg, gstatus,gok, gLabel,ggLabel);
 479 //    TLorentzVector ggprimary =GetMCAnalysisUtils()->GetMother(ggLabel  ,GetReader(),ggpdg,ggstatus,gok);
 480 //    printf("Max index %d; mother: Label %d; PDG %d; E %2.2f - grand mother label %d; PDG %d; E %2.2f- great grand mother label %d; PDG %d; E %2.2f\n",
 481 //           i,mcLabel1,pdg,primary.E(), gLabel,gpdg,gprimary.E(), ggLabel,ggpdg,ggprimary.E());
 482 //  }
 483
 484   // Compare ancestors of all local maxima
 485   Int_t nmaxima = nMax;
 486   if(nMax==1) nmaxima = nc ;
 487   for(Int_t i = 0; i < nmaxima-1; i++)
 488   {
 489     Int_t mcLabel1 = GetEMCALCells()->GetCellMCLabel(list[i]);
 490     GetReader()->RemapMCLabelForAODs(mcLabel1);
 491
 492     for(Int_t j = i+1; j < nmaxima; j++)
 493     {
 494       Int_t mcLabel2 = GetEMCALCells()->GetCellMCLabel(list[j]);
 495       GetReader()->RemapMCLabelForAODs(mcLabel2);
 496
 497       if(mcLabel1 < 0 || mcLabel2 < 0 )
 498       {
 499         //printf("\t i %d label %d - j %d label %d; skip!\n",i,mcLabel1,j,mcLabel2);
 500         continue;
 501       }
 502       ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(mcLabel1,mcLabel2,
 503                                                            GetReader(),ancPDG,ancStatus,momentum,prodVertex);
 504       if(ancPDG==111)
 505       {
 506         if((i==imax && j==imax2) ||  (j==imax && i==imax2))
 507           high = kTRUE;
 508         else
 509           low = kTRUE;
 510       }
 511       else if(ancPDG==22 || TMath::Abs(ancPDG)==11)
 512       {
 513         // If both bits are set, it could be that one of the maxima had a conversion
 514         // reset the bit in this case
 515         if(high && low)
 516         {
 517           //printf("\t Reset low bit\n");
 518           low = kFALSE;
 519         }
 520       }
 521
 522       Bool_t ok  =kFALSE;
 523       Int_t pdg = -22222, status = -1;
 524       TLorentzVector primary  =GetMCAnalysisUtils()->GetMother(ancLabel,GetReader(), pdg, status, ok);
 525
 526       //printf("\t i %d label %d - j %d label %d; ancestor label %d, PDG %d-%d; E %2.2f; high %d, any %d \n",i,mcLabel1,j,mcLabel2, ancLabel, ancPDG,pdg, primary.E(), high, low);
 527
 528     }
 529   }
 530
 531   Float_t en = cluster->E();
 532
 533   //printf("nMax %d; Match MC? %d; high %d; low %d\n",nMax,matchHighLMAndHighMC,high,low);
 534
 535   if(matchHighLMAndHighMC)
 536   {
 537     if     (high && !low)  fhMCPi0HighNLMPair->Fill(en,nMax);
 538     else if(low  && !high) fhMCPi0LowNLMPair ->Fill(en,nMax);
 539     else if(low  &&  high) fhMCPi0AnyNLMPair ->Fill(en,nMax);
 540     else                   fhMCPi0NoneNLMPair->Fill(en,nMax);
 541   }
 542   else
 543   {
 544     if     (high && !low)  fhMCPi0HighNLMPairNoMCMatch->Fill(en,nMax);
 545     else if(low  && !high) fhMCPi0LowNLMPairNoMCMatch ->Fill(en,nMax);
 546     else if(low  &&  high) fhMCPi0AnyNLMPairNoMCMatch ->Fill(en,nMax);
 547     else                   fhMCPi0NoneNLMPairNoMCMatch->Fill(en,nMax);
 548   }
 549 }
 550
 551 //___________________________________________________________________________________________________________________
 552 void AliAnaInsideClusterInvariantMass::FillAngleHistograms(const Int_t nMax, const Bool_t matched,
 553                                                            const Float_t en, const Float_t angle, const Float_t mass)
 554 {
 555   // Fill histograms related to opening angle
 556
 557   if(!fFillAngleHisto) return;
 558
 559   if     (nMax==1)
 560   {
 561     fhAnglePairNLocMax1[matched]->Fill(en,angle);
 562     if( en > fHistoECut ) fhAnglePairMassNLocMax1[matched]->Fill(mass,angle);
 563   }
 564   else if(nMax==2)
 565   {
 566     fhAnglePairNLocMax2[matched]->Fill(en,angle);
 567     if( en > fHistoECut ) fhAnglePairMassNLocMax2[matched]->Fill(mass,angle);
 568   }
 569   else if(nMax >2)
 570   {
 571     fhAnglePairNLocMaxN[matched]->Fill(en,angle);
 572     if( en > fHistoECut ) fhAnglePairMassNLocMaxN[matched]->Fill(mass,angle);
 573   }
 574
 575 }
 576
 577 //__________________________________________________________________________________________________________________________________________
 578 void AliAnaInsideClusterInvariantMass::FillEBinHistograms(const Int_t   ebin     , const Int_t   nMax, const Int_t mcindex,
 579                                                           const Float_t splitFrac, const Float_t mass, const Float_t asym, const Float_t l0)
 580 {
 581   // Fill some histograms integrating in few energy bins
 582
 583   if(ebin < 0 || !fFillEbinHisto) return ;
 584
 585   if     (nMax==1)
 586   {
 587     fhMassSplitEFractionNLocMax1Ebin[0][ebin]->Fill(splitFrac,  mass);
 588     if(IsDataMC())fhMassSplitEFractionNLocMax1Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
 589
 590     fhMassM02NLocMax1Ebin    [ebin]->Fill(l0  ,  mass );
 591     fhMassAsyNLocMax1Ebin    [ebin]->Fill(asym,  mass );
 592   }
 593   else if(nMax==2)
 594   {
 595     fhMassSplitEFractionNLocMax2Ebin[0][ebin]->Fill(splitFrac,  mass);
 596     if(IsDataMC())fhMassSplitEFractionNLocMax2Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
 597
 598     fhMassM02NLocMax2Ebin    [ebin]->Fill(l0  ,  mass );
 599     fhMassAsyNLocMax2Ebin    [ebin]->Fill(asym,  mass );
 600   }
 601   else if(nMax > 2 )
 602   {
 603     fhMassSplitEFractionNLocMaxNEbin[0][ebin]->Fill(splitFrac,  mass);
 604     if(IsDataMC())fhMassSplitEFractionNLocMaxNEbin[mcindex][ebin]->Fill(splitFrac,  mass);
 605
 606     fhMassM02NLocMaxNEbin    [ebin]->Fill(l0  ,  mass );
 607     fhMassAsyNLocMaxNEbin    [ebin]->Fill(asym,  mass );
 608   }
 609
 610 }
 611
 612 //_____________________________________________________________________________________________________________________
 613 void AliAnaInsideClusterInvariantMass::FillMCHistograms(const Float_t en,        const Float_t e1  , const Float_t e2,
 614                                                         const Int_t ebin,        const Int_t mcindex,
 615                                                         const Float_t l0,        const Float_t mass,
 616                                                         const Int_t nMax,        const Bool_t  matched,
 617                                                         const Float_t splitFrac, const Float_t asym,
 618                                                         const Float_t eprim,     const Float_t asymGen)
 619 {
 620   // Fill histograms needing some MC input
 621
 622   if(!IsDataMC() || !fFillMCHisto) return;
 623
 624   Float_t efrac      = eprim/en;
 625   Float_t efracSplit = 0;
 626   if(e1+e2 > 0) efracSplit = eprim/(e1+e2);
 627
 628   //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",
 629   //       e1,e2,eprim,en,splitFrac,efrac,efracSplit);
 630
 631   if(ebin >= 0 && fFillEbinHisto)
 632   {
 633     if( !matched ) fhMCGenFracNLocMaxEbin       [mcindex][ebin]->Fill(efrac,nMax);
 634     else           fhMCGenFracNLocMaxEbinMatched[mcindex][ebin]->Fill(efrac,nMax);
 635   }
 636
 637   if     (nMax==1)
 638   {
 639     fhMCGenFracNLocMax1      [mcindex][matched]->Fill(en     ,  efrac );
 640     fhMCGenSplitEFracNLocMax1[mcindex][matched]->Fill(en     ,  efracSplit );
 641     fhMCGenEvsSplitENLocMax1 [mcindex][matched]->Fill(eprim  ,  e1+e2);
 642
 643     if( en > fHistoECut )
 644     {
 645       fhMCGenEFracvsSplitEFracNLocMax1[mcindex][matched]->Fill(efrac,splitFrac );
 646
 647       if(!matched && ebin >= 0 && fFillEbinHisto)
 648       {
 649         fhM02MCGenFracNLocMax1Ebin [mcindex][ebin]->Fill(efrac  ,  l0    );
 650         fhMassMCGenFracNLocMax1Ebin[mcindex][ebin]->Fill(efrac  ,  mass  );
 651
 652         fhMCAsymM02NLocMax1MCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
 653         fhAsyMCGenRecoNLocMax1EbinPi0[ebin]->Fill(asym,  asymGen );
 654       }
 655     }
 656   }
 657   else if(nMax==2)
 658   {
 659     fhMCGenFracNLocMax2      [mcindex][matched]->Fill(en     ,  efrac );
 660     fhMCGenSplitEFracNLocMax2[mcindex][matched]->Fill(en     ,  efracSplit );
 661     fhMCGenEvsSplitENLocMax2 [mcindex][matched]->Fill(eprim  ,  e1+e2);
 662
 663     if( en > fHistoECut )
 664     {
 665       fhMCGenEFracvsSplitEFracNLocMax2[mcindex][matched]->Fill(efrac,splitFrac );
 666
 667       if(!matched && ebin >= 0 && fFillEbinHisto)
 668       {
 669         fhM02MCGenFracNLocMax2Ebin [mcindex][ebin]->Fill(efrac  ,  l0    );
 670         fhMassMCGenFracNLocMax2Ebin[mcindex][ebin]->Fill(efrac  ,  mass  );
 671
 672         fhMCAsymM02NLocMax2MCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
 673         fhAsyMCGenRecoNLocMax2EbinPi0[ebin]->Fill(asym,  asymGen );
 674       }
 675     }
 676
 677   }
 678   else if(nMax > 2 )
 679   {
 680     fhMCGenFracNLocMaxN      [mcindex][matched]->Fill(en     ,  efrac );
 681     fhMCGenSplitEFracNLocMaxN[mcindex][matched]->Fill(en     ,  efracSplit );
 682     fhMCGenEvsSplitENLocMaxN [mcindex][matched]->Fill(eprim  ,  e1+e2);
 683
 684     if( en > fHistoECut )
 685     {
 686       fhMCGenEFracvsSplitEFracNLocMaxN[mcindex][matched]->Fill(efrac,splitFrac );
 687
 688       if(!matched && ebin >= 0 && fFillEbinHisto)
 689       {
 690         fhM02MCGenFracNLocMaxNEbin [mcindex][ebin]->Fill(efrac  ,  l0    );
 691         fhMassMCGenFracNLocMaxNEbin[mcindex][ebin]->Fill(efrac  ,  mass  );
 692
 693         fhMCAsymM02NLocMaxNMCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
 694         fhAsyMCGenRecoNLocMaxNEbinPi0[ebin]->Fill(asym,  asymGen );
 695       }
 696     }
 697   }
 698 }
 699
 700 //_________________________________________________________________________________________________________________________
 701 void AliAnaInsideClusterInvariantMass::FillMCOverlapHistograms(const Float_t en,      const Float_t mass, const Float_t l0,
 702                                                                const Int_t   inlm,    const Int_t ebin, const Bool_t matched,
 703                                                                const Int_t   mcindex, const Int_t noverlaps)
 704 {
 705
 706   // Fill histograms for MC Overlaps
 707
 708   //printf("en %f,mass %f,l0 %f,inlm %d,ebin %d,matched %d,mcindex %d,noverlaps %d \n",en,mass,l0,inlm,ebin,matched,mcindex,noverlaps);
 709
 710   if(!fFillMCOverlapHisto || !IsDataMC()) return;
 711
 712   //printf("AliAnaInsideClusterInvariantMass::FillMCOverlapHistograms - NLM bin=%d, mcIndex %d, n Overlaps %d\n",inlm,mcindex,noverlaps);
 713
 714   if(!matched)
 715   {
 716     fhMCENOverlaps[inlm][mcindex]->Fill(en,noverlaps);
 717
 718     if     (noverlaps == 0)
 719     {
 720       fhMCEM02Overlap0[inlm][mcindex]->Fill(en, l0);
 721       fhMCEMassOverlap0[inlm][mcindex]->Fill(en, mass);
 722       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap0[inlm][ebin]->Fill(l0,mass);
 723     }
 724     else if(noverlaps == 1)
 725     {
 726       fhMCEM02Overlap1[inlm][mcindex]->Fill(en, l0);
 727       fhMCEMassOverlap1[inlm][mcindex]->Fill(en, mass);
 728       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap1[inlm][ebin]->Fill(l0,mass);
 729     }
 730     else if(noverlaps  > 1)
 731     {
 732       fhMCEM02OverlapN[inlm][mcindex]->Fill(en, l0);
 733       fhMCEMassOverlapN[inlm][mcindex]->Fill(en, mass);
 734       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02OverlapN[inlm][ebin]->Fill(l0,mass);
 735     }
 736     else
 737       printf("AliAnaInsideClusterInvariantMass::FillMCOverlapHistograms() - n overlaps = %d!!", noverlaps);
 738   }
 739   else if(fFillTMHisto)
 740   {
 741     fhMCENOverlapsMatch[inlm][mcindex]->Fill(en,noverlaps);
 742
 743     if     (noverlaps == 0)
 744     {
 745       fhMCEM02Overlap0Match[inlm][mcindex]->Fill(en, l0);
 746       fhMCEMassOverlap0Match[inlm][mcindex]->Fill(en, mass);
 747       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap0Match[inlm][ebin]->Fill(l0,mass);
 748     }
 749     else if(noverlaps == 1)
 750     {
 751       fhMCEM02Overlap1Match[inlm][mcindex]->Fill(en, l0);
 752       fhMCEMassOverlap1Match[inlm][mcindex]->Fill(en, mass);
 753       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap1Match[inlm][ebin]->Fill(l0,mass);
 754     }
 755     else if(noverlaps  > 1)
 756     {
 757       fhMCEM02OverlapNMatch[inlm][mcindex]->Fill(en, l0);
 758       fhMCEMassOverlapNMatch[inlm][mcindex]->Fill(en, mass);
 759       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02OverlapNMatch[inlm][ebin]->Fill(l0,mass);
 760     }
 761     else
 762         printf("AliAnaInsideClusterInvariantMass::FillMCOverlapHistograms() - n overlaps in matched = %d!!", noverlaps);
 763   }
 764 }
 765
 766 //______________________________________________________________________________________________________
 767 void AliAnaInsideClusterInvariantMass::FillSSExtraHistograms(AliVCluster  *cluster, const Int_t nMax,
 768                                                              const Bool_t  matched, const Int_t mcindex,
 769                                                              const Float_t mass   , const Int_t ebin)
 770 {
 771   // Fill optional histograms with more SS parameters
 772
 773   if(!fFillSSExtraHisto) return ;
 774
 775   Float_t en = cluster->E();
 776   Float_t nc = cluster->GetNCells();
 777
 778   // Get more Shower Shape parameters
 779   Float_t ll0  = 0., ll1  = 0.;
 780   Float_t disp= 0., dispEta = 0., dispPhi    = 0.;
 781   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
 782
 783   GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
 784                                                                                ll0, ll1, disp, dispEta, dispPhi, sEta, sPhi, sEtaPhi);
 785
 786   Float_t dispAsy = -1;
 787   if(dispEta+dispPhi >0 ) dispAsy = (dispPhi-dispEta) / (dispPhi+dispEta);
 788
 789
 790   if     (nMax==1)
 791   {
 792     fhNCellNLocMax1[0][matched]->Fill(en,nc) ;
 793     if(mcindex > 0 )  fhNCellNLocMax1[mcindex][matched]->Fill(en,nc) ;
 794
 795     if( en > fHistoECut )
 796     {
 797       fhMassDispEtaNLocMax1[0][matched]->Fill(dispEta,  mass );
 798       fhMassDispPhiNLocMax1[0][matched]->Fill(dispPhi,  mass );
 799       fhMassDispAsyNLocMax1[0][matched]->Fill(dispAsy,  mass );
 800
 801       if(IsDataMC())
 802       {
 803         fhMassDispEtaNLocMax1[mcindex][matched]->Fill(dispEta,  mass );
 804         fhMassDispPhiNLocMax1[mcindex][matched]->Fill(dispPhi,  mass );
 805         fhMassDispAsyNLocMax1[mcindex][matched]->Fill(dispAsy,  mass );
 806       }
 807     }
 808
 809     if(!matched && ebin >= 0 && fFillEbinHisto)
 810     {
 811       fhMassDispEtaNLocMax1Ebin[ebin]->Fill(dispEta,  mass );
 812       fhMassDispPhiNLocMax1Ebin[ebin]->Fill(dispPhi,  mass );
 813       fhMassDispAsyNLocMax1Ebin[ebin]->Fill(dispAsy,  mass );
 814     }
 815   }
 816   else if( nMax == 2  )
 817   {
 818     fhNCellNLocMax2[0][matched]->Fill(en,nc) ;
 819     if(mcindex > 0 )  fhNCellNLocMax2[mcindex][matched]->Fill(en,nc) ;
 820
 821     if( en > fHistoECut )
 822     {
 823       fhMassDispEtaNLocMax2[0][matched]->Fill(dispEta,  mass );
 824       fhMassDispPhiNLocMax2[0][matched]->Fill(dispPhi,  mass );
 825       fhMassDispAsyNLocMax2[0][matched]->Fill(dispAsy,  mass );
 826
 827       if(IsDataMC())
 828       {
 829         fhMassDispEtaNLocMax2[mcindex][matched]->Fill(dispEta,  mass );
 830         fhMassDispPhiNLocMax2[mcindex][matched]->Fill(dispPhi,  mass );
 831         fhMassDispAsyNLocMax2[mcindex][matched]->Fill(dispAsy,  mass );
 832       }
 833     }
 834
 835     if(!matched && ebin >= 0 && fFillEbinHisto)
 836     {
 837       fhMassDispEtaNLocMax2Ebin[ebin]->Fill(dispEta,  mass );
 838       fhMassDispPhiNLocMax2Ebin[ebin]->Fill(dispPhi,  mass );
 839       fhMassDispAsyNLocMax2Ebin[ebin]->Fill(dispAsy,  mass );
 840     }
 841
 842   }
 843   else if( nMax >= 3  )
 844   {
 845     fhNCellNLocMaxN[0][matched]->Fill(en,nc) ;
 846     if(mcindex > 0 )  fhNCellNLocMaxN[mcindex][matched]->Fill(en,nc) ;
 847
 848     if( en > fHistoECut )
 849     {
 850       fhMassDispEtaNLocMaxN[0][matched]->Fill(dispEta,  mass );
 851       fhMassDispPhiNLocMaxN[0][matched]->Fill(dispPhi,  mass );
 852       fhMassDispAsyNLocMaxN[0][matched]->Fill(dispAsy,  mass );
 853
 854       if(IsDataMC())
 855       {
 856         fhMassDispEtaNLocMaxN[mcindex][matched]->Fill(dispEta,  mass );
 857         fhMassDispPhiNLocMaxN[mcindex][matched]->Fill(dispPhi,  mass );
 858         fhMassDispAsyNLocMaxN[mcindex][matched]->Fill(dispAsy,  mass );
 859       }
 860     }
 861
 862     if(!matched && ebin >= 0 && fFillEbinHisto)
 863     {
 864       fhMassDispEtaNLocMaxNEbin[ebin]->Fill(dispEta,  mass );
 865       fhMassDispPhiNLocMaxNEbin[ebin]->Fill(dispPhi,  mass );
 866       fhMassDispAsyNLocMaxNEbin[ebin]->Fill(dispAsy,  mass );
 867     }
 868
 869   }
 870
 871 }
 872
 873 //__________________________________________________________________________________________________
 874 void AliAnaInsideClusterInvariantMass::FillSSWeightHistograms(AliVCluster *clus,  const Int_t nlm,
 875                                                               const Int_t absId1, const Int_t absId2)
 876 {
 877   // Calculate weights and fill histograms
 878
 879   if(!fFillSSWeightHisto) return;
 880
 881   AliVCaloCells* cells = 0;
 882   if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
 883   else                        cells = GetPHOSCells();
 884
 885   // First recalculate energy in case non linearity was applied
 886   Float_t  energy = 0;
 887   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 888   {
 889
 890     Int_t id       = clus->GetCellsAbsId()[ipos];
 891
 892     //Recalibrate cell energy if needed
 893     Float_t amp = cells->GetCellAmplitude(id);
 894     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
 895
 896     energy    += amp;
 897
 898   } // energy loop
 899
 900   if(energy <=0 )
 901   {
 902     printf("AliAnaInsideClusterInvatiantMass::WeightHistograms()- Wrong calculated energy %f\n",energy);
 903     return;
 904   }
 905
 906   //Get amplitude of  main local maxima, recalibrate if needed
 907   Float_t amp1 = cells->GetCellAmplitude(absId1);
 908   GetCaloUtils()->RecalibrateCellAmplitude(amp1,fCalorimeter, absId1);
 909   Float_t amp2 = cells->GetCellAmplitude(absId2);
 910   GetCaloUtils()->RecalibrateCellAmplitude(amp2,fCalorimeter, absId2);
 911
 912   if(amp1 < amp2)        printf("Bad local maxima E ordering : id1 E %f, id2 E %f\n ",amp1,amp2);
 913   if(amp1==0 || amp2==0) printf("Null E local maxima : id1 E %f, id2 E %f\n "        ,amp1,amp2);
 914
 915   if(amp1>0)fhPi0CellEMaxEMax2Frac   [nlm]->Fill(energy,amp2/amp1);
 916   fhPi0CellEMaxClusterFrac [nlm]->Fill(energy,amp1/energy);
 917   fhPi0CellEMax2ClusterFrac[nlm]->Fill(energy,amp2/energy);
 918
 919   //Get the ratio and log ratio to all cells in cluster
 920   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 921   {
 922     Int_t id       = clus->GetCellsAbsId()[ipos];
 923
 924     //Recalibrate cell energy if needed
 925     Float_t amp = cells->GetCellAmplitude(id);
 926     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
 927
 928     if(amp > 0)fhPi0CellE       [nlm]->Fill(energy,amp);
 929     fhPi0CellEFrac   [nlm]->Fill(energy,amp/energy);
 930     fhPi0CellLogEFrac[nlm]->Fill(energy,TMath::Log(amp/energy));
 931
 932     if     (id!=absId1 && id!=absId2)
 933     {
 934       if(amp1>0)fhPi0CellEMaxFrac [nlm]->Fill(energy,amp/amp1);
 935       if(amp2>0)fhPi0CellEMax2Frac[nlm]->Fill(energy,amp/amp2);
 936     }
 937
 938   }
 939
 940   //Recalculate shower shape for different W0
 941   if(fCalorimeter=="EMCAL")
 942   {
 943     Float_t l0org = clus->GetM02();
 944     Float_t l1org = clus->GetM20();
 945     Float_t dorg  = clus->GetDispersion();
 946     Float_t w0org =  GetCaloUtils()->GetEMCALRecoUtils()->GetW0();
 947
 948     for(Int_t iw = 0; iw < fSSWeightN; iw++)
 949     {
 950       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(fSSWeight[iw]);
 951       //GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
 952
 953       Float_t l0   = 0., l1   = 0.;
 954       Float_t disp = 0., dEta = 0., dPhi    = 0.;
 955       Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
 956
 957       RecalculateClusterShowerShapeParametersWithCellCut(GetEMCALGeometry(), cells, clus,l0,l1,disp,
 958                                                          dEta, dPhi, sEta, sPhi, sEtaPhi,0);
 959
 960
 961       fhM02WeightPi0[nlm][iw]->Fill(energy,clus->GetM02());
 962
 963     } // w0 loop
 964
 965     // Set the original values back
 966     clus->SetM02(l0org);
 967     clus->SetM20(l1org);
 968     clus->SetDispersion(dorg);
 969     GetCaloUtils()->GetEMCALRecoUtils()->SetW0(w0org);
 970
 971     for(Int_t iec = 0; iec < fSSECellCutN; iec++)
 972     {
 973       Float_t l0   = 0., l1   = 0.;
 974       Float_t disp = 0., dEta = 0., dPhi    = 0.;
 975       Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
 976
 977       RecalculateClusterShowerShapeParametersWithCellCut(GetEMCALGeometry(), cells, clus,l0,l1,disp,
 978                                                          dEta, dPhi, sEta, sPhi, sEtaPhi,fSSECellCut[iec]);
 979
 980       //printf("E %f, l0 org %f, l0 new %f, slope %f\n",clus->E(),l0org,l0,fSSECellCut[iec]);
 981       fhM02ECellCutPi0[nlm][iec]->Fill(energy,l0);
 982
 983     } // w0 loop
 984
 985   }// EMCAL
 986 }
 987
 988 //________________________________________________________________________________________
 989 void  AliAnaInsideClusterInvariantMass::FillTrackMatchingHistograms(AliVCluster * cluster,
 990                                                                     const Int_t nMax,
 991                                                                     const Int_t mcindex)
 992 {
 993   // Fill histograms related to track matching
 994
 995   if(!fFillTMResidualHisto || !fFillTMHisto) return;
 996
 997   Float_t dZ  = cluster->GetTrackDz();
 998   Float_t dR  = cluster->GetTrackDx();
 999   Float_t en  = cluster->E();
1000
1001   if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1002   {
1003     dR = 2000., dZ = 2000.;
1004     GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1005   }
1006
1007   //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
1008
1009   if(TMath::Abs(dR) < 999)
1010   {
1011     if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1[0]->Fill(en,dR); }
1012     else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2[0]->Fill(en,dR); }
1013     else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxN[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxN[0]->Fill(en,dR); }
1014
1015     if(IsDataMC())
1016     {
1017       if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1[mcindex]->Fill(en,dR); }
1018       else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2[mcindex]->Fill(en,dR); }
1019       else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxN[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxN[mcindex]->Fill(en,dR); }
1020     }
1021
1022     AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
1023
1024     Bool_t positive = kFALSE;
1025     if(track) positive = (track->Charge()>0);
1026
1027     if(track)
1028     {
1029       if(positive)
1030       {
1031         if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1Pos[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1Pos[0]->Fill(en,dR); }
1032         else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2Pos[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2Pos[0]->Fill(en,dR); }
1033         else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxNPos[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxNPos[0]->Fill(en,dR); }
1034
1035         if(IsDataMC())
1036         {
1037           if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1Pos[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1Pos[mcindex]->Fill(en,dR); }
1038           else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2Pos[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2Pos[mcindex]->Fill(en,dR); }
1039           else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxNPos[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxNPos[mcindex]->Fill(en,dR); }
1040         }
1041       }
1042       else
1043       {
1044         if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1Neg[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1Neg[0]->Fill(en,dR); }
1045         else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2Neg[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2Neg[0]->Fill(en,dR); }
1046         else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxNNeg[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxNNeg[0]->Fill(en,dR); }
1047
1048         if(IsDataMC())
1049         {
1050           if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1Neg[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1Neg[mcindex]->Fill(en,dR); }
1051           else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2Neg[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2Neg[mcindex]->Fill(en,dR); }
1052           else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxNNeg[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxNNeg[mcindex]->Fill(en,dR); }
1053         }
1054       }
1055
1056     }// track exists
1057
1058   }
1059 }
1060
1061 //_______________________________________________________________
1062 TObjString *  AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
1063 {
1064         //Save parameters used for analysis
1065   TString parList ; //this will be list of parameters used for this analysis.
1066   const Int_t buffersize = 255;
1067   char onePar[buffersize] ;
1068
1069   snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
1070   parList+=onePar ;
1071
1072   snprintf(onePar,buffersize,"Calorimeter: %s\n",        fCalorimeter.Data()) ;
1073   parList+=onePar ;
1074   snprintf(onePar,buffersize,"fNLocMaxCutE =%2.2f \n",    GetCaloUtils()->GetLocalMaximaCutE()) ;
1075   parList+=onePar ;
1076   snprintf(onePar,buffersize,"fNLocMaxCutEDiff =%2.2f \n",GetCaloUtils()->GetLocalMaximaCutEDiff()) ;
1077   parList+=onePar ;
1078   snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n",    fM02MinCut, fM02MaxCut) ;
1079   parList+=onePar ;
1080   snprintf(onePar,buffersize,"fMinNCells =%d \n",        fMinNCells) ;
1081   parList+=onePar ;
1082   snprintf(onePar,buffersize,"fMinBadDist =%1.1f \n",    fMinBadDist) ;
1083   parList+=onePar ;
1084   if(fFillSSWeightHisto)
1085   {
1086     snprintf(onePar,buffersize," N w %d - N e cut %d \n",fSSWeightN,fSSECellCutN);
1087     parList+=onePar ;
1088   }
1089
1090   return new TObjString(parList) ;
1091
1092 }
1093
1094 //________________________________________________________________
1095 TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
1096 {
1097   // Create histograms to be saved in output file and
1098   // store them in outputContainer
1099   TList * outputContainer = new TList() ;
1100   outputContainer->SetName("InsideClusterHistos") ;
1101
1102   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
1103   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
1104   Int_t mbins    = GetHistogramRanges()->GetHistoMassBins();         Float_t mmax   = GetHistogramRanges()->GetHistoMassMax();         Float_t mmin   = GetHistogramRanges()->GetHistoMassMin();
1105   Int_t ncbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   ncmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   ncmin  = GetHistogramRanges()->GetHistoNClusterCellMin();
1106   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();          Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
1107   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
1108
1109   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
1110   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
1111   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
1112   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
1113   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
1114   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
1115
1116   TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron","#pi^{0} (#gamma->e^{#pm})"};
1117   TString pname[] ={"","Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron","Pi0Conv"};
1118
1119   Int_t n = 1;
1120
1121   if(IsDataMC()) n = 8;
1122
1123   Int_t nMaxBins = 10;
1124
1125   TString sMatched[] = {"","Matched"};
1126
1127   Int_t nMatched = 2;
1128   if(!fFillTMHisto) nMatched = 1;
1129
1130   for(Int_t i = 0; i < n; i++)
1131   {
1132     for(Int_t j = 0; j < nMatched; j++)
1133     {
1134
1135       fhMassNLocMax1[i][j]  = new TH2F(Form("hMassNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1136                                        Form("Invariant mass of splitted cluster with NLM=1 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
1137                                        nptbins,ptmin,ptmax,mbins,mmin,mmax);
1138       fhMassNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
1139       fhMassNLocMax1[i][j]->SetXTitle("E (GeV)");
1140       outputContainer->Add(fhMassNLocMax1[i][j]) ;
1141
1142       fhMassNLocMax2[i][j]  = new TH2F(Form("hMassNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1143                                        Form("Invariant mass of splitted cluster with NLM=2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
1144                                        nptbins,ptmin,ptmax,mbins,mmin,mmax);
1145       fhMassNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
1146       fhMassNLocMax2[i][j]->SetXTitle("E (GeV)");
1147       outputContainer->Add(fhMassNLocMax2[i][j]) ;
1148
1149       fhMassNLocMaxN[i][j]  = new TH2F(Form("hMassNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1150                                        Form("Invariant mass of splitted cluster with NLM>2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
1151                                        nptbins,ptmin,ptmax,mbins,mmin,mmax);
1152       fhMassNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
1153       fhMassNLocMaxN[i][j]->SetXTitle("E (GeV)");
1154       outputContainer->Add(fhMassNLocMaxN[i][j]) ;
1155
1156       if(i==0 && j==0)
1157       {
1158         fhMassSplitECutNLocMax1  = new TH2F("hMassSplitECutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, (E1+E2)/E cut, no TM",
1159                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
1160         fhMassSplitECutNLocMax1->SetYTitle("M (GeV/c^{2})");
1161         fhMassSplitECutNLocMax1->SetXTitle("E (GeV)");
1162         outputContainer->Add(fhMassSplitECutNLocMax1) ;
1163
1164         fhMassSplitECutNLocMax2  = new TH2F("hMassSplitECutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, (E1+E2)/E cut, no TM",
1165                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
1166         fhMassSplitECutNLocMax2->SetYTitle("M (GeV/c^{2})");
1167         fhMassSplitECutNLocMax2->SetXTitle("E (GeV)");
1168         outputContainer->Add(fhMassSplitECutNLocMax2) ;
1169
1170         fhMassSplitECutNLocMaxN  = new TH2F("hMassSplitECutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, (E1+E2)/E cut, no TM",
1171                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
1172         fhMassSplitECutNLocMaxN->SetYTitle("M (GeV/c^{2})");
1173         fhMassSplitECutNLocMaxN->SetXTitle("E (GeV)");
1174         outputContainer->Add(fhMassSplitECutNLocMaxN) ;
1175
1176         fhMassM02CutNLocMax1  = new TH2F("hMassM02CutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, M02 cut, no TM",
1177                                          nptbins,ptmin,ptmax,mbins,mmin,mmax);
1178         fhMassM02CutNLocMax1->SetYTitle("M (GeV/c^{2})");
1179         fhMassM02CutNLocMax1->SetXTitle("E (GeV)");
1180         outputContainer->Add(fhMassM02CutNLocMax1) ;
1181
1182         fhMassM02CutNLocMax2  = new TH2F("hMassM02CutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, M02 cut, no TM",
1183                                          nptbins,ptmin,ptmax,mbins,mmin,mmax);
1184         fhMassM02CutNLocMax2->SetYTitle("M (GeV/c^{2})");
1185         fhMassM02CutNLocMax2->SetXTitle("E (GeV)");
1186         outputContainer->Add(fhMassM02CutNLocMax2) ;
1187
1188         fhMassM02CutNLocMaxN  = new TH2F("hMassM02CutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, M02 cut, no TM",
1189                                          nptbins,ptmin,ptmax,mbins,mmin,mmax);
1190         fhMassM02CutNLocMaxN->SetYTitle("M (GeV/c^{2})");
1191         fhMassM02CutNLocMaxN->SetXTitle("E (GeV)");
1192         outputContainer->Add(fhMassM02CutNLocMaxN) ;
1193
1194       }
1195
1196       fhMassAfterCutsNLocMax1[i][j]     = new TH2F(Form("hMassAfterCutsNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1197                                                    Form("Mass vs E, %s %s, for N Local max = 1, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
1198                                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
1199       fhMassAfterCutsNLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1200       fhMassAfterCutsNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1201       outputContainer->Add(fhMassAfterCutsNLocMax1[i][j]) ;
1202
1203       fhMassAfterCutsNLocMax2[i][j]     = new TH2F(Form("hMassAfterCutsNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1204                                                    Form("Mass vs E, %s %s, for N Local max = 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
1205                                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
1206       fhMassAfterCutsNLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1207       fhMassAfterCutsNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1208       outputContainer->Add(fhMassAfterCutsNLocMax2[i][j]) ;
1209
1210
1211       fhMassAfterCutsNLocMaxN[i][j]     = new TH2F(Form("hMassAfterCutsNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1212                                                    Form("Mass vs E, %s %s, for N Local max > 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
1213                                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
1214       fhMassAfterCutsNLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1215       fhMassAfterCutsNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1216       outputContainer->Add(fhMassAfterCutsNLocMaxN[i][j]) ;
1217
1218       fhSplitEFractionAfterCutsNLocMax1[i][j]     = new TH2F(Form("hSplitEFractionAfterCutsNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1219                                                              Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 1, M02 and Asy cut on, %s %s",ptype[i].Data(),sMatched[j].Data()),
1220                                                              nptbins,ptmin,ptmax,120,0,1.2);
1221       fhSplitEFractionAfterCutsNLocMax1[i][j]   ->SetXTitle("E_{cluster} (GeV)");
1222       fhSplitEFractionAfterCutsNLocMax1[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1223       outputContainer->Add(fhSplitEFractionAfterCutsNLocMax1[i][j]) ;
1224
1225       fhSplitEFractionAfterCutsNLocMax2[i][j]     = new TH2F(Form("hSplitEFractionAfterCutsNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1226                                                              Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 2, M02 and Asy cut on, %s %s",ptype[i].Data(),sMatched[j].Data()),
1227                                                              nptbins,ptmin,ptmax,120,0,1.2);
1228       fhSplitEFractionAfterCutsNLocMax2[i][j]   ->SetXTitle("E_{cluster} (GeV)");
1229       fhSplitEFractionAfterCutsNLocMax2[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1230       outputContainer->Add(fhSplitEFractionAfterCutsNLocMax2[i][j]) ;
1231
1232       fhSplitEFractionAfterCutsNLocMaxN[i][j]    = new TH2F(Form("hSplitEFractionAfterCutsNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1233                                                             Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  > 2, M02 and Asy cut on, %s %s",ptype[i].Data(),sMatched[j].Data()),
1234                                                             nptbins,ptmin,ptmax,120,0,1.2);
1235       fhSplitEFractionAfterCutsNLocMaxN[i][j]   ->SetXTitle("E_{cluster} (GeV)");
1236       fhSplitEFractionAfterCutsNLocMaxN[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1237       outputContainer->Add(fhSplitEFractionAfterCutsNLocMaxN[i][j]) ;
1238
1239
1240       fhMassM02NLocMax1[i][j]  = new TH2F(Form("hMassM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1241                                           Form("Invariant mass of splitted cluster with NLM=1, #lambda_{0}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
1242                                           ssbins,ssmin,ssmax,mbins,mmin,mmax);
1243       fhMassM02NLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
1244       fhMassM02NLocMax1[i][j]->SetXTitle("#lambda_{0}^{2}");
1245       outputContainer->Add(fhMassM02NLocMax1[i][j]) ;
1246
1247       fhMassM02NLocMax2[i][j]  = new TH2F(Form("hMassM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1248                                           Form("Invariant mass of splitted cluster with NLM=2, #lambda_{0}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
1249                                           ssbins,ssmin,ssmax,mbins,mmin,mmax);
1250       fhMassM02NLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
1251       fhMassM02NLocMax2[i][j]->SetXTitle("#lambda_{0}^{2}");
1252       outputContainer->Add(fhMassM02NLocMax2[i][j]) ;
1253
1254       fhMassM02NLocMaxN[i][j]  = new TH2F(Form("hMassM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1255                                           Form("Invariant mass of splitted cluster with NLM>2, vs #lambda_{0}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
1256                                           ssbins,ssmin,ssmax,mbins,mmin,mmax);
1257       fhMassM02NLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
1258       fhMassM02NLocMaxN[i][j]->SetXTitle("#lambda_{0}^{2}");
1259       outputContainer->Add(fhMassM02NLocMaxN[i][j]) ;
1260
1261
1262       fhAsymNLocMax1[i][j]  = new TH2F(Form("hAsymNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1263                                     Form("Asymmetry of NLM=1  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
1264                                     nptbins,ptmin,ptmax,200,-1,1);
1265       fhAsymNLocMax1[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
1266       fhAsymNLocMax1[i][j]->SetXTitle("E (GeV)");
1267       outputContainer->Add(fhAsymNLocMax1[i][j]) ;
1268
1269       fhAsymNLocMax2[i][j]  = new TH2F(Form("hAsymNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1270                                     Form("Asymmetry of NLM=2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
1271                                     nptbins,ptmin,ptmax,200,-1,1);
1272       fhAsymNLocMax2[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
1273       fhAsymNLocMax2[i][j]->SetXTitle("E (GeV)");
1274       outputContainer->Add(fhAsymNLocMax2[i][j]) ;
1275
1276       fhAsymNLocMaxN[i][j]  = new TH2F(Form("hAsymNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1277                                     Form("Asymmetry of NLM>2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
1278                                     nptbins,ptmin,ptmax,200,-1,1);
1279       fhAsymNLocMaxN[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
1280       fhAsymNLocMaxN[i][j]->SetXTitle("E (GeV)");
1281       outputContainer->Add(fhAsymNLocMaxN[i][j]) ;
1282
1283       if(i==0 && j==0)
1284       {
1285         fhAsymM02CutNLocMax1  = new TH2F("hAsymM02CutNLocMax1","Asymmetry of NLM=1  vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
1286         fhAsymM02CutNLocMax1->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
1287         fhAsymM02CutNLocMax1->SetXTitle("E (GeV)");
1288         outputContainer->Add(fhAsymM02CutNLocMax1) ;
1289
1290         fhAsymM02CutNLocMax2  = new TH2F("hAsymM02CutNLocMax2","Asymmetry of NLM=2  vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
1291         fhAsymM02CutNLocMax2->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
1292         fhAsymM02CutNLocMax2->SetXTitle("E (GeV)");
1293         outputContainer->Add(fhAsymM02CutNLocMax2) ;
1294
1295         fhAsymM02CutNLocMaxN  = new TH2F("hAsymM02CutNLocMaxN","Asymmetry of NLM>2  vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
1296         fhAsymM02CutNLocMaxN->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
1297         fhAsymM02CutNLocMaxN->SetXTitle("E (GeV)");
1298         outputContainer->Add(fhAsymM02CutNLocMaxN) ;
1299       }
1300
1301       if(fFillSSExtraHisto)
1302       {
1303         fhMassDispEtaNLocMax1[i][j]  = new TH2F(Form("hMassDispEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1304                                                 Form("Invariant mass of splitted cluster with NLM=1, #sigma_{#eta #eta}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
1305                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
1306         fhMassDispEtaNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
1307         fhMassDispEtaNLocMax1[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
1308         outputContainer->Add(fhMassDispEtaNLocMax1[i][j]) ;
1309
1310         fhMassDispEtaNLocMax2[i][j]  = new TH2F(Form("hMassDispEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1311                                                 Form("Invariant mass of splitted cluster with NLM=2 #sigma_{#eta #eta}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
1312                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
1313         fhMassDispEtaNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
1314         fhMassDispEtaNLocMax2[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
1315         outputContainer->Add(fhMassDispEtaNLocMax2[i][j]) ;
1316
1317         fhMassDispEtaNLocMaxN[i][j]  = new TH2F(Form("hMassDispEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1318                                                 Form("Invariant mass of splitted cluster with NLM>2, #sigma_{#eta #eta}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
1319                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
1320         fhMassDispEtaNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
1321         fhMassDispEtaNLocMaxN[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
1322         outputContainer->Add(fhMassDispEtaNLocMaxN[i][j]) ;
1323
1324         fhMassDispPhiNLocMax1[i][j]  = new TH2F(Form("hMassDispPhiNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1325                                                 Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
1326                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
1327         fhMassDispPhiNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
1328         fhMassDispPhiNLocMax1[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
1329         outputContainer->Add(fhMassDispPhiNLocMax1[i][j]) ;
1330
1331         fhMassDispPhiNLocMax2[i][j]  = new TH2F(Form("hMassDispPhiNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1332                                                 Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
1333                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
1334         fhMassDispPhiNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
1335         fhMassDispPhiNLocMax2[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
1336         outputContainer->Add(fhMassDispPhiNLocMax2[i][j]) ;
1337
1338         fhMassDispPhiNLocMaxN[i][j]  = new TH2F(Form("hMassDispPhiNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1339                                                 Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
1340                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax);
1341         fhMassDispPhiNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
1342         fhMassDispPhiNLocMaxN[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
1343         outputContainer->Add(fhMassDispPhiNLocMaxN[i][j]) ;
1344
1345         fhMassDispAsyNLocMax1[i][j]  = new TH2F(Form("hMassDispAsyNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1346                                                 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 > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
1347                                                 200,-1,1,mbins,mmin,mmax);
1348         fhMassDispAsyNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
1349         fhMassDispAsyNLocMax1[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1350         outputContainer->Add(fhMassDispAsyNLocMax1[i][j]) ;
1351
1352         fhMassDispAsyNLocMax2[i][j]  = new TH2F(Form("hMassDispAsyNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1353                                                 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 > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
1354                                                 200,-1,1,mbins,mmin,mmax);
1355         fhMassDispAsyNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
1356         fhMassDispAsyNLocMax2[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1357         outputContainer->Add(fhMassDispAsyNLocMax2[i][j]) ;
1358
1359         fhMassDispAsyNLocMaxN[i][j]  = new TH2F(Form("hMassDispAsyNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1360                                                 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()),
1361                                                 200,-1,1,mbins,mmin,mmax);
1362         fhMassDispAsyNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
1363         fhMassDispAsyNLocMaxN[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1364         outputContainer->Add(fhMassDispAsyNLocMaxN[i][j]) ;
1365       }
1366
1367       fhNLocMax[i][j]     = new TH2F(Form("hNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
1368                                      Form("Number of local maxima in cluster %s %s",ptype[i].Data(),sMatched[j].Data()),
1369                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
1370       fhNLocMax[i][j]   ->SetYTitle("N maxima");
1371       fhNLocMax[i][j]   ->SetXTitle("E (GeV)");
1372       outputContainer->Add(fhNLocMax[i][j]) ;
1373
1374       fhNLocMaxM02Cut[i][j] = new TH2F(Form("hNLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
1375                                        Form("Number of local maxima in cluster %s for %2.2f < M02 < %2.2f",ptype[i].Data(),fM02MinCut,fM02MaxCut),
1376                                        nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
1377       fhNLocMaxM02Cut[i][j]->SetYTitle("N maxima");
1378       fhNLocMaxM02Cut[i][j]->SetXTitle("E (GeV)");
1379       outputContainer->Add(fhNLocMaxM02Cut[i][j]) ;
1380
1381
1382       fhM02NLocMax1[i][j]     = new TH2F(Form("hM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1383                                          Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
1384                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1385       fhM02NLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1386       fhM02NLocMax1[i][j]   ->SetXTitle("E (GeV)");
1387       outputContainer->Add(fhM02NLocMax1[i][j]) ;
1388
1389       fhM02NLocMax2[i][j]     = new TH2F(Form("hM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1390                                          Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1391                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1392       fhM02NLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1393       fhM02NLocMax2[i][j]   ->SetXTitle("E (GeV)");
1394       outputContainer->Add(fhM02NLocMax2[i][j]) ;
1395
1396       fhM02NLocMaxN[i][j]    = new TH2F(Form("hM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1397                                         Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1398                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1399       fhM02NLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1400       fhM02NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1401       outputContainer->Add(fhM02NLocMaxN[i][j]) ;
1402
1403
1404       fhSplitEFractionNLocMax1[i][j]     = new TH2F(Form("hSplitEFractionNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1405                                          Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
1406                                          nptbins,ptmin,ptmax,120,0,1.2);
1407       fhSplitEFractionNLocMax1[i][j]   ->SetXTitle("E_{cluster} (GeV)");
1408       fhSplitEFractionNLocMax1[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1409       outputContainer->Add(fhSplitEFractionNLocMax1[i][j]) ;
1410
1411       fhSplitEFractionNLocMax2[i][j]     = new TH2F(Form("hSplitEFractionNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1412                                          Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1413                                          nptbins,ptmin,ptmax,120,0,1.2);
1414       fhSplitEFractionNLocMax2[i][j]   ->SetXTitle("E_{cluster} (GeV)");
1415       fhSplitEFractionNLocMax2[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1416       outputContainer->Add(fhSplitEFractionNLocMax2[i][j]) ;
1417
1418       fhSplitEFractionNLocMaxN[i][j]    = new TH2F(Form("hSplitEFractionNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1419                                         Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1420                                         nptbins,ptmin,ptmax,120,0,1.2);
1421       fhSplitEFractionNLocMaxN[i][j]   ->SetXTitle("E_{cluster} (GeV)");
1422       fhSplitEFractionNLocMaxN[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1423       outputContainer->Add(fhSplitEFractionNLocMaxN[i][j]) ;
1424
1425
1426       if(i > 0 && fFillMCHisto) // skip first entry in array, general case not filled
1427       {
1428         fhMCGenFracNLocMax1[i][j]     = new TH2F(Form("hMCGenFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1429                                                  Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
1430                                                  nptbins,ptmin,ptmax,200,0,2);
1431         fhMCGenFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
1432         fhMCGenFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1433         outputContainer->Add(fhMCGenFracNLocMax1[i][j]) ;
1434
1435         fhMCGenFracNLocMax2[i][j]     = new TH2F(Form("hMCGenFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1436                                                  Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1437                                                  nptbins,ptmin,ptmax,200,0,2);
1438         fhMCGenFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
1439         fhMCGenFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1440         outputContainer->Add(fhMCGenFracNLocMax2[i][j]) ;
1441
1442
1443         fhMCGenFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1444                                                 Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1445                                                 nptbins,ptmin,ptmax,200,0,2);
1446         fhMCGenFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
1447         fhMCGenFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1448         outputContainer->Add(fhMCGenFracNLocMaxN[i][j]) ;
1449
1450         fhMCGenSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1451                                                  Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
1452                                                  nptbins,ptmin,ptmax,200,0,2);
1453         fhMCGenSplitEFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
1454         fhMCGenSplitEFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1455         outputContainer->Add(fhMCGenSplitEFracNLocMax1[i][j]) ;
1456
1457         fhMCGenSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1458                                                  Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1459                                                  nptbins,ptmin,ptmax,200,0,2);
1460         fhMCGenSplitEFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
1461         fhMCGenSplitEFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1462         outputContainer->Add(fhMCGenSplitEFracNLocMax2[i][j]) ;
1463
1464
1465         fhMCGenSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1466                                                 Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1467                                                 nptbins,ptmin,ptmax,200,0,2);
1468         fhMCGenSplitEFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
1469         fhMCGenSplitEFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1470         outputContainer->Add(fhMCGenSplitEFracNLocMaxN[i][j]) ;
1471
1472         fhMCGenEFracvsSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1473                                                        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()),
1474                                                        200,0,2,200,0,2);
1475         fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
1476         fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1477         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax1[i][j]) ;
1478
1479         fhMCGenEFracvsSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1480                                                        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()),
1481                                                        200,0,2,200,0,2);
1482         fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
1483         fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1484         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax2[i][j]) ;
1485
1486
1487         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1488                                                       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()),
1489                                                       200,0,2,200,0,2);
1490         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
1491         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1492         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMaxN[i][j]) ;
1493
1494
1495         fhMCGenEvsSplitENLocMax1[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1496                                                              Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
1497                                                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1498         fhMCGenEvsSplitENLocMax1[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
1499         fhMCGenEvsSplitENLocMax1[i][j]   ->SetXTitle("E_{gen} (GeV)");
1500         outputContainer->Add(fhMCGenEvsSplitENLocMax1[i][j]) ;
1501
1502         fhMCGenEvsSplitENLocMax2[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1503                                                              Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1504                                                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1505         fhMCGenEvsSplitENLocMax2[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
1506         fhMCGenEvsSplitENLocMax2[i][j]   ->SetXTitle("E_{gen} (GeV)");
1507         outputContainer->Add(fhMCGenEvsSplitENLocMax2[i][j]) ;
1508
1509
1510         fhMCGenEvsSplitENLocMaxN[i][j]    = new TH2F(Form("hMCGenEvsSplitENLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1511                                                             Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1512                                                             nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1513         fhMCGenEvsSplitENLocMaxN[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
1514         fhMCGenEvsSplitENLocMaxN[i][j]   ->SetXTitle("E_{gen} (GeV)");
1515         outputContainer->Add(fhMCGenEvsSplitENLocMaxN[i][j]) ;
1516
1517       }
1518
1519       if(fFillSSExtraHisto)
1520       {
1521         fhNCellNLocMax1[i][j]  = new TH2F(Form("hNCellNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1522                                           Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
1523                                           nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
1524         fhNCellNLocMax1[i][j] ->SetYTitle("N cells");
1525         fhNCellNLocMax1[i][j] ->SetXTitle("E (GeV)");
1526         outputContainer->Add(fhNCellNLocMax1[i][j]) ;
1527
1528         fhNCellNLocMax2[i][j]     = new TH2F(Form("hNCellNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1529                                              Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1530                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
1531         fhNCellNLocMax2[i][j]   ->SetYTitle("N cells");
1532         fhNCellNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1533         outputContainer->Add(fhNCellNLocMax2[i][j]) ;
1534
1535
1536         fhNCellNLocMaxN[i][j]     = new TH2F(Form("hNCellNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1537                                              Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
1538                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
1539         fhNCellNLocMaxN[i][j]   ->SetYTitle("N cells");
1540         fhNCellNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1541         outputContainer->Add(fhNCellNLocMaxN[i][j]) ;
1542       }
1543
1544
1545       // E vs centrality
1546
1547       fhCentralityPi0NLocMax1[i][j]  = new TH2F(Form("hCentralityPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1548                                                 Form("E vs Centrality, selected pi0 cluster with NLM=1, %s",ptype[i].Data()),
1549                                                 nptbins,ptmin,ptmax,100,0,100);
1550       fhCentralityPi0NLocMax1[i][j]->SetYTitle("Centrality");
1551       fhCentralityPi0NLocMax1[i][j]->SetXTitle("E (GeV)");
1552       outputContainer->Add(fhCentralityPi0NLocMax1[i][j]) ;
1553
1554       fhCentralityPi0NLocMax2[i][j]  = new TH2F(Form("hCentralityPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1555                                                 Form("E vs Centrality, selected pi0 cluster with NLM=2, %s",ptype[i].Data()),
1556                                                 nptbins,ptmin,ptmax,100,0,100);
1557       fhCentralityPi0NLocMax2[i][j]->SetYTitle("Centrality");
1558       fhCentralityPi0NLocMax2[i][j]->SetXTitle("E (GeV)");
1559       outputContainer->Add(fhCentralityPi0NLocMax2[i][j]) ;
1560
1561       fhCentralityPi0NLocMaxN[i][j]  = new TH2F(Form("hCentralityPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1562                                                 Form("E vs Centrality, selected pi0 cluster with NLM>1, %s",ptype[i].Data()),
1563                                                 nptbins,ptmin,ptmax,100,0,100);
1564       fhCentralityPi0NLocMaxN[i][j]->SetYTitle("Centrality");
1565       fhCentralityPi0NLocMaxN[i][j]->SetXTitle("E (GeV)");
1566       outputContainer->Add(fhCentralityPi0NLocMaxN[i][j]) ;
1567
1568       fhCentralityEtaNLocMax1[i][j]  = new TH2F(Form("hCentralityEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1569                                                 Form("E vs Centrality, selected pi0 cluster with NLM=1, %s",ptype[i].Data()),
1570                                                 nptbins,ptmin,ptmax,100,0,100);
1571       fhCentralityEtaNLocMax1[i][j]->SetYTitle("Centrality");
1572       fhCentralityEtaNLocMax1[i][j]->SetXTitle("E (GeV)");
1573       outputContainer->Add(fhCentralityEtaNLocMax1[i][j]) ;
1574
1575       fhCentralityEtaNLocMax2[i][j]  = new TH2F(Form("hCentralityEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1576                                                 Form("E vs Centrality, selected pi0 cluster with NLM=2, %s",ptype[i].Data()),
1577                                                 nptbins,ptmin,ptmax,100,0,100);
1578       fhCentralityEtaNLocMax2[i][j]->SetYTitle("Centrality");
1579       fhCentralityEtaNLocMax2[i][j]->SetXTitle("E (GeV)");
1580       outputContainer->Add(fhCentralityEtaNLocMax2[i][j]) ;
1581
1582       fhCentralityEtaNLocMaxN[i][j]  = new TH2F(Form("hCentralityEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1583                                                 Form("E vs Centrality, selected pi0 cluster with NLM>1, %s",ptype[i].Data()),
1584                                                 nptbins,ptmin,ptmax,100,0,100);
1585       fhCentralityEtaNLocMaxN[i][j]->SetYTitle("Centrality");
1586       fhCentralityEtaNLocMaxN[i][j]->SetXTitle("E (GeV)");
1587       outputContainer->Add(fhCentralityEtaNLocMaxN[i][j]) ;
1588
1589
1590       fhM02Pi0NLocMax1[i][j]     = new TH2F(Form("hM02Pi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1591                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 1",
1592                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1593                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1594       fhM02Pi0NLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1595       fhM02Pi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
1596       outputContainer->Add(fhM02Pi0NLocMax1[i][j]) ;
1597
1598       fhM02EtaNLocMax1[i][j]     = new TH2F(Form("hM02EtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1599                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
1600                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1601                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1602       fhM02EtaNLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1603       fhM02EtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1604       outputContainer->Add(fhM02EtaNLocMax1[i][j]) ;
1605
1606       fhM02ConNLocMax1[i][j]    = new TH2F(Form("hM02ConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1607                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
1608                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1609                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1610       fhM02ConNLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1611       fhM02ConNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1612       outputContainer->Add(fhM02ConNLocMax1[i][j]) ;
1613
1614       fhM02Pi0NLocMax2[i][j]     = new TH2F(Form("hM02Pi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1615                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 2",
1616                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1617                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1618       fhM02Pi0NLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1619       fhM02Pi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
1620       outputContainer->Add(fhM02Pi0NLocMax2[i][j]) ;
1621
1622       fhM02EtaNLocMax2[i][j]     = new TH2F(Form("hM02EtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1623                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
1624                                                GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1625                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1626       fhM02EtaNLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1627       fhM02EtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1628       outputContainer->Add(fhM02EtaNLocMax2[i][j]) ;
1629
1630       fhM02ConNLocMax2[i][j]    = new TH2F(Form("hM02ConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1631                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
1632                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1633                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1634       fhM02ConNLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1635       fhM02ConNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1636       outputContainer->Add(fhM02ConNLocMax2[i][j]) ;
1637
1638       fhM02Pi0NLocMaxN[i][j]     = new TH2F(Form("hM02Pi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1639                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max > 2",
1640                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1641                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1642       fhM02Pi0NLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1643       fhM02Pi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1644       outputContainer->Add(fhM02Pi0NLocMaxN[i][j]) ;
1645
1646       fhM02EtaNLocMaxN[i][j]     = new TH2F(Form("hM02EtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1647                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max > 2",
1648                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1649                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1650       fhM02EtaNLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1651       fhM02EtaNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1652       outputContainer->Add(fhM02EtaNLocMaxN[i][j]) ;
1653
1654       fhM02ConNLocMaxN[i][j]    = new TH2F(Form("hM02ConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1655                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
1656                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1657                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1658       fhM02ConNLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1659       fhM02ConNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1660       outputContainer->Add(fhM02ConNLocMaxN[i][j]) ;
1661
1662
1663       fhMassPi0NLocMax1[i][j]     = new TH2F(Form("hMassPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1664                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 1",
1665                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1666                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
1667       fhMassPi0NLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1668       fhMassPi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
1669       outputContainer->Add(fhMassPi0NLocMax1[i][j]) ;
1670
1671
1672       fhMassEtaNLocMax1[i][j]     = new TH2F(Form("hMassEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1673                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
1674                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1675                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
1676       fhMassEtaNLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1677       fhMassEtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1678       outputContainer->Add(fhMassEtaNLocMax1[i][j]) ;
1679
1680       fhMassConNLocMax1[i][j]    = new TH2F(Form("hMassConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1681                                           Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
1682                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1683                                           nptbins,ptmin,ptmax,mbins,mmin,mmax);
1684       fhMassConNLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1685       fhMassConNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1686       outputContainer->Add(fhMassConNLocMax1[i][j]) ;
1687
1688       fhMassPi0NLocMax2[i][j]     = new TH2F(Form("hMassPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1689                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 2",
1690                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1691                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
1692       fhMassPi0NLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1693       fhMassPi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
1694       outputContainer->Add(fhMassPi0NLocMax2[i][j]) ;
1695
1696
1697       fhMassEtaNLocMax2[i][j]     = new TH2F(Form("hMassEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1698                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
1699                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1700                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
1701       fhMassEtaNLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1702       fhMassEtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1703       outputContainer->Add(fhMassEtaNLocMax2[i][j]) ;
1704
1705       fhMassConNLocMax2[i][j]    = new TH2F(Form("hMassConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1706                                           Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
1707                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1708                                           nptbins,ptmin,ptmax,mbins,mmin,mmax);
1709       fhMassConNLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1710       fhMassConNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1711       outputContainer->Add(fhMassConNLocMax2[i][j]) ;
1712
1713       fhMassPi0NLocMaxN[i][j]     = new TH2F(Form("hMassPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1714                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max > 2",
1715                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1716                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
1717       fhMassPi0NLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1718       fhMassPi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1719       outputContainer->Add(fhMassPi0NLocMaxN[i][j]) ;
1720
1721       fhMassEtaNLocMaxN[i][j]     = new TH2F(Form("hMassEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1722                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max > 2",
1723                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1724                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
1725       fhMassEtaNLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1726       fhMassEtaNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1727       outputContainer->Add(fhMassEtaNLocMaxN[i][j]) ;
1728
1729       fhMassConNLocMaxN[i][j]    = new TH2F(Form("hMassConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1730                                           Form("Mass vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
1731                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1732                                           nptbins,ptmin,ptmax,mbins,mmin,mmax);
1733       fhMassConNLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1734       fhMassConNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1735       outputContainer->Add(fhMassConNLocMaxN[i][j]) ;
1736
1737
1738       fhAsyPi0NLocMax1[i][j]     = new TH2F(Form("hAsyPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1739                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 1",
1740                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1741                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
1742       fhAsyPi0NLocMax1[i][j]   ->SetYTitle("Asymmetry");
1743       fhAsyPi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
1744       outputContainer->Add(fhAsyPi0NLocMax1[i][j]) ;
1745
1746       fhAsyEtaNLocMax1[i][j]     = new TH2F(Form("hAsyEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1747                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
1748                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1749                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
1750       fhAsyEtaNLocMax1[i][j]   ->SetYTitle("Asymmetry");
1751       fhAsyEtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1752       outputContainer->Add(fhAsyEtaNLocMax1[i][j]) ;
1753
1754       fhAsyConNLocMax1[i][j]    = new TH2F(Form("hAsyConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1755                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
1756                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1757                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
1758       fhAsyConNLocMax1[i][j]   ->SetYTitle("Asymmetry");
1759       fhAsyConNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1760       outputContainer->Add(fhAsyConNLocMax1[i][j]) ;
1761
1762       fhAsyPi0NLocMax2[i][j]     = new TH2F(Form("hAsyPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1763                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 2",
1764                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1765                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
1766       fhAsyPi0NLocMax2[i][j]   ->SetYTitle("Asymmetry");
1767       fhAsyPi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
1768       outputContainer->Add(fhAsyPi0NLocMax2[i][j]) ;
1769
1770       fhAsyEtaNLocMax2[i][j]     = new TH2F(Form("hAsyEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1771                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
1772                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1773                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
1774       fhAsyEtaNLocMax2[i][j]   ->SetYTitle("Asymmetry");
1775       fhAsyEtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1776       outputContainer->Add(fhAsyEtaNLocMax2[i][j]) ;
1777
1778       fhAsyConNLocMax2[i][j]    = new TH2F(Form("hAsyConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1779                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
1780                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1781                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
1782       fhAsyConNLocMax2[i][j]   ->SetYTitle("Asymmetry");
1783       fhAsyConNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1784       outputContainer->Add(fhAsyConNLocMax2[i][j]) ;
1785
1786       fhAsyPi0NLocMaxN[i][j]     = new TH2F(Form("hAsyPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1787                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max > 2",
1788                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1789                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
1790       fhAsyPi0NLocMaxN[i][j]   ->SetYTitle("Asymmetry");
1791       fhAsyPi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1792       outputContainer->Add(fhAsyPi0NLocMaxN[i][j]) ;
1793
1794       fhAsyEtaNLocMaxN[i][j]     = new TH2F(Form("hAsyEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1795                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max > 2",
1796                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1797                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
1798       fhAsyEtaNLocMaxN[i][j]   ->SetYTitle("Asymmetry");
1799       fhAsyEtaNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1800       outputContainer->Add(fhAsyEtaNLocMaxN[i][j]) ;
1801
1802       fhAsyConNLocMaxN[i][j]    = new TH2F(Form("hAsyConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1803                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
1804                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1805                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
1806       fhAsyConNLocMaxN[i][j]   ->SetYTitle("Asymmetry");
1807       fhAsyConNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1808       outputContainer->Add(fhAsyConNLocMaxN[i][j]) ;
1809
1810     } // matched, not matched
1811
1812     if(fFillEbinHisto)
1813     {
1814       for(Int_t j = 0; j < 4; j++)
1815       {
1816
1817         fhMassSplitEFractionNLocMax1Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax1%sEbin%d",pname[i].Data(),j),
1818                                                            Form("Invariant mass of 2 highest energy cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
1819                                                            120,0,1.2,mbins,mmin,mmax);
1820         fhMassSplitEFractionNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
1821         fhMassSplitEFractionNLocMax1Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
1822         outputContainer->Add(fhMassSplitEFractionNLocMax1Ebin[i][j]) ;
1823
1824         fhMassSplitEFractionNLocMax2Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax2%sEbin%d",pname[i].Data(),j),
1825                                                            Form("Invariant mass of 2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
1826                                                            120,0,1.2,mbins,mmin,mmax);
1827         fhMassSplitEFractionNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
1828         fhMassSplitEFractionNLocMax2Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
1829         outputContainer->Add(fhMassSplitEFractionNLocMax2Ebin[i][j]) ;
1830
1831         fhMassSplitEFractionNLocMaxNEbin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMaxN%sEbin%d",pname[i].Data(),j),
1832                                                            Form("Invariant mass of N>2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
1833                                                            120,0,1.2,mbins,mmin,mmax);
1834         fhMassSplitEFractionNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
1835         fhMassSplitEFractionNLocMaxNEbin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
1836         outputContainer->Add(fhMassSplitEFractionNLocMaxNEbin[i][j]) ;
1837
1838         if(i>0 && fFillMCHisto) // skip first entry in array, general case not filled
1839         {
1840           fhMCGenFracNLocMaxEbin[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%d",pname[i].Data(),j),
1841                                                    Form("NLM vs E, %s, E bin %d",ptype[i].Data(),j),
1842                                                    200,0,2,nMaxBins,0,nMaxBins);
1843           fhMCGenFracNLocMaxEbin[i][j]->SetYTitle("NLM");
1844           fhMCGenFracNLocMaxEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
1845           outputContainer->Add(fhMCGenFracNLocMaxEbin[i][j]) ;
1846
1847           fhMCGenFracNLocMaxEbinMatched[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%dMatched",pname[i].Data(),j),
1848                                                           Form("NLM vs E, %s, E bin %d, matched to a track",ptype[i].Data(),j),
1849                                                           200,0,2,nMaxBins,0,nMaxBins);
1850           fhMCGenFracNLocMaxEbinMatched[i][j]->SetYTitle("NLM");
1851           fhMCGenFracNLocMaxEbinMatched[i][j]->SetXTitle("E_{gen} / E_{reco}");
1852           outputContainer->Add(fhMCGenFracNLocMaxEbinMatched[i][j]) ;
1853
1854           fhMassMCGenFracNLocMax1Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
1855                                                         Form("Invariant mass of 2 highest energy cells vs E, %s, E bin %d",ptype[i].Data(),j),
1856                                                         200,0,2,mbins,mmin,mmax);
1857           fhMassMCGenFracNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
1858           fhMassMCGenFracNLocMax1Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
1859           outputContainer->Add(fhMassMCGenFracNLocMax1Ebin[i][j]) ;
1860
1861           fhMassMCGenFracNLocMax2Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
1862                                                         Form("Invariant mass of 2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
1863                                                         200,0,2,mbins,mmin,mmax);
1864           fhMassMCGenFracNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
1865           fhMassMCGenFracNLocMax2Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
1866           outputContainer->Add(fhMassMCGenFracNLocMax2Ebin[i][j]) ;
1867
1868           fhMassMCGenFracNLocMaxNEbin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
1869                                                         Form("Invariant mass of N>2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
1870                                                         200,0,2,mbins,mmin,mmax);
1871           fhMassMCGenFracNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
1872           fhMassMCGenFracNLocMaxNEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
1873           outputContainer->Add(fhMassMCGenFracNLocMaxNEbin[i][j]) ;
1874
1875           fhM02MCGenFracNLocMax1Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
1876                                                           Form("#lambda_{0}^{2} vs E for N max  = 1 %s, E bin %d",ptype[i].Data(), j),
1877                                                           200,0,2,ssbins,ssmin,ssmax);
1878           fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1879           fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1880           outputContainer->Add(fhM02MCGenFracNLocMax1Ebin[i][j]) ;
1881
1882           fhM02MCGenFracNLocMax2Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
1883                                                           Form("#lambda_{0}^{2} vs E for N max  = 2 %s, E bin %d",ptype[i].Data(),j),
1884                                                           200,0,2,ssbins,ssmin,ssmax);
1885           fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1886           fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1887           outputContainer->Add(fhM02MCGenFracNLocMax2Ebin[i][j]) ;
1888
1889           fhM02MCGenFracNLocMaxNEbin[i][j]    = new TH2F(Form("hM02MCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
1890                                                          Form("#lambda_{0}^{2} vs E for N max  > 2 %s, E bin %d",ptype[i].Data(),j),
1891                                                          200,0,2,ssbins,ssmin,ssmax);
1892           fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1893           fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1894           outputContainer->Add(fhM02MCGenFracNLocMaxNEbin[i][j]) ;
1895         }
1896       }
1897     }
1898   } // MC particle list
1899
1900   // E vs Event plane angle
1901
1902   fhEventPlanePi0NLocMax1  = new TH2F("hEventPlanePi0NLocMax1","E vs Event Plane Angle, selected pi0 cluster with NLM=1",
1903                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1904   fhEventPlanePi0NLocMax1->SetYTitle("Event Plane Angle (rad)");
1905   fhEventPlanePi0NLocMax1->SetXTitle("E (GeV)");
1906   outputContainer->Add(fhEventPlanePi0NLocMax1) ;
1907
1908   fhEventPlanePi0NLocMax2  = new TH2F("hEventPlanePi0NLocMax2","E vs Event Plane Angle, selected pi0 cluster with NLM=2",
1909                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1910   fhEventPlanePi0NLocMax2->SetYTitle("Event Plane Angle (rad)");
1911   fhEventPlanePi0NLocMax2->SetXTitle("E (GeV)");
1912   outputContainer->Add(fhEventPlanePi0NLocMax2) ;
1913
1914   fhEventPlanePi0NLocMaxN  = new TH2F("hEventPlanePi0NLocMaxN","E vs Event Plane Angle, selected pi0 cluster with NLM>1",
1915                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1916   fhEventPlanePi0NLocMaxN->SetYTitle("Event Plane Angle (rad)");
1917   fhEventPlanePi0NLocMaxN->SetXTitle("E (GeV)");
1918   outputContainer->Add(fhEventPlanePi0NLocMaxN) ;
1919
1920   fhEventPlaneEtaNLocMax1  = new TH2F("hEventPlaneEtaNLocMax1","E vs Event Plane Angle, selected pi0 cluster with NLM=1",
1921                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1922   fhEventPlaneEtaNLocMax1->SetYTitle("Event Plane Angle (rad)");
1923   fhEventPlaneEtaNLocMax1->SetXTitle("E (GeV)");
1924   outputContainer->Add(fhEventPlaneEtaNLocMax1) ;
1925
1926   fhEventPlaneEtaNLocMax2  = new TH2F("hEventPlaneEtaNLocMax2","E vs Event Plane Angle, selected pi0 cluster with NLM=2",
1927                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1928   fhEventPlaneEtaNLocMax2->SetYTitle("Event Plane Angle (rad)");
1929   fhEventPlaneEtaNLocMax2->SetXTitle("E (GeV)");
1930   outputContainer->Add(fhEventPlaneEtaNLocMax2) ;
1931
1932   fhEventPlaneEtaNLocMaxN  = new TH2F("hEventPlaneEtaNLocMaxN","E vs Event Plane Angle, selected pi0 cluster with NLM>1",
1933                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1934   fhEventPlaneEtaNLocMaxN->SetYTitle("Event Plane Angle (rad)");
1935   fhEventPlaneEtaNLocMaxN->SetXTitle("E (GeV)");
1936   outputContainer->Add(fhEventPlaneEtaNLocMaxN) ;
1937
1938   if(fFillEbinHisto)
1939   {
1940     for(Int_t i = 0; i < 4; i++)
1941     {
1942       fhMassM02NLocMax1Ebin[i]  = new TH2F(Form("hMassM02NLocMax1Ebin%d",i),
1943                                            Form("Invariant mass of split clusters vs #lambda_{0}^{2}, NLM=1, E bin %d",i),
1944                                            ssbins,ssmin,ssmax,mbins,mmin,mmax);
1945       fhMassM02NLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1946       fhMassM02NLocMax1Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1947       outputContainer->Add(fhMassM02NLocMax1Ebin[i]) ;
1948
1949       fhMassM02NLocMax2Ebin[i]  = new TH2F(Form("hMassM02NLocMax2Ebin%d",i),
1950                                            Form("Invariant mass of split clusters vs #lambda_{0}^{2}, NLM=2, E bin %d",i),
1951                                            ssbins,ssmin,ssmax,mbins,mmin,mmax);
1952       fhMassM02NLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1953       fhMassM02NLocMax2Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1954       outputContainer->Add(fhMassM02NLocMax2Ebin[i]) ;
1955
1956       fhMassM02NLocMaxNEbin[i]  = new TH2F(Form("hMassM02NLocMaxNEbin%d",i),
1957                                            Form("Invariant mass of split clusters vs vs #lambda_{0}^{2}, NLM>2, E bin %d",i),
1958                                            ssbins,ssmin,ssmax,mbins,mmin,mmax);
1959       fhMassM02NLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1960       fhMassM02NLocMaxNEbin[i]->SetXTitle("#lambda_{0}^{2}");
1961       outputContainer->Add(fhMassM02NLocMaxNEbin[i]) ;
1962
1963
1964       fhMassAsyNLocMax1Ebin[i]  = new TH2F(Form("hMassAsyNLocMax1Ebin%d",i),
1965                                            Form("Invariant mass of split clusters vs split asymmetry, NLM=1, E bin %d",i),
1966                                            200,-1,1,mbins,mmin,mmax);
1967       fhMassAsyNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1968       fhMassAsyNLocMax1Ebin[i]->SetXTitle("asymmetry");
1969       outputContainer->Add(fhMassAsyNLocMax1Ebin[i]) ;
1970
1971       fhMassAsyNLocMax2Ebin[i]  = new TH2F(Form("hMassAsyNLocMax2Ebin%d",i),
1972                                            Form("Invariant mass of split clusters vs split asymmetry, NLM=2, E bin %d",i),
1973                                            200,-1,1,mbins,mmin,mmax);
1974       fhMassAsyNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1975       fhMassAsyNLocMax2Ebin[i]->SetXTitle("asymmetry");
1976       outputContainer->Add(fhMassAsyNLocMax2Ebin[i]) ;
1977
1978       fhMassAsyNLocMaxNEbin[i]  = new TH2F(Form("hMassAsyNLocMaxNEbin%d",i),
1979                                            Form("Invariant mass of split clusters vs split asymmetry, NLM>2, E bin %d",i),
1980                                            200,-1,1,mbins,mmin,mmax);
1981       fhMassAsyNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1982       fhMassAsyNLocMaxNEbin[i]->SetXTitle("asymmetry");
1983       outputContainer->Add(fhMassAsyNLocMaxNEbin[i]) ;
1984
1985
1986       if(IsDataMC() && fFillMCHisto)
1987       {
1988         fhMCAsymM02NLocMax1MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax1MCPi0Ebin%d",i),
1989                                                     Form("Asymmetry of MC #pi^{0} vs #lambda_{0}^{2}, NLM=1, E bin %d",i),
1990                                                     ssbins,ssmin,ssmax,100,0,1);
1991         fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
1992         fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1993         outputContainer->Add(fhMCAsymM02NLocMax1MCPi0Ebin[i]) ;
1994
1995         fhMCAsymM02NLocMax2MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax2MCPi0Ebin%d",i),
1996                                                     Form("Asymmetry of MC #pi^{0} vs #lambda_{0}^{2}, NLM=2, E bin %d",i),
1997                                                     ssbins,ssmin,ssmax,100,0,1);
1998         fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
1999         fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
2000         outputContainer->Add(fhMCAsymM02NLocMax2MCPi0Ebin[i]) ;
2001
2002         fhMCAsymM02NLocMaxNMCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMaxNMCPi0Ebin%d",i),
2003                                                     Form("Asymmetry of MC #pi^{0} vs #lambda_{0}^{2}, NLM>2, E bin %d",i),
2004                                                     ssbins,ssmin,ssmax,100,0,1);
2005         fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetYTitle("Decay asymmetry");
2006         fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
2007         outputContainer->Add(fhMCAsymM02NLocMaxNMCPi0Ebin[i]) ;
2008
2009
2010         fhAsyMCGenRecoNLocMax1EbinPi0[i]  = new TH2F(Form("hAsyMCGenRecoNLocMax1Ebin%dPi0",i),
2011                                                      Form("Generated vs reconstructed asymmetry of split clusters from pi0, NLM=1, E bin %d",i),
2012                                                      200,-1,1,200,-1,1);
2013         fhAsyMCGenRecoNLocMax1EbinPi0[i]->SetYTitle("M (GeV/c^{2})");
2014         fhAsyMCGenRecoNLocMax1EbinPi0[i]->SetXTitle("asymmetry");
2015         outputContainer->Add(fhAsyMCGenRecoNLocMax1EbinPi0[i]) ;
2016
2017         fhAsyMCGenRecoNLocMax2EbinPi0[i]  = new TH2F(Form("hAsyMCGenRecoNLocMax2Ebin%dPi0",i),
2018                                                      Form("Generated vs reconstructed asymmetry of split clusters from pi0, NLM=2, E bin %d",i),
2019                                                      200,-1,1,200,-1,1);
2020         fhAsyMCGenRecoNLocMax2EbinPi0[i]->SetYTitle("M (GeV/c^{2})");
2021         fhAsyMCGenRecoNLocMax2EbinPi0[i]->SetXTitle("asymmetry");
2022         outputContainer->Add(fhAsyMCGenRecoNLocMax2EbinPi0[i]) ;
2023
2024         fhAsyMCGenRecoNLocMaxNEbinPi0[i]  = new TH2F(Form("hAsyMCGenRecoNLocMaxNEbin%dPi0",i),
2025                                                      Form("Generated vs reconstructed asymmetry of split clusters from pi0, NLM>2, E bin %d",i),
2026                                                      200,-1,1,200,-1,1);
2027         fhAsyMCGenRecoNLocMaxNEbinPi0[i]->SetYTitle("M (GeV/c^{2})");
2028         fhAsyMCGenRecoNLocMaxNEbinPi0[i]->SetXTitle("asymmetry");
2029         outputContainer->Add(fhAsyMCGenRecoNLocMaxNEbinPi0[i]) ;
2030       }
2031
2032       if(fFillSSExtraHisto)
2033       {
2034         fhMassDispEtaNLocMax1Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax1Ebin%d",i),
2035                                                  Form("Invariant mass of 2 highest energy cells #sigma_{#eta #eta}^{2}, E bin %d",i),
2036                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax);
2037         fhMassDispEtaNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
2038         fhMassDispEtaNLocMax1Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
2039         outputContainer->Add(fhMassDispEtaNLocMax1Ebin[i]) ;
2040
2041         fhMassDispEtaNLocMax2Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax2Ebin%d",i),
2042                                                  Form("Invariant mass of 2 local maxima cells #sigma_{#eta #eta}^{2}, E bin %d",i),
2043                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax);
2044         fhMassDispEtaNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
2045         fhMassDispEtaNLocMax2Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
2046         outputContainer->Add(fhMassDispEtaNLocMax2Ebin[i]) ;
2047
2048         fhMassDispEtaNLocMaxNEbin[i]  = new TH2F(Form("hMassDispEtaNLocMaxNEbin%d",i),
2049                                                  Form("Invariant mass of N>2 local maxima cells vs #sigma_{#eta #eta}^{2}, E bin %d",i),
2050                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax);
2051         fhMassDispEtaNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
2052         fhMassDispEtaNLocMaxNEbin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
2053         outputContainer->Add(fhMassDispEtaNLocMaxNEbin[i]) ;
2054
2055         fhMassDispPhiNLocMax1Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax1Ebin%d",i),
2056                                                  Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E bin %d",i),
2057                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax);
2058         fhMassDispPhiNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
2059         fhMassDispPhiNLocMax1Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
2060         outputContainer->Add(fhMassDispPhiNLocMax1Ebin[i]) ;
2061
2062         fhMassDispPhiNLocMax2Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax2Ebin%d",i),
2063                                                  Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E bin %d",i),
2064                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax);
2065         fhMassDispPhiNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
2066         fhMassDispPhiNLocMax2Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
2067         outputContainer->Add(fhMassDispPhiNLocMax2Ebin[i]) ;
2068
2069         fhMassDispPhiNLocMaxNEbin[i]  = new TH2F(Form("hMassDispPhiNLocMaxNEbin%d",i),
2070                                                  Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, E bin %d",i),
2071                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax);
2072         fhMassDispPhiNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
2073         fhMassDispPhiNLocMaxNEbin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
2074         outputContainer->Add(fhMassDispPhiNLocMaxNEbin[i]) ;
2075
2076         fhMassDispAsyNLocMax1Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax1Ebin%d",i),
2077                                                  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),
2078                                                  200,-1,1,mbins,mmin,mmax);
2079         fhMassDispAsyNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
2080         fhMassDispAsyNLocMax1Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
2081         outputContainer->Add(fhMassDispAsyNLocMax1Ebin[i]) ;
2082
2083         fhMassDispAsyNLocMax2Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax2Ebin%d",i),
2084                                                  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),
2085                                                  200,-1,1,mbins,mmin,mmax);
2086         fhMassDispAsyNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
2087         fhMassDispAsyNLocMax2Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
2088         outputContainer->Add(fhMassDispAsyNLocMax2Ebin[i]) ;
2089
2090         fhMassDispAsyNLocMaxNEbin[i]  = new TH2F(Form("hMassDispAsyNLocMaxNEbin%d",i),
2091                                                  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),
2092                                                  200,-1,1,mbins,mmin,mmax);
2093         fhMassDispAsyNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
2094         fhMassDispAsyNLocMaxNEbin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
2095         outputContainer->Add(fhMassDispAsyNLocMaxNEbin[i]) ;
2096       }
2097     }
2098   }
2099
2100   if(IsDataMC() && fFillMCHisto)
2101   {
2102     fhMCGenSplitEFracAfterCutsNLocMax1MCPi0     = new TH2F("hMCGenSplitEFracAfterCutsNLocMax1MCPi0",
2103                                                            "E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 1 MC Pi0, after M02 and Asym cut",
2104                                                            nptbins,ptmin,ptmax,200,0,2);
2105     fhMCGenSplitEFracAfterCutsNLocMax1MCPi0   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
2106     fhMCGenSplitEFracAfterCutsNLocMax1MCPi0   ->SetXTitle("E (GeV)");
2107     outputContainer->Add(fhMCGenSplitEFracAfterCutsNLocMax1MCPi0) ;
2108
2109     fhMCGenSplitEFracAfterCutsNLocMax2MCPi0    = new TH2F("hMCGenSplitEFracAfterCutsNLocMax2MCPi0",
2110                                                           "E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 2 MC Pi0, after M02 and Asym cut",
2111                                                           nptbins,ptmin,ptmax,200,0,2);
2112     fhMCGenSplitEFracAfterCutsNLocMax2MCPi0  ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
2113     fhMCGenSplitEFracAfterCutsNLocMax2MCPi0  ->SetXTitle("E (GeV)");
2114     outputContainer->Add(fhMCGenSplitEFracAfterCutsNLocMax2MCPi0) ;
2115
2116
2117     fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0    = new TH2F("hMCGenSplitEFracAfterCutsNLocMaxNMCPi0",
2118                                                           "E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  > 2 MC Pi0, after M02 and Asym cut",
2119                                                           nptbins,ptmin,ptmax,200,0,2);
2120     fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0  ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
2121     fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0  ->SetXTitle("E (GeV)");
2122     outputContainer->Add(fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0) ;
2123
2124     fhMCGenFracAfterCutsNLocMax1MCPi0     = new TH2F("hMCGenFracAfterCutsNLocMax1MCPi0",
2125                                                      "E_{gen} / E_{reco} vs E_{reco} for N max  = 1 MC Pi0, after M02 and Asym cut",
2126                                                      nptbins,ptmin,ptmax,200,0,2);
2127     fhMCGenFracAfterCutsNLocMax1MCPi0   ->SetYTitle("E_{gen} / E_{reco}");
2128     fhMCGenFracAfterCutsNLocMax1MCPi0   ->SetXTitle("E (GeV)");
2129     outputContainer->Add(fhMCGenFracAfterCutsNLocMax1MCPi0) ;
2130
2131     fhMCGenFracAfterCutsNLocMax2MCPi0    = new TH2F("hMCGenFracAfterCutsNLocMax2MCPi0",
2132                                                     " E_{gen} / E_{reco} vs E_{reco} for N max  = 2 MC Pi0, after M02 and Asym cut",
2133                                                     nptbins,ptmin,ptmax,200,0,2);
2134     fhMCGenFracAfterCutsNLocMax2MCPi0   ->SetYTitle("E_{gen} / E_{reco}");
2135     fhMCGenFracAfterCutsNLocMax2MCPi0   ->SetXTitle("E (GeV)");
2136     outputContainer->Add(fhMCGenFracAfterCutsNLocMax2MCPi0) ;
2137
2138
2139     fhMCGenFracAfterCutsNLocMaxNMCPi0   = new TH2F("hMCGenFracAfterCutsNLocMaxNMCPi0",
2140                                                    " E_{gen} / E_{reco}  vs E_{reco} for N max  > 2 MC Pi0, after M02 and Asym cut",
2141                                                    nptbins,ptmin,ptmax,200,0,2);
2142     fhMCGenFracAfterCutsNLocMaxNMCPi0   ->SetYTitle("E_{gen} / E_{reco}");
2143     fhMCGenFracAfterCutsNLocMaxNMCPi0   ->SetXTitle("E (GeV)");
2144     outputContainer->Add(fhMCGenFracAfterCutsNLocMaxNMCPi0) ;
2145
2146   }
2147
2148   if(fFillTMResidualHisto && fFillTMHisto)
2149   {
2150     for(Int_t i = 0; i < n; i++)
2151     {
2152
2153       fhTrackMatchedDEtaNLocMax1[i]  = new TH2F
2154       (Form("hTrackMatchedDEtaNLocMax1%s",pname[i].Data()),
2155        Form("d#eta of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
2156        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2157       fhTrackMatchedDEtaNLocMax1[i]->SetYTitle("d#eta");
2158       fhTrackMatchedDEtaNLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
2159
2160       fhTrackMatchedDPhiNLocMax1[i]  = new TH2F
2161       (Form("hTrackMatchedDPhiNLocMax1%s",pname[i].Data()),
2162        Form("d#phi of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
2163        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2164       fhTrackMatchedDPhiNLocMax1[i]->SetYTitle("d#phi (rad)");
2165       fhTrackMatchedDPhiNLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
2166
2167       outputContainer->Add(fhTrackMatchedDEtaNLocMax1[i]) ;
2168       outputContainer->Add(fhTrackMatchedDPhiNLocMax1[i]) ;
2169
2170       fhTrackMatchedDEtaNLocMax2[i]  = new TH2F
2171       (Form("hTrackMatchedDEtaNLocMax2%s",pname[i].Data()),
2172        Form("d#eta of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
2173        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2174       fhTrackMatchedDEtaNLocMax2[i]->SetYTitle("d#eta");
2175       fhTrackMatchedDEtaNLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
2176
2177       fhTrackMatchedDPhiNLocMax2[i]  = new TH2F
2178       (Form("hTrackMatchedDPhiNLocMax2%s",pname[i].Data()),
2179        Form("d#phi of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
2180        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2181       fhTrackMatchedDPhiNLocMax2[i]->SetYTitle("d#phi (rad)");
2182       fhTrackMatchedDPhiNLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
2183
2184       outputContainer->Add(fhTrackMatchedDEtaNLocMax2[i]) ;
2185       outputContainer->Add(fhTrackMatchedDPhiNLocMax2[i]) ;
2186
2187       fhTrackMatchedDEtaNLocMaxN[i]  = new TH2F
2188       (Form("hTrackMatchedDEtaNLocMaxN%s",pname[i].Data()),
2189        Form("d#eta of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
2190        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2191       fhTrackMatchedDEtaNLocMaxN[i]->SetYTitle("d#eta");
2192       fhTrackMatchedDEtaNLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
2193
2194       fhTrackMatchedDPhiNLocMaxN[i]  = new TH2F
2195       (Form("hTrackMatchedDPhiNLocMaxN%s",pname[i].Data()),
2196        Form("d#phi of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
2197        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2198       fhTrackMatchedDPhiNLocMaxN[i]->SetYTitle("d#phi (rad)");
2199       fhTrackMatchedDPhiNLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
2200
2201       outputContainer->Add(fhTrackMatchedDEtaNLocMaxN[i]) ;
2202       outputContainer->Add(fhTrackMatchedDPhiNLocMaxN[i]) ;
2203
2204       fhTrackMatchedDEtaNLocMax1Pos[i]  = new TH2F
2205       (Form("hTrackMatchedDEtaNLocMax1Pos%s",pname[i].Data()),
2206        Form("d#eta of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
2207        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2208       fhTrackMatchedDEtaNLocMax1Pos[i]->SetYTitle("d#eta");
2209       fhTrackMatchedDEtaNLocMax1Pos[i]->SetXTitle("E_{cluster} (GeV)");
2210
2211       fhTrackMatchedDPhiNLocMax1Pos[i]  = new TH2F
2212       (Form("hTrackMatchedDPhiNLocMax1Pos%s",pname[i].Data()),
2213        Form("d#phi of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
2214        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2215       fhTrackMatchedDPhiNLocMax1Pos[i]->SetYTitle("d#phi (rad)");
2216       fhTrackMatchedDPhiNLocMax1Pos[i]->SetXTitle("E_{cluster} (GeV)");
2217
2218       outputContainer->Add(fhTrackMatchedDEtaNLocMax1Pos[i]) ;
2219       outputContainer->Add(fhTrackMatchedDPhiNLocMax1Pos[i]) ;
2220
2221       fhTrackMatchedDEtaNLocMax2Pos[i]  = new TH2F
2222       (Form("hTrackMatchedDEtaNLocMax2Pos%s",pname[i].Data()),
2223        Form("d#eta of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
2224        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2225       fhTrackMatchedDEtaNLocMax2Pos[i]->SetYTitle("d#eta");
2226       fhTrackMatchedDEtaNLocMax2Pos[i]->SetXTitle("E_{cluster} (GeV)");
2227
2228       fhTrackMatchedDPhiNLocMax2Pos[i]  = new TH2F
2229       (Form("hTrackMatchedDPhiNLocMax2Pos%s",pname[i].Data()),
2230        Form("d#phi of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
2231        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2232       fhTrackMatchedDPhiNLocMax2Pos[i]->SetYTitle("d#phi (rad)");
2233       fhTrackMatchedDPhiNLocMax2Pos[i]->SetXTitle("E_{cluster} (GeV)");
2234
2235       outputContainer->Add(fhTrackMatchedDEtaNLocMax2Pos[i]) ;
2236       outputContainer->Add(fhTrackMatchedDPhiNLocMax2Pos[i]) ;
2237
2238       fhTrackMatchedDEtaNLocMaxNPos[i]  = new TH2F
2239       (Form("hTrackMatchedDEtaNLocMaxNPos%s",pname[i].Data()),
2240        Form("d#eta of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
2241        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2242       fhTrackMatchedDEtaNLocMaxNPos[i]->SetYTitle("d#eta");
2243       fhTrackMatchedDEtaNLocMaxNPos[i]->SetXTitle("E_{cluster} (GeV)");
2244
2245       fhTrackMatchedDPhiNLocMaxNPos[i]  = new TH2F
2246       (Form("hTrackMatchedDPhiNLocMaxNPos%s",pname[i].Data()),
2247        Form("d#phi of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
2248        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2249       fhTrackMatchedDPhiNLocMaxNPos[i]->SetYTitle("d#phi (rad)");
2250       fhTrackMatchedDPhiNLocMaxNPos[i]->SetXTitle("E_{cluster} (GeV)");
2251
2252       outputContainer->Add(fhTrackMatchedDEtaNLocMaxNPos[i]) ;
2253       outputContainer->Add(fhTrackMatchedDPhiNLocMaxNPos[i]) ;
2254
2255       fhTrackMatchedDEtaNLocMax1Neg[i]  = new TH2F
2256       (Form("hTrackMatchedDEtaNLocMax1Neg%s",pname[i].Data()),
2257        Form("d#eta of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
2258        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2259       fhTrackMatchedDEtaNLocMax1Neg[i]->SetYTitle("d#eta");
2260       fhTrackMatchedDEtaNLocMax1Neg[i]->SetXTitle("E_{cluster} (GeV)");
2261
2262       fhTrackMatchedDPhiNLocMax1Neg[i]  = new TH2F
2263       (Form("hTrackMatchedDPhiNLocMax1Neg%s",pname[i].Data()),
2264        Form("d#phi of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
2265        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2266       fhTrackMatchedDPhiNLocMax1Neg[i]->SetYTitle("d#phi (rad)");
2267       fhTrackMatchedDPhiNLocMax1Neg[i]->SetXTitle("E_{cluster} (GeV)");
2268
2269       outputContainer->Add(fhTrackMatchedDEtaNLocMax1Neg[i]) ;
2270       outputContainer->Add(fhTrackMatchedDPhiNLocMax1Neg[i]) ;
2271
2272       fhTrackMatchedDEtaNLocMax2Neg[i]  = new TH2F
2273       (Form("hTrackMatchedDEtaNLocMax2Neg%s",pname[i].Data()),
2274        Form("d#eta of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
2275        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2276       fhTrackMatchedDEtaNLocMax2Neg[i]->SetYTitle("d#eta");
2277       fhTrackMatchedDEtaNLocMax2Neg[i]->SetXTitle("E_{cluster} (GeV)");
2278
2279       fhTrackMatchedDPhiNLocMax2Neg[i]  = new TH2F
2280       (Form("hTrackMatchedDPhiNLocMax2Neg%s",pname[i].Data()),
2281        Form("d#phi of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
2282        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2283       fhTrackMatchedDPhiNLocMax2Neg[i]->SetYTitle("d#phi (rad)");
2284       fhTrackMatchedDPhiNLocMax2Neg[i]->SetXTitle("E_{cluster} (GeV)");
2285
2286       outputContainer->Add(fhTrackMatchedDEtaNLocMax2Neg[i]) ;
2287       outputContainer->Add(fhTrackMatchedDPhiNLocMax2Neg[i]) ;
2288
2289       fhTrackMatchedDEtaNLocMaxNNeg[i]  = new TH2F
2290       (Form("hTrackMatchedDEtaNLocMaxNNeg%s",pname[i].Data()),
2291        Form("d#eta of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
2292        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2293       fhTrackMatchedDEtaNLocMaxNNeg[i]->SetYTitle("d#eta");
2294       fhTrackMatchedDEtaNLocMaxNNeg[i]->SetXTitle("E_{cluster} (GeV)");
2295
2296       fhTrackMatchedDPhiNLocMaxNNeg[i]  = new TH2F
2297       (Form("hTrackMatchedDPhiNLocMaxNNeg%s",pname[i].Data()),
2298        Form("d#phi of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
2299        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2300       fhTrackMatchedDPhiNLocMaxNNeg[i]->SetYTitle("d#phi (rad)");
2301       fhTrackMatchedDPhiNLocMaxNNeg[i]->SetXTitle("E_{cluster} (GeV)");
2302
2303       outputContainer->Add(fhTrackMatchedDEtaNLocMaxNNeg[i]) ;
2304       outputContainer->Add(fhTrackMatchedDPhiNLocMaxNNeg[i]) ;
2305
2306     }
2307   }
2308
2309   if(fFillAngleHisto)
2310   {
2311     for(Int_t j = 0; j < nMatched; j++)
2312     {
2313
2314       fhAnglePairNLocMax1[j]  = new TH2F(Form("hAnglePairNLocMax1%s",sMatched[j].Data()),
2315                                         Form("Opening angle of 2 highest energy cells vs pair Energy, %s",sMatched[j].Data()),
2316                                         nptbins,ptmin,ptmax,200,0,0.2);
2317       fhAnglePairNLocMax1[j]->SetYTitle("#alpha (rad)");
2318       fhAnglePairNLocMax1[j]->SetXTitle("E (GeV)");
2319       outputContainer->Add(fhAnglePairNLocMax1[j]) ;
2320
2321       fhAnglePairNLocMax2[j]  = new TH2F(Form("hAnglePairNLocMax2%s",sMatched[j].Data()),
2322                                         Form("Opening angle of 2 local maxima cells vs Energy, %s",sMatched[j].Data()),
2323                                         nptbins,ptmin,ptmax,200,0,0.2);
2324       fhAnglePairNLocMax2[j]->SetYTitle("#alpha (rad)");
2325       fhAnglePairNLocMax2[j]->SetXTitle("E (GeV)");
2326       outputContainer->Add(fhAnglePairNLocMax2[j]) ;
2327
2328       fhAnglePairNLocMaxN[j]  = new TH2F(Form("hAnglePairNLocMaxN%s",sMatched[j].Data()),
2329                                         Form("Opening angle of N>2 local maxima cells vs Energy, %s",sMatched[j].Data()),
2330                                         nptbins,ptmin,ptmax,200,0,0.2);
2331       fhAnglePairNLocMaxN[j]->SetYTitle("#alpha (rad)");
2332       fhAnglePairNLocMaxN[j]->SetXTitle("E (GeV)");
2333       outputContainer->Add(fhAnglePairNLocMaxN[j]) ;
2334
2335       fhAnglePairMassNLocMax1[j]  = new TH2F(Form("hAnglePairMassNLocMax1%s",sMatched[j].Data()),
2336                                             Form("Opening angle of 2 highest energy cells vs Mass for E > 12 GeV, %s",sMatched[j].Data()),
2337                                             mbins,mmin,mmax,200,0,0.2);
2338       fhAnglePairMassNLocMax1[j]->SetXTitle("M (GeV/c^{2})");
2339       fhAnglePairMassNLocMax1[j]->SetYTitle("#alpha (rad)");
2340       outputContainer->Add(fhAnglePairMassNLocMax1[j]) ;
2341
2342       fhAnglePairMassNLocMax2[j]  = new TH2F(Form("hAnglePairMassNLocMax2%s",sMatched[j].Data()),
2343                                             Form("Opening angle of 2 local maxima cells vs Mass for E > 12 GeV, %s",sMatched[j].Data()),
2344                                             mbins,mmin,mmax,200,0,0.2);
2345       fhAnglePairMassNLocMax2[j]->SetXTitle("M (GeV/c^{2})");
2346       fhAnglePairMassNLocMax2[j]->SetYTitle("#alpha (rad)");
2347       outputContainer->Add(fhAnglePairMassNLocMax2[j]) ;
2348
2349       fhAnglePairMassNLocMaxN[j]  = new TH2F(Form("hAnglePairMassNLocMaxN%s",sMatched[j].Data()),
2350                                             Form("Opening angle of N>2 local maxima cells vs Mass for E > 12 GeV, %s",sMatched[j].Data()),
2351                                             mbins,mmin,mmax,200,0,0.2);
2352       fhAnglePairMassNLocMaxN[j]->SetXTitle("M (GeV/c^{2})");
2353       fhAnglePairMassNLocMaxN[j]->SetYTitle("#alpha (rad)");
2354       outputContainer->Add(fhAnglePairMassNLocMaxN[j]) ;
2355
2356     }
2357   }
2358
2359   for(Int_t j = 0; j < nMatched; j++)
2360   {
2361     fhSplitEFractionvsAsyNLocMax1[j]     = new TH2F(Form("hSplitEFractionvsAsyNLocMax1%s",sMatched[j].Data()),
2362                                                     Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  = 1, E>12, %s",sMatched[j].Data()),
2363                                                     100,-1,1,120,0,1.2);
2364     fhSplitEFractionvsAsyNLocMax1[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
2365     fhSplitEFractionvsAsyNLocMax1[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
2366     outputContainer->Add(fhSplitEFractionvsAsyNLocMax1[j]) ;
2367
2368     fhSplitEFractionvsAsyNLocMax2[j]     = new TH2F(Form("hSplitEFractionvsAsyNLocMax2%s",sMatched[j].Data()),
2369                                                     Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  = 2,E>12, %s",sMatched[j].Data()),
2370                                                     100,-1,1,120,0,1.2);
2371     fhSplitEFractionvsAsyNLocMax2[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
2372     fhSplitEFractionvsAsyNLocMax2[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
2373     outputContainer->Add(fhSplitEFractionvsAsyNLocMax2[j]) ;
2374
2375     fhSplitEFractionvsAsyNLocMaxN[j]    = new TH2F(Form("hSplitEFractionvsAsyNLocMaxN%s",sMatched[j].Data()),
2376                                                    Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  > 2, E>12, %s",sMatched[j].Data()),
2377                                                    100,-1,1,120,0,1.2);
2378     fhSplitEFractionvsAsyNLocMaxN[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
2379     fhSplitEFractionvsAsyNLocMaxN[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
2380     outputContainer->Add(fhSplitEFractionvsAsyNLocMaxN[j]) ;
2381   }
2382
2383
2384   fhClusterEtaPhiNLocMax1  = new TH2F
2385   ("hClusterEtaPhiNLocMax1","Neutral Clusters with E > 8 GeV, NLM = 1: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
2386   fhClusterEtaPhiNLocMax1->SetYTitle("#phi (rad)");
2387   fhClusterEtaPhiNLocMax1->SetXTitle("#eta");
2388   outputContainer->Add(fhClusterEtaPhiNLocMax1) ;
2389
2390   fhClusterEtaPhiNLocMax2  = new TH2F
2391   ("hClusterEtaPhiNLocMax2","Neutral Clusters with E > 8 GeV, NLM = 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
2392   fhClusterEtaPhiNLocMax2->SetYTitle("#phi (rad)");
2393   fhClusterEtaPhiNLocMax2->SetXTitle("#eta");
2394   outputContainer->Add(fhClusterEtaPhiNLocMax2) ;
2395
2396   fhClusterEtaPhiNLocMaxN  = new TH2F
2397   ("hClusterEtaPhiNLocMaxN","Neutral Clusters with E > 8 GeV, NLM > 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
2398   fhClusterEtaPhiNLocMaxN->SetYTitle("#phi (rad)");
2399   fhClusterEtaPhiNLocMaxN->SetXTitle("#eta");
2400   outputContainer->Add(fhClusterEtaPhiNLocMaxN) ;
2401
2402   fhPi0EtaPhiNLocMax1  = new TH2F
2403   ("hPi0EtaPhiNLocMax1","Selected #pi^{0}'s with E > 8 GeV, NLM = 1: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
2404   fhPi0EtaPhiNLocMax1->SetYTitle("#phi (rad)");
2405   fhPi0EtaPhiNLocMax1->SetXTitle("#eta");
2406   outputContainer->Add(fhPi0EtaPhiNLocMax1) ;
2407
2408   fhPi0EtaPhiNLocMax2  = new TH2F
2409   ("hPi0EtaPhiNLocMax2","Selected #pi^{0}'s with E > 8 GeV, NLM = 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
2410   fhPi0EtaPhiNLocMax2->SetYTitle("#phi (rad)");
2411   fhPi0EtaPhiNLocMax2->SetXTitle("#eta");
2412   outputContainer->Add(fhPi0EtaPhiNLocMax2) ;
2413
2414   fhPi0EtaPhiNLocMaxN  = new TH2F
2415   ("hPi0EtaPhiNLocMaxN","Selected #pi^{0}'s with E > 8 GeV, NLM > 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
2416   fhPi0EtaPhiNLocMaxN->SetYTitle("#phi (rad)");
2417   fhPi0EtaPhiNLocMaxN->SetXTitle("#eta");
2418   outputContainer->Add(fhPi0EtaPhiNLocMaxN) ;
2419
2420   fhEtaEtaPhiNLocMax1  = new TH2F
2421   ("hEtaEtaPhiNLocMax1","Selected #eta's with E > 8 GeV, NLM = 1: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
2422   fhEtaEtaPhiNLocMax1->SetYTitle("#phi (rad)");
2423   fhEtaEtaPhiNLocMax1->SetXTitle("#eta");
2424   outputContainer->Add(fhEtaEtaPhiNLocMax1) ;
2425
2426   fhEtaEtaPhiNLocMax2  = new TH2F
2427   ("hEtaEtaPhiNLocMax2","Selected #eta's with E > 8 GeV, NLM = 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
2428   fhEtaEtaPhiNLocMax2->SetYTitle("#phi (rad)");
2429   fhEtaEtaPhiNLocMax2->SetXTitle("#eta");
2430   outputContainer->Add(fhEtaEtaPhiNLocMax2) ;
2431
2432   fhEtaEtaPhiNLocMaxN  = new TH2F
2433   ("hEtaEtaPhiNLocMaxN","Selected #eta's with E > 8 GeV, NLM > 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
2434   fhEtaEtaPhiNLocMaxN->SetYTitle("#phi (rad)");
2435   fhEtaEtaPhiNLocMaxN->SetXTitle("#eta");
2436   outputContainer->Add(fhEtaEtaPhiNLocMaxN) ;
2437
2438   TString snlm[] = {"1","2","N"};
2439
2440   if(fFillSSWeightHisto)
2441   {
2442     for(Int_t nlm = 0; nlm < 3; nlm++)
2443     {
2444       fhPi0CellE[nlm]  = new TH2F(Form("hPi0CellENLocMax%s",snlm[nlm].Data()),
2445                                   Form("Selected #pi^{0}'s, NLM = %s: cluster E vs cell E",snlm[nlm].Data()),
2446                                   nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
2447       fhPi0CellE[nlm]->SetYTitle("E_{cell}");
2448       fhPi0CellE[nlm]->SetXTitle("E_{cluster}");
2449       outputContainer->Add(fhPi0CellE[nlm]) ;
2450
2451       fhPi0CellEFrac[nlm]  = new TH2F(Form("hPi0CellEFracNLocMax%s",snlm[nlm].Data()),
2452                                       Form("Selected #pi^{0}'s, NLM = %s: cluster E vs cell E / cluster E",snlm[nlm].Data()),
2453                                       nptbins,ptmin,ptmax, 100,0,1);
2454       fhPi0CellEFrac[nlm]->SetYTitle("E_{cell} / E_{cluster}");
2455       fhPi0CellEFrac[nlm]->SetXTitle("E_{cluster}");
2456       outputContainer->Add(fhPi0CellEFrac[nlm]) ;
2457
2458       fhPi0CellLogEFrac[nlm]  = new TH2F(Form("hPi0CellLogEFracNLocMax%s",snlm[nlm].Data()),
2459                                       Form("Selected #pi^{0}'s, NLM = %s: cluster E vs Log(cell E / cluster E)",snlm[nlm].Data()),
2460                                       nptbins,ptmin,ptmax, 100,-10,0);
2461       fhPi0CellLogEFrac[nlm]->SetYTitle("Log(E_{cell} / E_{cluster})");
2462       fhPi0CellLogEFrac[nlm]->SetXTitle("E_{cluster}");
2463       outputContainer->Add(fhPi0CellLogEFrac[nlm]) ;
2464
2465
2466       fhPi0CellEMaxEMax2Frac[nlm]  = new TH2F(Form("hPi0CellEMaxEMax2FracNLocMax%s",snlm[nlm].Data()),
2467                                       Form("Selected #pi^{0}'s, NLM = %s: cluster E vs 2nd loc. max. E / 1st loc. max.  E",snlm[nlm].Data()),
2468                                       nptbins,ptmin,ptmax, 100,0,1);
2469       fhPi0CellEMaxEMax2Frac[nlm]->SetYTitle("E_{Loc Max 2} / E_{Loc Max 1}");
2470       fhPi0CellEMaxEMax2Frac[nlm]->SetXTitle("E_{cluster}");
2471       outputContainer->Add(fhPi0CellEMaxEMax2Frac[nlm]) ;
2472
2473       fhPi0CellEMaxClusterFrac[nlm]  = new TH2F(Form("hPi0CellEMaxClusterFracNLocMax%s",snlm[nlm].Data()),
2474                                               Form("Selected #pi^{0}'s, NLM = %s: cluster E vs 1st loc. max. E / E cluster",snlm[nlm].Data()),
2475                                               nptbins,ptmin,ptmax, 100,0,1);
2476       fhPi0CellEMaxClusterFrac[nlm]->SetYTitle("E_{Loc Max 1} / E_{cluster}");
2477       fhPi0CellEMaxClusterFrac[nlm]->SetXTitle("E_{cluster}");
2478       outputContainer->Add(fhPi0CellEMaxClusterFrac[nlm]) ;
2479
2480       fhPi0CellEMax2ClusterFrac[nlm]  = new TH2F(Form("hPi0CellEMax2ClusterFracNLocMax%s",snlm[nlm].Data()),
2481                                                 Form("Selected #pi^{0}'s, NLM = %s: cluster E vs 2nd loc. max. E / E cluster",snlm[nlm].Data()),
2482                                                 nptbins,ptmin,ptmax, 100,0,1);
2483       fhPi0CellEMax2ClusterFrac[nlm]->SetYTitle("E_{Loc Max 2} / E_{cluster}");
2484       fhPi0CellEMax2ClusterFrac[nlm]->SetXTitle("E_{cluster}");
2485       outputContainer->Add(fhPi0CellEMax2ClusterFrac[nlm]) ;
2486
2487       fhPi0CellEMaxFrac[nlm]  = new TH2F(Form("hPi0CellEMaxFracNLocMax%s",snlm[nlm].Data()),
2488                                                 Form("Selected #pi^{0}'s, NLM = %s: cluster E vs 1st loc. max. E / E cell i",snlm[nlm].Data()),
2489                                                 nptbins,ptmin,ptmax, 100,0,1);
2490       fhPi0CellEMaxFrac[nlm]->SetYTitle("E_{Loc Max 1} / E_{cell i}");
2491       fhPi0CellEMaxFrac[nlm]->SetXTitle("E_{cluster}");
2492       outputContainer->Add(fhPi0CellEMaxFrac[nlm]) ;
2493
2494       fhPi0CellEMax2Frac[nlm]  = new TH2F(Form("hPi0CellEMax2FracNLocMax%s",snlm[nlm].Data()),
2495                                                  Form("Selected #pi^{0}'s, NLM = %s: cluster E vs 2nd loc. max. E / E cell i",snlm[nlm].Data()),
2496                                                  nptbins,ptmin,ptmax, 200,0,2);
2497       fhPi0CellEMax2Frac[nlm]->SetYTitle("E_{Loc Max 2} / E_{cell i}");
2498       fhPi0CellEMax2Frac[nlm]->SetXTitle("E_{cluster}");
2499       outputContainer->Add(fhPi0CellEMax2Frac[nlm]) ;
2500
2501
2502       for(Int_t i = 0; i < fSSWeightN; i++)
2503       {
2504         fhM02WeightPi0[nlm][i]     = new TH2F(Form("hM02Pi0NLocMax%s_W%d",snlm[nlm].Data(),i),
2505                                               Form("#lambda_{0}^{2} vs E, with W0 = %2.2f, for N Local max = %s", fSSWeight[i], snlm[nlm].Data()),
2506                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2507         fhM02WeightPi0[nlm][i]   ->SetYTitle("#lambda_{0}^{2}");
2508         fhM02WeightPi0[nlm][i]   ->SetXTitle("E (GeV)");
2509         outputContainer->Add(fhM02WeightPi0[nlm][i]) ;
2510       }
2511
2512       for(Int_t i = 0; i < fSSECellCutN; i++)
2513       {
2514         fhM02ECellCutPi0[nlm][i]     = new TH2F(Form("hM02Pi0NLocMax%s_Ecell%d",snlm[nlm].Data(),i),
2515                                               Form("#lambda_{0}^{2} vs E, with Ecell > %2.2f, for N Local max = %s", fSSECellCut[i], snlm[nlm].Data()),
2516                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2517         fhM02ECellCutPi0[nlm][i]   ->SetYTitle("#lambda_{0}^{2}");
2518         fhM02ECellCutPi0[nlm][i]   ->SetXTitle("E (GeV)");
2519         outputContainer->Add(fhM02ECellCutPi0[nlm][i]) ;
2520       }
2521
2522     }
2523   }
2524
2525   Int_t tdbins   = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax  = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin  = GetHistogramRanges()->GetHistoDiffTimeMin();
2526
2527   fhPi0EPairDiffTimeNLM1 = new TH2F("hPi0EPairDiffTimeNLocMax1","cluster pair time difference vs E, selected #pi, NLM=1",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2528   fhPi0EPairDiffTimeNLM1->SetXTitle("E_{pair} (GeV)");
2529   fhPi0EPairDiffTimeNLM1->SetYTitle("#Delta t (ns)");
2530   outputContainer->Add(fhPi0EPairDiffTimeNLM1);
2531
2532   fhPi0EPairDiffTimeNLM2 = new TH2F("hPi0EPairDiffTimeNLocMax2","cluster pair time difference vs E, selected #pi, NLM=2",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2533   fhPi0EPairDiffTimeNLM2->SetXTitle("E_{pair} (GeV)");
2534   fhPi0EPairDiffTimeNLM2->SetYTitle("#Delta t (ns)");
2535   outputContainer->Add(fhPi0EPairDiffTimeNLM2);
2536
2537   fhPi0EPairDiffTimeNLMN = new TH2F("hPi0EPairDiffTimeNLocMaxN","cluster pair time difference vs E, selected #pi, NLM>2",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2538   fhPi0EPairDiffTimeNLMN->SetXTitle("E_{pair} (GeV)");
2539   fhPi0EPairDiffTimeNLMN->SetYTitle("#Delta t (ns)");
2540   outputContainer->Add(fhPi0EPairDiffTimeNLMN);
2541
2542   fhEtaEPairDiffTimeNLM1 = new TH2F("hEtaEPairDiffTimeNLocMax1","cluster pair time difference vs E, selected #eta, NLM=1",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2543   fhEtaEPairDiffTimeNLM1->SetXTitle("E_{pair} (GeV)");
2544   fhEtaEPairDiffTimeNLM1->SetYTitle("#Delta t (ns)");
2545   outputContainer->Add(fhEtaEPairDiffTimeNLM1);
2546
2547   fhEtaEPairDiffTimeNLM2 = new TH2F("hEtaEPairDiffTimeNLocMax2","cluster pair time difference vs E, selected #eta, NLM=2",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2548   fhEtaEPairDiffTimeNLM2->SetXTitle("E_{pair} (GeV)");
2549   fhEtaEPairDiffTimeNLM2->SetYTitle("#Delta t (ns)");
2550   outputContainer->Add(fhEtaEPairDiffTimeNLM2);
2551
2552   fhEtaEPairDiffTimeNLMN = new TH2F("hEtaEPairDiffTimeNLocMaxN","cluster pair time difference vs E, selected #eta, NLM>2",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2553   fhEtaEPairDiffTimeNLMN->SetXTitle("E_{pair} (GeV)");
2554   fhEtaEPairDiffTimeNLMN->SetYTitle("#Delta t (ns)");
2555   outputContainer->Add(fhEtaEPairDiffTimeNLMN);
2556
2557
2558   if(IsDataMC() && fFillMCOverlapHisto)
2559   {
2560     for(Int_t i = 1; i < n; i++)
2561     {
2562       for(Int_t j = 0; j < 3; j++)
2563       {
2564         fhMCENOverlaps[j][i]     = new TH2F(Form("hMCENOverlapsNLocMax%s%s",snlm[j].Data(),pname[i].Data()),
2565                                             Form("# overlaps vs E for NLM=%s, %s",snlm[j].Data(),ptype[i].Data()),
2566                                             nptbins,ptmin,ptmax,10,0,10);
2567         fhMCENOverlaps[j][i]   ->SetYTitle("# overlaps");
2568         fhMCENOverlaps[j][i]   ->SetXTitle("E (GeV)");
2569         outputContainer->Add(fhMCENOverlaps[j][i]) ;
2570
2571         fhMCEM02Overlap0[j][i]     = new TH2F(Form("hMCEM02Overlap0NLocMax%s%s",snlm[j].Data(),pname[i].Data()),
2572                                               Form("Overlap 0, #lambda_{0}^{2} vs E for NLM=%s, %s",snlm[j].Data(),ptype[i].Data()),
2573                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2574         fhMCEM02Overlap0[j][i]   ->SetYTitle("#lambda_{0}^{2}");
2575         fhMCEM02Overlap0[j][i]   ->SetXTitle("E (GeV)");
2576         outputContainer->Add(fhMCEM02Overlap0[j][i]) ;
2577
2578         fhMCEM02Overlap1[j][i]     = new TH2F(Form("hMCEM02Overlap1NLocMax%s%s",snlm[j].Data(), pname[i].Data()),
2579                                               Form("Overlap 1, #lambda_{0}^{2} vs E for NLM=%s, %s",snlm[j].Data(),ptype[i].Data()),
2580                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2581         fhMCEM02Overlap1[j][i]   ->SetYTitle("#lambda_{0}^{2}");
2582         fhMCEM02Overlap1[j][i]   ->SetXTitle("E (GeV)");
2583         outputContainer->Add(fhMCEM02Overlap1[j][i]) ;
2584
2585         fhMCEM02OverlapN[j][i]     = new TH2F(Form("hMCEM02OverlapNNLocMax%s%s",snlm[j].Data(), pname[i].Data()),
2586                                               Form("Overlap N, #lambda_{0}^{2} vs E for NLM=%s %s",snlm[j].Data(),ptype[i].Data()),
2587                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2588         fhMCEM02OverlapN[j][i]   ->SetYTitle("#lambda_{0}^{2}");
2589         fhMCEM02OverlapN[j][i]   ->SetXTitle("E (GeV)");
2590         outputContainer->Add(fhMCEM02OverlapN[j][i]) ;
2591
2592         fhMCEMassOverlap0[j][i]     = new TH2F(Form("hMCEMassOverlap0NLocMax%s%s",snlm[j].Data(),pname[i].Data()),
2593                                                Form("Overlap 0, Mass vs E for NLM=%s, %s",snlm[j].Data(),ptype[i].Data()),
2594                                                nptbins,ptmin,ptmax,mbins,mmin,mmax);
2595         fhMCEMassOverlap0[j][i]   ->SetYTitle("Mass (GeV/c^{2}");
2596         fhMCEMassOverlap0[j][i]   ->SetXTitle("E (GeV)");
2597         outputContainer->Add(fhMCEMassOverlap0[j][i]) ;
2598
2599         fhMCEMassOverlap1[j][i]     = new TH2F(Form("hMCEMassOverlap1NLocMax%s%s",snlm[j].Data(), pname[i].Data()),
2600                                                Form("Overalap 1, Mass vs E for NLM=%s, %s",snlm[j].Data(),ptype[i].Data()),
2601                                                nptbins,ptmin,ptmax,mbins,mmin,mmax);
2602         fhMCEMassOverlap1[j][i]   ->SetYTitle("Mass (GeV/c^{2}");
2603         fhMCEMassOverlap1[j][i]   ->SetXTitle("E (GeV)");
2604         outputContainer->Add(fhMCEMassOverlap1[j][i]) ;
2605
2606         fhMCEMassOverlapN[j][i]     = new TH2F(Form("hMCEMassOverlapNNLocMax%s%s",snlm[j].Data(), pname[i].Data()),
2607                                                Form("Overlap N, Mass vs E for NLM=%s %s",snlm[j].Data(),ptype[i].Data()),
2608                                                nptbins,ptmin,ptmax,mbins,mmin,mmax);
2609         fhMCEMassOverlapN[j][i]   ->SetYTitle("Mass (GeV/c^{2})");
2610         fhMCEMassOverlapN[j][i]   ->SetXTitle("E (GeV)");
2611         outputContainer->Add(fhMCEMassOverlapN[j][i]) ;
2612
2613
2614         if(i < 5)
2615         {
2616           fhMCPi0MassM02Overlap0[j][i-1]  = new TH2F(Form("hMCPi0MassM02Overlap0NLocMax%sEbin%d",snlm[j].Data(),i-1),
2617                                                    Form("Overlap 0, Mass vs #lambda_{0}^{2}, NLM=%s, E bin %d",snlm[j].Data(),i-1),
2618                                                    ssbins,ssmin,ssmax,mbins,mmin,mmax);
2619           fhMCPi0MassM02Overlap0[j][i-1]->SetYTitle("M (GeV/c^{2})");
2620           fhMCPi0MassM02Overlap0[j][i-1]->SetXTitle("#lambda_{0}^{2}");
2621           outputContainer->Add(fhMCPi0MassM02Overlap0[j][i-1]) ;
2622
2623           fhMCPi0MassM02Overlap1[j][i-1]  = new TH2F(Form("hMCPi0MassM02Overlap1NLocMax%sEbin%d",snlm[j].Data(),i-1),
2624                                                    Form("Overlap 1, Mass vs #lambda_{0}^{2}, NLM=%s, E bin %d",snlm[j].Data(),i-1),
2625                                                    ssbins,ssmin,ssmax,mbins,mmin,mmax);
2626           fhMCPi0MassM02Overlap1[j][i-1]->SetYTitle("M (GeV/c^{2})");
2627           fhMCPi0MassM02Overlap1[j][i-1]->SetXTitle("#lambda_{0}^{2}");
2628           outputContainer->Add(fhMCPi0MassM02Overlap1[j][i-1]) ;
2629
2630           fhMCPi0MassM02OverlapN[j][i-1]  = new TH2F(Form("hMCPi0MassM02OverlapNNLocMax%sEbin%d",snlm[j].Data(),i-1),
2631                                                    Form("Overlap N, Mass vs #lambda_{0}^{2}, NLM=%s, E bin %d",snlm[j].Data(),i-1),
2632                                                    ssbins,ssmin,ssmax,mbins,mmin,mmax);
2633           fhMCPi0MassM02OverlapN[j][i-1]->SetYTitle("M (GeV/c^{2})");
2634           fhMCPi0MassM02OverlapN[j][i-1]->SetXTitle("#lambda_{0}^{2}");
2635           outputContainer->Add(fhMCPi0MassM02OverlapN[j][i-1]) ;
2636         }
2637
2638         if(fFillTMHisto)
2639         {
2640           fhMCENOverlapsMatch[j][i]     = new TH2F(Form("hMCENOverlapsNLocMax%s%sMatched",snlm[j].Data(),pname[i].Data()),
2641                                                    Form("# overlaps vs E for NLM=%s, %s",snlm[j].Data(),ptype[i].Data()),
2642                                                    nptbins,ptmin,ptmax,10,0,10);
2643           fhMCENOverlapsMatch[j][i]   ->SetYTitle("# overlaps");
2644           fhMCENOverlapsMatch[j][i]   ->SetXTitle("E (GeV)");
2645           outputContainer->Add(fhMCENOverlapsMatch[j][i]) ;
2646
2647           fhMCEM02Overlap0Match[j][i]     = new TH2F(Form("hMCEM02Overlap0NLocMax%s%sMatched",snlm[j].Data(),pname[i].Data()),
2648                                                      Form("#lambda_{0}^{2} vs E for NLM=%s, %s, Track Matched",snlm[j].Data(),ptype[i].Data()),
2649                                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2650           fhMCEM02Overlap0Match[j][i]   ->SetYTitle("#lambda_{0}^{2}");
2651           fhMCEM02Overlap0Match[j][i]   ->SetXTitle("E (GeV)");
2652           outputContainer->Add(fhMCEM02Overlap0Match[j][i]) ;
2653
2654           fhMCEM02Overlap1Match[j][i]     = new TH2F(Form("hMCEM02Overlap1NLocMax%s%sMatched",snlm[j].Data(), pname[i].Data()),
2655                                                      Form("#lambda_{0}^{2} vs E for NLM=%s, %s, Track Matched",snlm[j].Data(),ptype[i].Data()),
2656                                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2657           fhMCEM02Overlap1Match[j][i]   ->SetYTitle("#lambda_{0}^{2}");
2658           fhMCEM02Overlap1Match[j][i]   ->SetXTitle("E (GeV)");
2659           outputContainer->Add(fhMCEM02Overlap1Match[j][i]) ;
2660
2661           fhMCEM02OverlapNMatch[j][i]     = new TH2F(Form("hMCEM02OverlapNNLocMax%s%sMatched",snlm[j].Data(), pname[i].Data()),
2662                                                      Form("#lambda_{0}^{2} vs E for NLM=%s, %s, Track Matched",snlm[j].Data(),ptype[i].Data()),
2663                                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2664           fhMCEM02OverlapNMatch[j][i]   ->SetYTitle("#lambda_{0}^{2}");
2665           fhMCEM02OverlapNMatch[j][i]   ->SetXTitle("E (GeV)");
2666           outputContainer->Add(fhMCEM02OverlapNMatch[j][i]) ;
2667
2668           fhMCEMassOverlap0Match[j][i]     = new TH2F(Form("hMCEMassOverlap0NLocMax%s%sMatched",snlm[j].Data(),pname[i].Data()),
2669                                                       Form("Mass vs E for NLM=%s, %s, Track Matched",snlm[j].Data(),ptype[i].Data()),
2670                                                       nptbins,ptmin,ptmax,mbins,mmin,mmax);
2671           fhMCEMassOverlap0Match[j][i]   ->SetYTitle("Mass (GeV/c^{2}");
2672           fhMCEMassOverlap0Match[j][i]   ->SetXTitle("E (GeV)");
2673           outputContainer->Add(fhMCEMassOverlap0Match[j][i]) ;
2674
2675           fhMCEMassOverlap1Match[j][i]     = new TH2F(Form("hMCEMassOverlap1NLocMax%s%sMatched",snlm[j].Data(), pname[i].Data()),
2676                                                       Form("Mass vs E for NLM=%s, %s, Track Matched",snlm[j].Data(),ptype[i].Data()),
2677                                                       nptbins,ptmin,ptmax,mbins,mmin,mmax);
2678           fhMCEMassOverlap1Match[j][i]   ->SetYTitle("Mass (GeV/c^{2}");
2679           fhMCEMassOverlap1Match[j][i]   ->SetXTitle("E (GeV)");
2680           outputContainer->Add(fhMCEMassOverlap1Match[j][i]) ;
2681
2682           fhMCEMassOverlapNMatch[j][i]     = new TH2F(Form("hMCEMassOverlapNNLocMax%s%sMatched",snlm[j].Data(), pname[i].Data()),
2683                                                       Form("Mass vs E for NLM=%s, %s, Track Matched",snlm[j].Data(),ptype[i].Data()),
2684                                                       nptbins,ptmin,ptmax,mbins,mmin,mmax);
2685           fhMCEMassOverlapNMatch[j][i]   ->SetYTitle("Mass (GeV/c^{2}");
2686           fhMCEMassOverlapNMatch[j][i]   ->SetXTitle("E (GeV)");
2687           outputContainer->Add(fhMCEMassOverlapNMatch[j][i]) ;
2688
2689           if(i < 5)
2690           {
2691             fhMCPi0MassM02Overlap0Match[j][i-1]  = new TH2F(Form("hMCPi0MassM02Overlap0NLocMax%sEbin%dMatched",snlm[j].Data(),i-1),
2692                                                      Form("Overlap 0, Mass vs #lambda_{0}^{2}, NLM=%s, E bin %d, Track Matched",snlm[j].Data(),i-1),
2693                                                      ssbins,ssmin,ssmax,mbins,mmin,mmax);
2694             fhMCPi0MassM02Overlap0Match[j][i-1]->SetYTitle("M (GeV/c^{2})");
2695             fhMCPi0MassM02Overlap0Match[j][i-1]->SetXTitle("#lambda_{0}^{2}");
2696             outputContainer->Add(fhMCPi0MassM02Overlap0Match[j][i-1]) ;
2697
2698             fhMCPi0MassM02Overlap1Match[j][i-1]  = new TH2F(Form("hMCPi0MassM02Overlap1NLocMax%sEbin%dMatched",snlm[j].Data(),i-1),
2699                                                      Form("Overlap 1, Mass vs #lambda_{0}^{2}, NLM=%s, E bin %d, Track Matched",snlm[j].Data(),i-1),
2700                                                      ssbins,ssmin,ssmax,mbins,mmin,mmax);
2701             fhMCPi0MassM02Overlap1Match[j][i-1]->SetYTitle("M (GeV/c^{2})");
2702             fhMCPi0MassM02Overlap1Match[j][i-1]->SetXTitle("#lambda_{0}^{2}");
2703             outputContainer->Add(fhMCPi0MassM02Overlap1Match[j][i-1]) ;
2704
2705             fhMCPi0MassM02OverlapNMatch[j][i-1]  = new TH2F(Form("hMCPi0MassM02OverlapNNLocMax%sEbin%dMatched",snlm[j].Data(),i-1),
2706                                                      Form("Overlap N, Mass vs #lambda_{0}^{2}, NLM=%s, E bin %d, Track Matched",snlm[j].Data(),i-1),
2707                                                      ssbins,ssmin,ssmax,mbins,mmin,mmax);
2708             fhMCPi0MassM02OverlapNMatch[j][i-1]->SetYTitle("M (GeV/c^{2})");
2709             fhMCPi0MassM02OverlapNMatch[j][i-1]->SetXTitle("#lambda_{0}^{2}");
2710             outputContainer->Add(fhMCPi0MassM02OverlapNMatch[j][i-1]) ;
2711
2712           }
2713
2714         }
2715       }
2716     }
2717
2718     fhMCPi0HighNLMPair    = new TH2F("hMCPi0HighNLMPair","NLM vs E for merged pi0 cluster, high energy NLM pair are decays",
2719                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2720     fhMCPi0HighNLMPair   ->SetYTitle("N maxima");
2721     fhMCPi0HighNLMPair   ->SetXTitle("E (GeV)");
2722     outputContainer->Add(fhMCPi0HighNLMPair) ;
2723
2724     fhMCPi0LowNLMPair     = new TH2F("hMCPi0LowNLMPair","NLM vs E for merged pi0 cluster, lower energy NLM pair are decays",
2725                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2726     fhMCPi0LowNLMPair   ->SetYTitle("N maxima");
2727     fhMCPi0LowNLMPair   ->SetXTitle("E (GeV)");
2728     outputContainer->Add(fhMCPi0LowNLMPair) ;
2729
2730     fhMCPi0AnyNLMPair     = new TH2F("hMCPi0AnyNLMPair","NLM vs E for merged pi0 cluster, both high and other energy NLM pair are decays",
2731                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2732     fhMCPi0AnyNLMPair   ->SetYTitle("N maxima");
2733     fhMCPi0AnyNLMPair   ->SetXTitle("E (GeV)");
2734     outputContainer->Add(fhMCPi0AnyNLMPair) ;
2735
2736     fhMCPi0NoneNLMPair     = new TH2F("hMCPi0NoneNLMPair","NLM vs E for merged pi0 cluster, no NLM pair are decays",
2737                                       nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2738     fhMCPi0NoneNLMPair   ->SetYTitle("N maxima");
2739     fhMCPi0NoneNLMPair   ->SetXTitle("E (GeV)");
2740     outputContainer->Add(fhMCPi0NoneNLMPair) ;
2741
2742
2743     fhMCPi0HighNLMPairNoMCMatch    = new TH2F("hMCPi0HighNLMPairNoMCMatch","NLM vs E for merged pi0 cluster, high energy NLM pair are decays",
2744                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2745     fhMCPi0HighNLMPairNoMCMatch   ->SetYTitle("N maxima");
2746     fhMCPi0HighNLMPairNoMCMatch   ->SetXTitle("E (GeV)");
2747     outputContainer->Add(fhMCPi0HighNLMPairNoMCMatch) ;
2748
2749     fhMCPi0LowNLMPairNoMCMatch     = new TH2F("hMCPi0LowNLMPairNoMCMatch","NLM vs E for merged pi0 cluster, lower energy NLM pair are decays",
2750                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2751     fhMCPi0LowNLMPairNoMCMatch   ->SetYTitle("N maxima");
2752     fhMCPi0LowNLMPairNoMCMatch   ->SetXTitle("E (GeV)");
2753     outputContainer->Add(fhMCPi0LowNLMPairNoMCMatch) ;
2754
2755     fhMCPi0AnyNLMPairNoMCMatch     = new TH2F("hMCPi0AnyNLMPairNoMCMatch","NLM vs E for merged pi0 cluster, both high and other energy NLM pair are decays",
2756                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2757     fhMCPi0AnyNLMPairNoMCMatch   ->SetYTitle("N maxima");
2758     fhMCPi0AnyNLMPairNoMCMatch   ->SetXTitle("E (GeV)");
2759     outputContainer->Add(fhMCPi0AnyNLMPairNoMCMatch) ;
2760
2761     fhMCPi0NoneNLMPairNoMCMatch     = new TH2F("hMCPi0NoneNLMPairNoMCMatch","NLM vs E for merged pi0 cluster, no NLM pair are decays",
2762                                       nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2763     fhMCPi0NoneNLMPairNoMCMatch   ->SetYTitle("N maxima");
2764     fhMCPi0NoneNLMPairNoMCMatch   ->SetXTitle("E (GeV)");
2765     outputContainer->Add(fhMCPi0NoneNLMPairNoMCMatch) ;
2766
2767
2768     fhMCEOverlapType = new TH2F("hMCEOverlapType","Kind of overlap particle, neutral clusters",
2769                                 nptbins,ptmin,ptmax,5,0,5);
2770     //fhMCEOverlapType   ->SetYTitle("Overlap Type");
2771     fhMCEOverlapType->GetYaxis()->SetBinLabel(1 ,"#gamma");
2772     fhMCEOverlapType->GetYaxis()->SetBinLabel(2 ,"e^{#pm}");
2773     fhMCEOverlapType->GetYaxis()->SetBinLabel(3 ,"hadron^{#pm}");
2774     fhMCEOverlapType->GetYaxis()->SetBinLabel(4 ,"hadron^{0}");
2775     fhMCEOverlapType->GetYaxis()->SetBinLabel(5 ,"??");
2776     fhMCEOverlapType   ->SetXTitle("Cluster E (GeV)");
2777     outputContainer->Add(fhMCEOverlapType) ;
2778
2779     fhMCEOverlapTypeMatch = new TH2F("hMCEOverlapTypeMatched","Kind of overlap particle, charged clusters",
2780                                 nptbins,ptmin,ptmax,5,0,5);
2781     //fhMCEOverlapTypeMatch   ->SetYTitle("Overlap Type");
2782     fhMCEOverlapTypeMatch->GetYaxis()->SetBinLabel(1 ,"#gamma");
2783     fhMCEOverlapTypeMatch->GetYaxis()->SetBinLabel(2 ,"e^{#pm}");
2784     fhMCEOverlapTypeMatch->GetYaxis()->SetBinLabel(3 ,"hadron^{#pm}");
2785     fhMCEOverlapTypeMatch->GetYaxis()->SetBinLabel(4 ,"hadron^{0}");
2786     fhMCEOverlapTypeMatch->GetYaxis()->SetBinLabel(5 ,"??");
2787     fhMCEOverlapTypeMatch->SetXTitle("Cluster E (GeV)");
2788     outputContainer->Add(fhMCEOverlapTypeMatch) ;
2789
2790   }// MC analysis, check overlaps
2791
2792   return outputContainer ;
2793
2794 }
2795
2796 //_____________________________________________________________________________
2797 void AliAnaInsideClusterInvariantMass::GetMCIndex(AliVCluster* cluster, Int_t & mcindex)
2798 {
2799
2800   // Assign mc index depending on MC bit set, to be used in histograms arrays
2801
2802   if(!IsDataMC()) return;
2803
2804   Int_t tag     = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader());
2805
2806   if      ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) &&
2807            !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPi0;
2808   else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )      mcindex = kmcPi0Conv;
2809   else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )      mcindex = kmcEta;
2810   else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2811            !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPhoton;
2812   else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
2813             GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
2814   else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))   mcindex = kmcElectron;
2815   else                                                                                mcindex = kmcHadron;
2816
2817 }
2818
2819
2820 //______________________________________________________________________________________________________________________
2821 void AliAnaInsideClusterInvariantMass::GetMCPrimaryKine(AliVCluster* cluster, const Int_t mcindex, const Bool_t matched,
2822                                                         Float_t & eprim, Float_t & asymGen, Int_t & noverlaps )
2823 {
2824   // Check origin of the candidates, get primary kinematics if overlapped meson decay
2825
2826   if(!IsDataMC()) return ;
2827
2828   Bool_t ok      = kFALSE;
2829   Int_t  mcLabel = cluster->GetLabel();
2830
2831   TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
2832   eprim = primary.E();
2833
2834   if(mcindex == kmcPi0 || mcindex == kmcEta || mcindex == kmcPi0Conv)
2835   {
2836     if(mcindex == kmcPi0 || mcindex == kmcPi0Conv)
2837     {
2838       asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,111,GetReader(),ok));
2839       TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok);
2840       if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
2841     }
2842     else
2843     {
2844       asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,221,GetReader(),ok));
2845       TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok);
2846       if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
2847     }
2848   }
2849
2850   if(!fFillMCOverlapHisto) return;
2851
2852   // Compare the primary depositing more energy with the rest,
2853   // if no photon/electron (conversion) or neutral meson as comon ancestor, consider it as other particle contributing
2854   Int_t ancPDG = 0, ancStatus = -1;
2855   TLorentzVector momentum; TVector3 prodVertex;
2856   Int_t ancLabel = 0;
2857   noverlaps = 0;
2858   for (UInt_t ilab = 1; ilab < cluster->GetNLabels(); ilab++ )
2859   {
2860     ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab],
2861                                                          GetReader(),ancPDG,ancStatus,momentum,prodVertex);
2862
2863     if(ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG!=111 && ancPDG!=221)
2864     {
2865       noverlaps++;
2866
2867       // What is the origin of the overlap?
2868       Bool_t  mOK = 0,      gOK = 0;
2869       Int_t   mpdg = -999999,  gpdg = -1;
2870       Int_t   mstatus = -1, gstatus = -1;
2871       Int_t gLabel = -1, ggLabel = -1;
2872       TLorentzVector mother      = GetMCAnalysisUtils()->GetMother     (cluster->GetLabels()[ilab],GetReader(),mpdg,mstatus,mOK);
2873       TLorentzVector grandmother = GetMCAnalysisUtils()->GetGrandMother(cluster->GetLabels()[ilab],GetReader(),gpdg,gstatus,gOK, gLabel,ggLabel);
2874
2875       //printf("Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
2876
2877       if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
2878           ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
2879             gLabel >=0 )
2880       {
2881         Int_t label = gLabel;
2882         while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
2883         {
2884           mpdg=gpdg;
2885           grandmother = GetMCAnalysisUtils()->GetGrandMother(label,GetReader(),gpdg,gstatus,ok, gLabel,ggLabel);
2886           label=gLabel;
2887         }
2888       }
2889
2890       //printf("; Final PDG %d\n",mpdg);
2891       Float_t histobin = -1;
2892       if     (mpdg==22)      histobin = 0.5;
2893       else if(TMath::Abs(mpdg)==11) histobin = 1.5;
2894       else if(mpdg==-999999) histobin = 4.5;
2895       else {
2896         Double_t charge = TDatabasePDG::Instance()->GetParticle(mpdg)->Charge();
2897         if(TMath::Abs(charge) > 0 ) histobin = 2.5;
2898         else                        histobin = 3.5;
2899         //printf("charge %f\n",charge);
2900       }
2901
2902       //printf("pdg = %d, histobin %2.1f\n",mpdg,histobin);
2903       if(histobin > 0)
2904       {
2905         if(matched)fhMCEOverlapType     ->Fill(cluster->E(),histobin);
2906         else       fhMCEOverlapTypeMatch->Fill(cluster->E(),histobin);
2907       }
2908     }
2909   }
2910
2911 }
2912
2913 //___________________________________________
2914 void AliAnaInsideClusterInvariantMass::Init()
2915 {
2916   //Init
2917   //Do some checks
2918   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD())
2919   {
2920     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
2921     abort();
2922   }
2923   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD())
2924   {
2925     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
2926     abort();
2927   }
2928
2929   if( GetReader()->GetDataType() == AliCaloTrackReader::kMC )
2930   {
2931     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
2932     abort();
2933
2934   }
2935
2936 }
2937
2938 //_____________________________________________________
2939 void AliAnaInsideClusterInvariantMass::InitParameters()
2940 {
2941   //Initialize the parameters of the analysis.
2942   AddToHistogramsName("AnaPi0InsideClusterInvariantMass_");
2943
2944   fCalorimeter = "EMCAL" ;
2945
2946   fM02MinCut   = 0.26 ;
2947   fM02MaxCut   = 10 ;
2948
2949   fMinNCells   = 4 ;
2950   fMinBadDist  = 2 ;
2951
2952   fHistoECut   = 8 ;
2953
2954   fSSWeightN   = 5;
2955   fSSWeight[0] = 4.6;  fSSWeight[1] = 4.7; fSSWeight[2] = 4.8; fSSWeight[3] = 4.9; fSSWeight[4] = 5.0;
2956   fSSWeight[5] = 5.1;  fSSWeight[6] = 5.2; fSSWeight[7] = 5.3; fSSWeight[8] = 5.4; fSSWeight[9] = 5.5;
2957
2958   fSSECellCutN   = 10;
2959   fSSECellCut[0] = 0.16;  fSSECellCut[1] = 0.18; fSSECellCut[2] = 0.2; fSSECellCut[3] = 0.22; fSSECellCut[4] = 0.24;
2960   fSSECellCut[5] = 0.26;  fSSECellCut[6] = 0.28; fSSECellCut[7] = 0.3; fSSECellCut[8] = 0.32; fSSECellCut[9] = 0.34;
2961
2962 }
2963
2964
2965 //__________________________________________________________________
2966 void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
2967 {
2968   //Search for pi0 in fCalorimeter with shower shape analysis
2969
2970   TObjArray * pl       = 0x0;
2971   AliVCaloCells* cells = 0x0;
2972
2973   //Select the Calorimeter of the photon
2974   if(fCalorimeter == "PHOS")
2975   {
2976     pl    = GetPHOSClusters();
2977     cells = GetPHOSCells();
2978   }
2979   else if (fCalorimeter == "EMCAL")
2980   {
2981     pl    = GetEMCALClusters();
2982     cells = GetEMCALCells();
2983   }
2984
2985   if(!pl || !cells)
2986   {
2987     Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2988     return;
2989   }
2990
2991         if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
2992
2993   for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++)
2994   {
2995     AliVCluster * cluster = (AliVCluster*) (pl->At(icluster));
2996
2997     //-------------------------------------------
2998     // Get cluster parameters, do some rejection
2999     //-------------------------------------------
3000
3001     Float_t en = cluster->E();
3002     Float_t l0 = cluster->GetM02();
3003     Int_t   nc = cluster->GetNCells();
3004     Float_t bd = cluster->GetDistanceToBadChannel() ;
3005
3006     //If too small or big E or low number of cells, or close to a bad channel skip it
3007
3008     if( en < GetMinEnergy() || en > GetMaxEnergy() || nc < fMinNCells || bd < fMinBadDist) continue ;
3009
3010     // Track-cluster matching
3011
3012     Bool_t  matched   = IsTrackMatched(cluster,GetReader()->GetInputEvent());
3013     if(!fFillTMHisto && matched) continue ;
3014
3015     // Get cluster angles
3016
3017     TLorentzVector lv;
3018     cluster->GetMomentum(lv, GetVertex(0));
3019     Float_t eta = lv.Eta();
3020     Float_t phi = lv.Phi();
3021     if(phi<0) phi=+TMath::TwoPi();
3022
3023     //printf("en %2.2f, GetMinEnergy() %2.2f, GetMaxEnergy() %2.2f, nc %d, fMinNCells %d,  bd %2.2f, fMinBadDist %2.2f\n",
3024     //       en,GetMinEnergy(), GetMaxEnergy(), nc, fMinNCells, bd, fMinBadDist);
3025
3026     // Get PID, N local maximum, *** split cluster ***
3027
3028     Int_t    nMax = 0;
3029     Double_t mass = 0., angle = 0.;
3030     TLorentzVector    l1, l2;
3031     Int_t    absId1 = -1; Int_t absId2 = -1;
3032
3033     Int_t pidTag = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(cluster,cells,GetCaloUtils(),
3034                                                                                GetVertex(0), nMax, mass, angle,
3035                                                                                l1,l2,absId1,absId2);
3036     if (nMax <= 0)
3037     {
3038       if(GetDebug() > 0 )
3039         printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found! It did not pass CaloPID selection criteria \n");
3040
3041       return;
3042     }
3043
3044     // Set some index for array histograms
3045
3046     Int_t inlm = -1;
3047     if     (nMax == 1) inlm = 0;
3048     else if(nMax == 2) inlm = 1;
3049     else if(nMax >  2) inlm = 2;
3050     else printf("Wrong N local maximum -> %d, n cells in cluster %d \n",nMax,nc);
3051
3052     // Get sub-cluster parameters
3053
3054     Float_t e1 = l1.Energy();
3055     Float_t e2 = l2.Energy();
3056
3057     Double_t tof1  = cells->GetCellTime(absId1);
3058     GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3059     tof1*=1.e9;
3060
3061     Double_t tof2  = cells->GetCellTime(absId2);
3062     GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3063     tof2*=1.e9;
3064
3065     Double_t t12diff = tof1-tof2;
3066
3067     Float_t splitFrac = (e1+e2)/en;
3068
3069     Float_t asym = -10;
3070     if(e1+e2>0) asym = (e1-e2)/(e1+e2);
3071
3072     //
3073
3074     Int_t ebin = -1;
3075     if(en > 8  && en <= 12) ebin = 0;
3076     if(en > 12 && en <= 16) ebin = 1;
3077     if(en > 16 && en <= 20) ebin = 2;
3078     if(en > 20)             ebin = 3;
3079
3080     // MC indexes
3081
3082     Int_t mcindex   = -1;
3083     //Int_t mcLabel   = cluster->GetLabel();
3084     GetMCIndex(cluster,mcindex);
3085
3086     // MC primary kine, generation fractions
3087
3088     Float_t eprim   = -1;
3089     Float_t asymGen = -2;
3090     Int_t noverlaps =  0;
3091     GetMCPrimaryKine(cluster,mcindex,matched,eprim,asymGen,noverlaps);
3092
3093     //
3094
3095     // For cluster with MC pi0 and more than 1 maxima
3096     CheckLocalMaximaMCOrigin(cluster, mcindex);
3097
3098     //----------------
3099     // Fill histograms
3100     //----------------
3101
3102     fhNLocMax[0][matched]->Fill(en,nMax);
3103     if(IsDataMC())
3104       fhNLocMax[mcindex][matched]->Fill(en,nMax);
3105
3106     if     ( nMax == 1  )
3107     {
3108       fhM02NLocMax1[0][matched]->Fill(en,l0) ;
3109       fhSplitEFractionNLocMax1[0][matched]->Fill(en,splitFrac) ;
3110
3111       if(IsDataMC())
3112       {
3113         fhM02NLocMax1[mcindex][matched]->Fill(en,l0) ;
3114         fhSplitEFractionNLocMax1[mcindex][matched]->Fill(en,splitFrac) ;
3115       }
3116
3117       if(en > fHistoECut)
3118       {
3119         fhMassM02NLocMax1[0][matched]->Fill(l0, mass);
3120         if( IsDataMC() ) fhMassM02NLocMax1[mcindex][matched]->Fill(l0, mass);
3121
3122         fhSplitEFractionvsAsyNLocMax1[matched]->Fill(asym,splitFrac) ;
3123         if(!matched)fhClusterEtaPhiNLocMax1->Fill(eta,phi);
3124       }
3125
3126
3127     }
3128     else if( nMax == 2  )
3129     {
3130       fhM02NLocMax2[0][matched]->Fill(en,l0) ;
3131       fhSplitEFractionNLocMax2[0][matched]->Fill(en,splitFrac) ;
3132
3133       if(IsDataMC())
3134       {
3135         fhM02NLocMax2[mcindex][matched]->Fill(en,l0) ;
3136         fhSplitEFractionNLocMax2[mcindex][matched]->Fill(en,splitFrac) ;
3137       }
3138
3139       if(en > fHistoECut)
3140       {
3141         fhMassM02NLocMax2[0][matched]->Fill(l0,  mass );
3142         if( IsDataMC() ) fhMassM02NLocMax2[mcindex][matched]->Fill(l0,mass);
3143
3144         fhSplitEFractionvsAsyNLocMax2[matched]->Fill(asym,splitFrac) ;
3145         if(!matched)fhClusterEtaPhiNLocMax2->Fill(eta,phi);
3146       }
3147     }
3148     else if( nMax >= 3  )
3149     {
3150       fhM02NLocMaxN[0][matched]->Fill(en,l0) ;
3151       fhSplitEFractionNLocMaxN[0][matched]->Fill(en,splitFrac) ;
3152
3153       if(IsDataMC())
3154       {
3155         fhM02NLocMaxN[mcindex][matched]->Fill(en,l0) ;
3156         fhSplitEFractionNLocMaxN[mcindex][matched]->Fill(en,splitFrac) ;
3157       }
3158
3159       if(en > fHistoECut)
3160       {
3161
3162         fhMassM02NLocMaxN[0][matched]->Fill(l0,mass);
3163         if( IsDataMC() ) fhMassM02NLocMaxN[mcindex][matched]->Fill(l0,mass);
3164
3165         fhSplitEFractionvsAsyNLocMaxN[matched]->Fill(asym,splitFrac) ;
3166         if(!matched)fhClusterEtaPhiNLocMaxN->Fill(eta,phi);
3167       }
3168     }
3169
3170     //
3171
3172     FillTrackMatchingHistograms(cluster,nMax,mcindex);
3173
3174     //
3175
3176     FillSSExtraHistograms(cluster, nMax, matched,mcindex,mass,ebin)  ;
3177
3178     //
3179
3180     FillMCHistograms(en,e1,e2,ebin,mcindex,l0,mass,
3181                      nMax,matched,splitFrac, asym, eprim,asymGen);
3182
3183     //
3184
3185     FillMCOverlapHistograms(en,mass,l0,inlm,ebin,matched,mcindex,noverlaps);
3186
3187     //
3188
3189     if(!matched) FillEBinHistograms(ebin,nMax,mcindex,splitFrac,mass,asym,l0);
3190
3191     //---------------------------------------------------------------------
3192     // From here only if M02 is large but not too large
3193     // or more strict cuts, fill histograms
3194     //---------------------------------------------------------------------
3195
3196     if( l0 < fM02MinCut || l0 > fM02MaxCut ) continue ;
3197
3198     FillAngleHistograms(matched,nMax,en,angle,mass);
3199
3200     Bool_t m02OK = GetCaloPID()->IsInPi0M02Range(en,l0,nMax);
3201     Bool_t asyOK = GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax);
3202
3203     Float_t efrac      = eprim/en;
3204     Float_t efracSplit = 0;
3205     if(e1+e2 > 0) efracSplit = eprim/(e1+e2);
3206
3207     Float_t cent = GetEventCentrality();
3208     Float_t evp  = GetEventPlaneAngle();
3209
3210     fhNLocMaxM02Cut[0][matched]->Fill(en,nMax);
3211     if(IsDataMC()) fhNLocMaxM02Cut[mcindex][matched]->Fill(en,nMax);
3212
3213     Float_t splitFracMin = GetCaloPID()->GetSplitEnergyFractionMinimum(inlm) ;
3214
3215     if     (nMax==1)
3216     {
3217       fhMassNLocMax1[0][matched]->Fill(en,mass );
3218       fhAsymNLocMax1[0][matched]->Fill(en,asym );
3219
3220       // Effect of cuts in mass histograms
3221
3222       if(!matched)
3223       {
3224         if(m02OK)
3225         {
3226           fhMassM02CutNLocMax1->Fill(en,mass);
3227           fhAsymM02CutNLocMax1->Fill(en,asym );
3228           if(splitFrac > splitFracMin) fhMassSplitECutNLocMax1->Fill(en,mass );
3229         } // m02
3230       } // split frac
3231
3232       if(m02OK && asyOK)
3233       {
3234         fhSplitEFractionAfterCutsNLocMax1[0][matched]->Fill(en,splitFrac);
3235         if(splitFrac > splitFracMin) fhMassAfterCutsNLocMax1[0][matched]->Fill(en,mass);
3236
3237         if(!matched && IsDataMC() && fFillMCHisto && mcindex==kmcPi0)
3238         {
3239           fhMCGenFracAfterCutsNLocMax1MCPi0      ->Fill(en   ,  efrac     );
3240           fhMCGenSplitEFracAfterCutsNLocMax1MCPi0->Fill(en   ,  efracSplit);
3241         }
3242       }
3243
3244       if     (pidTag==AliCaloPID::kPhoton)
3245       {
3246         fhM02ConNLocMax1 [0][matched]->Fill(en,l0);
3247         fhMassConNLocMax1[0][matched]->Fill(en,mass);
3248         fhAsyConNLocMax1 [0][matched]->Fill(en,asym);
3249       }
3250       else if(pidTag==AliCaloPID::kPi0   )
3251       {
3252         fhM02Pi0NLocMax1[0][matched]->Fill(en,l0);
3253         fhMassPi0NLocMax1[0][matched]->Fill(en,mass);
3254         fhAsyPi0NLocMax1[0][matched]->Fill(en,asym);
3255         fhCentralityPi0NLocMax1[0][matched]->Fill(en,cent) ;
3256
3257         if(!matched)
3258         {
3259           fhEventPlanePi0NLocMax1->Fill(en,evp) ;
3260           if(en > fHistoECut)fhPi0EtaPhiNLocMax1->Fill(eta,phi);
3261           FillSSWeightHistograms(cluster, 0, absId1, absId2);
3262           fhPi0EPairDiffTimeNLM1->Fill(e1+e2,t12diff);
3263         }
3264       }
3265       else if(pidTag==AliCaloPID::kEta)
3266       {
3267         fhM02EtaNLocMax1 [0][matched]->Fill(en,l0);
3268         fhMassEtaNLocMax1[0][matched]->Fill(en,mass);
3269         fhAsyEtaNLocMax1 [0][matched]->Fill(en,asym);
3270         fhCentralityEtaNLocMax1[0][matched]->Fill(en,cent) ;
3271
3272         if(!matched)
3273         {
3274           fhEventPlaneEtaNLocMax1->Fill(en,evp) ;
3275           if(en > fHistoECut)fhEtaEtaPhiNLocMax1->Fill(eta,phi);
3276           fhEtaEPairDiffTimeNLM1->Fill(e1+e2,t12diff);
3277         }
3278       }
3279
3280     }
3281     else if(nMax==2)
3282     {
3283       fhMassNLocMax2[0][matched]->Fill(en,mass );
3284       fhAsymNLocMax2[0][matched]->Fill(en,asym );
3285
3286       // Effect of cuts in mass histograms
3287
3288       if(!matched)
3289       {
3290         if(m02OK)
3291         {
3292           fhMassM02CutNLocMax2->Fill(en,mass);
3293           fhAsymM02CutNLocMax2->Fill(en,asym );
3294           if(splitFrac > splitFracMin) fhMassSplitECutNLocMax2->Fill(en,mass );
3295         } // m02
3296       } // split frac
3297
3298       if(m02OK && asyOK)
3299       {
3300         fhSplitEFractionAfterCutsNLocMax2[0][matched]->Fill(en,splitFrac);
3301         if(splitFrac >splitFracMin) fhMassAfterCutsNLocMax2[0][matched]->Fill(en,mass);
3302
3303         if(!matched && IsDataMC() && fFillMCHisto && mcindex==kmcPi0)
3304         {
3305
3306           fhMCGenFracAfterCutsNLocMax2MCPi0      ->Fill(en   ,  efrac     );
3307           fhMCGenSplitEFracAfterCutsNLocMax2MCPi0->Fill(en   ,  efracSplit);
3308         }
3309       }
3310
3311       if     (pidTag==AliCaloPID::kPhoton)
3312       {
3313         fhM02ConNLocMax2 [0][matched]->Fill(en,l0);
3314         fhMassConNLocMax2[0][matched]->Fill(en,mass);
3315         fhAsyConNLocMax2 [0][matched]->Fill(en,asym);
3316       }
3317       else if(pidTag==AliCaloPID::kPi0   )
3318       {
3319         fhM02Pi0NLocMax2 [0][matched]->Fill(en,l0);
3320         fhMassPi0NLocMax2[0][matched]->Fill(en,mass);
3321         fhAsyPi0NLocMax2 [0][matched]->Fill(en,asym);
3322         fhCentralityPi0NLocMax2[0][matched]->Fill(en,cent) ;
3323
3324         if(!matched)
3325         {
3326           fhEventPlanePi0NLocMax2->Fill(en,evp) ;
3327           if(en > fHistoECut)fhPi0EtaPhiNLocMax2->Fill(eta,phi);
3328           FillSSWeightHistograms(cluster, 1, absId1, absId2);
3329           fhPi0EPairDiffTimeNLM2->Fill(e1+e2,t12diff);
3330         }
3331       }
3332       else if(pidTag==AliCaloPID::kEta)
3333       {
3334         fhM02EtaNLocMax2 [0][matched]->Fill(en,l0);
3335         fhMassEtaNLocMax2[0][matched]->Fill(en,mass);
3336         fhAsyEtaNLocMax2 [0][matched]->Fill(en,asym);
3337         fhCentralityEtaNLocMax2[0][matched]->Fill(en,cent) ;
3338
3339         if(!matched)
3340         {
3341           fhEventPlaneEtaNLocMax2->Fill(en,evp) ;
3342           if(en > fHistoECut)fhEtaEtaPhiNLocMax2->Fill(eta,phi);
3343           fhEtaEPairDiffTimeNLM2->Fill(e1+e2,t12diff);
3344         }
3345       }
3346
3347     }
3348     else if(nMax >2)
3349     {
3350       fhMassNLocMaxN[0][matched]->Fill(en,mass);
3351       fhAsymNLocMaxN[0][matched]->Fill(en,asym);
3352
3353       // Effect of cuts in mass histograms
3354       if(!matched)
3355       {
3356         if(m02OK)
3357         {
3358           fhMassM02CutNLocMaxN->Fill(en,mass);
3359           fhAsymM02CutNLocMaxN->Fill(en,asym );
3360           if(splitFrac > splitFracMin)fhMassSplitECutNLocMaxN->Fill(en,mass );
3361         } // m02
3362       } // split frac
3363
3364       if(m02OK && asyOK)
3365       {
3366         fhSplitEFractionAfterCutsNLocMaxN[0][matched]->Fill(en,splitFrac);
3367         if(splitFrac > splitFracMin) fhMassAfterCutsNLocMaxN[0][matched]->Fill(en,mass);
3368
3369         if(!matched && IsDataMC() && fFillMCHisto && mcindex==kmcPi0)
3370         {
3371           fhMCGenFracAfterCutsNLocMaxNMCPi0      ->Fill(en   ,  efrac     );
3372           fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0->Fill(en   ,  efracSplit);
3373         }
3374       }
3375
3376       if     (pidTag==AliCaloPID::kPhoton)
3377       {
3378         fhM02ConNLocMaxN [0][matched]->Fill(en,l0);
3379         fhMassConNLocMaxN[0][matched]->Fill(en,mass);
3380         fhAsyConNLocMaxN [0][matched]->Fill(en,asym);
3381       }
3382       else if(pidTag==AliCaloPID::kPi0   )
3383       {
3384         fhM02Pi0NLocMaxN [0][matched]->Fill(en,l0);
3385         fhMassPi0NLocMaxN[0][matched]->Fill(en,mass);
3386         fhAsyPi0NLocMaxN [0][matched]->Fill(en,asym);
3387         fhCentralityPi0NLocMaxN[0][matched]->Fill(en,cent) ;
3388
3389         if(!matched)
3390         {
3391           fhEventPlanePi0NLocMaxN->Fill(en,evp) ;
3392           if(en > fHistoECut)fhPi0EtaPhiNLocMaxN->Fill(eta,phi);
3393           FillSSWeightHistograms(cluster, 2,  absId1, absId2);
3394           fhPi0EPairDiffTimeNLMN->Fill(e1+e2,t12diff);
3395         }
3396       }
3397       else if(pidTag==AliCaloPID::kEta)
3398       {
3399         fhM02EtaNLocMaxN [0][matched]->Fill(en,l0);
3400         fhMassEtaNLocMaxN[0][matched]->Fill(en,mass);
3401         fhAsyEtaNLocMaxN [0][matched]->Fill(en,asym);
3402         fhCentralityEtaNLocMaxN[0][matched]->Fill(en,cent) ;
3403
3404         if(!matched)
3405         {
3406           fhEventPlaneEtaNLocMaxN->Fill(en,evp) ;
3407           if(en > fHistoECut)fhEtaEtaPhiNLocMaxN->Fill(eta,phi);
3408           fhEtaEPairDiffTimeNLMN->Fill(e1+e2,t12diff);
3409         }
3410       }
3411     }
3412
3413     if(IsDataMC())
3414     {
3415       if     (nMax==1)
3416       {
3417         fhMassNLocMax1[mcindex][matched]->Fill(en,mass);
3418         fhAsymNLocMax1[mcindex][matched]->Fill(en,asym);
3419
3420         if(asyOK && m02OK)
3421         {
3422           fhSplitEFractionAfterCutsNLocMax1[mcindex][matched]->Fill(en,splitFrac);
3423           if(splitFrac > splitFracMin)
3424             fhMassAfterCutsNLocMax1[mcindex][matched]->Fill(en,mass);
3425         }
3426
3427         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMax1[mcindex][matched]->Fill(en,l0); fhMassConNLocMax1[mcindex][matched]->Fill(en,mass); fhAsyConNLocMax1[mcindex][matched]->Fill(en,asym); }
3428         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMax1[mcindex][matched]->Fill(en,l0); fhMassPi0NLocMax1[mcindex][matched]->Fill(en,mass); fhAsyPi0NLocMax1[mcindex][matched]->Fill(en,asym); }
3429         else if(pidTag==AliCaloPID::kEta   ) { fhM02EtaNLocMax1[mcindex][matched]->Fill(en,l0); fhMassEtaNLocMax1[mcindex][matched]->Fill(en,mass); fhAsyEtaNLocMax1[mcindex][matched]->Fill(en,asym); }
3430
3431         if     (pidTag==AliCaloPID::kPi0) fhCentralityPi0NLocMax1[mcindex][matched]->Fill(en,cent) ;
3432         else if(pidTag==AliCaloPID::kEta) fhCentralityEtaNLocMax1[mcindex][matched]->Fill(en,cent) ;
3433       }
3434       else if(nMax==2)
3435       {
3436         fhMassNLocMax2[mcindex][matched]->Fill(en,mass);
3437         fhAsymNLocMax2[mcindex][matched]->Fill(en,asym);
3438
3439         if(asyOK && m02OK)
3440         {
3441           fhSplitEFractionAfterCutsNLocMax2[mcindex][matched]->Fill(en,splitFrac);
3442           if(splitFrac >splitFracMin)
3443             fhMassAfterCutsNLocMax2[mcindex][matched]->Fill(en,mass);
3444         }
3445
3446         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMax2[mcindex][matched]->Fill(en,l0); fhMassConNLocMax2[mcindex][matched]->Fill(en,mass); fhAsyConNLocMax2[mcindex][matched]->Fill(en,asym); }
3447         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMax2[mcindex][matched]->Fill(en,l0); fhMassPi0NLocMax2[mcindex][matched]->Fill(en,mass); fhAsyPi0NLocMax2[mcindex][matched]->Fill(en,asym); }
3448         else if(pidTag==AliCaloPID::kEta   ) { fhM02EtaNLocMax2[mcindex][matched]->Fill(en,l0); fhMassEtaNLocMax2[mcindex][matched]->Fill(en,mass); fhAsyEtaNLocMax2[mcindex][matched]->Fill(en,asym); }
3449
3450         if     (pidTag==AliCaloPID::kPi0) fhCentralityPi0NLocMax2[mcindex][matched]->Fill(en,cent) ;
3451         else if(pidTag==AliCaloPID::kEta) fhCentralityEtaNLocMax2[mcindex][matched]->Fill(en,cent) ;
3452
3453       }
3454       else if(nMax >2)
3455       {
3456         fhMassNLocMaxN[mcindex][matched]->Fill(en,mass);
3457         fhAsymNLocMaxN[mcindex][matched]->Fill(en,asym);
3458
3459         if(asyOK && m02OK)
3460         {
3461           fhSplitEFractionAfterCutsNLocMaxN[mcindex][matched]->Fill(en,splitFrac);
3462           if(splitFrac > splitFracMin )
3463             fhMassAfterCutsNLocMaxN[mcindex][matched]->Fill(en,mass);
3464         }
3465
3466         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMaxN[mcindex][matched]->Fill(en,l0); fhMassConNLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyConNLocMaxN[mcindex][matched]->Fill(en,asym); }
3467         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMaxN[mcindex][matched]->Fill(en,l0); fhMassPi0NLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyPi0NLocMaxN[mcindex][matched]->Fill(en,asym); }
3468         else if(pidTag==AliCaloPID::kEta   ) { fhM02EtaNLocMaxN[mcindex][matched]->Fill(en,l0); fhMassEtaNLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyEtaNLocMaxN[mcindex][matched]->Fill(en,asym); }
3469
3470         if     (pidTag==AliCaloPID::kPi0) fhCentralityPi0NLocMaxN[mcindex][matched]->Fill(en,cent) ;
3471         else if(pidTag==AliCaloPID::kEta) fhCentralityEtaNLocMaxN[mcindex][matched]->Fill(en,cent) ;
3472       }
3473
3474     }//Work with MC truth first
3475
3476   }//loop
3477
3478   if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");
3479
3480 }
3481
3482 //______________________________________________________________________
3483 void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
3484 {
3485   //Print some relevant parameters set for the analysis
3486   if(! opt)
3487     return;
3488
3489   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3490   AliAnaCaloTrackCorrBaseClass::Print("");
3491   printf("Calorimeter     =     %s\n",  fCalorimeter.Data()) ;
3492   if(GetCaloUtils()) printf("Loc. Max. E > %2.2f\n",       GetCaloUtils()->GetLocalMaximaCutE());
3493   if(GetCaloUtils()) printf("Loc. Max. E Diff > %2.2f\n",  GetCaloUtils()->GetLocalMaximaCutEDiff());
3494   printf("Min. N Cells =%d \n",         fMinNCells) ;
3495   printf("Min. Dist. to Bad =%1.1f \n", fMinBadDist) ;
3496   printf("%2.2f < lambda_0^2 <%2.2f \n",fM02MinCut,fM02MaxCut);
3497   if(fFillSSWeightHisto) printf(" N w %d - N e cut %d \n",fSSWeightN,fSSECellCutN);
3498
3499   printf("    \n") ;
3500
3501 }
3502
3503
3504 //___________________________________________________________________________________________________________________
3505 void AliAnaInsideClusterInvariantMass::RecalculateClusterShowerShapeParametersWithCellCut(const AliEMCALGeometry * geom,
3506                                                                                           AliVCaloCells* cells,
3507                                                                                           AliVCluster * cluster,
3508                                                                                           Float_t & l0,   Float_t & l1,
3509                                                                                           Float_t & disp, Float_t & dEta, Float_t & dPhi,
3510                                                                                           Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi,
3511                                                                                           Float_t eCellMin)
3512 {
3513   // Calculates new center of gravity in the local EMCAL-module coordinates
3514   // and tranfers into global ALICE coordinates
3515   // Calculates Dispersion and main axis
3516
3517   if(!cluster)
3518   {
3519     AliInfo("Cluster pointer null!");
3520     return;
3521   }
3522
3523   Double_t eCell       = 0.;
3524   Float_t  fraction    = 1.;
3525   Float_t  recalFactor = 1.;
3526
3527   Int_t    iSupMod = -1;
3528   Int_t    iTower  = -1;
3529   Int_t    iIphi   = -1;
3530   Int_t    iIeta   = -1;
3531   Int_t    iphi    = -1;
3532   Int_t    ieta    = -1;
3533   Double_t etai    = -1.;
3534   Double_t phii    = -1.;
3535
3536   Int_t    nstat   = 0 ;
3537   Float_t  wtot    = 0.;
3538   Double_t w       = 0.;
3539   Double_t etaMean = 0.;
3540   Double_t phiMean = 0.;
3541
3542   //Loop on cells, calculate the cluster energy, in case a cut on cell energy is added
3543   // and to check if the cluster is between 2 SM in eta
3544   Int_t   iSM0   = -1;
3545   Bool_t  shared = kFALSE;
3546   Float_t energy = 0;
3547
3548   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
3549   {
3550     //Get from the absid the supermodule, tower and eta/phi numbers
3551     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
3552     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
3553
3554     //Check if there are cells of different SM
3555     if     (iDigit == 0   ) iSM0 = iSupMod;
3556     else if(iSupMod!= iSM0) shared = kTRUE;
3557
3558     //Get the cell energy, if recalibration is on, apply factors
3559     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
3560     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
3561
3562     if(GetCaloUtils()->GetEMCALRecoUtils()->IsRecalibrationOn())
3563     {
3564       recalFactor = GetCaloUtils()->GetEMCALRecoUtils()->GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
3565     }
3566
3567     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
3568
3569     if(eCell > eCellMin) energy += eCell;
3570
3571   }//cell loop
3572
3573   //Loop on cells, get weighted parameters
3574   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
3575   {
3576     //Get from the absid the supermodule, tower and eta/phi numbers
3577     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
3578     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
3579
3580     //Get the cell energy, if recalibration is on, apply factors
3581     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
3582     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
3583
3584     if(GetCaloUtils()->GetEMCALRecoUtils()->IsRecalibrationOn())
3585     {
3586       recalFactor = GetCaloUtils()->GetEMCALRecoUtils()->GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
3587     }
3588
3589     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
3590
3591     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
3592     // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
3593     if(shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
3594
3595     if(energy > 0 && eCell > eCellMin)
3596     {
3597       w  = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(eCell,energy);
3598
3599       //correct weight, ONLY in simulation
3600       w *= (1 - fWSimu * w );
3601
3602       etai=(Double_t)ieta;
3603       phii=(Double_t)iphi;
3604
3605       if(w > 0.0)
3606       {
3607         wtot += w ;
3608         nstat++;
3609         //Shower shape
3610         sEta     += w * etai * etai ;
3611         etaMean  += w * etai ;
3612         sPhi     += w * phii * phii ;
3613         phiMean  += w * phii ;
3614         sEtaPhi  += w * etai * phii ;
3615       }
3616     }
3617     else if(energy == 0 || (eCellMin <0.01 && eCell == 0)) AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, energy));
3618
3619   }//cell loop
3620
3621   //Normalize to the weight
3622   if (wtot > 0)
3623   {
3624     etaMean /= wtot ;
3625     phiMean /= wtot ;
3626   }
3627   else
3628     AliError(Form("Wrong weight %f\n", wtot));
3629
3630   //Calculate dispersion
3631   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
3632   {
3633     //Get from the absid the supermodule, tower and eta/phi numbers
3634     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
3635     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
3636
3637     //Get the cell energy, if recalibration is on, apply factors
3638     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
3639     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
3640     if (GetCaloUtils()->GetEMCALRecoUtils()->IsRecalibrationOn())
3641     {
3642       recalFactor = GetCaloUtils()->GetEMCALRecoUtils()->GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
3643     }
3644
3645     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
3646
3647     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
3648     // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
3649     if(shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
3650
3651     if(energy > 0 && eCell > eCellMin)
3652     {
3653       w  = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(eCell,energy);
3654
3655       //correct weight, ONLY in simulation
3656       w *= (1 - fWSimu * w );
3657
3658       etai=(Double_t)ieta;
3659       phii=(Double_t)iphi;
3660       if(w > 0.0)
3661       {
3662         disp +=  w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
3663         dEta +=  w * (etai-etaMean)*(etai-etaMean) ;
3664         dPhi +=  w * (phii-phiMean)*(phii-phiMean) ;
3665       }
3666     }
3667     else if(energy == 0 || (eCellMin <0.01 && eCell == 0)) AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, energy));
3668   }// cell loop
3669
3670   //Normalize to the weigth and set shower shape parameters
3671   if (wtot > 0 && nstat > 1)
3672   {
3673     disp    /= wtot ;
3674     dEta    /= wtot ;
3675     dPhi    /= wtot ;
3676     sEta    /= wtot ;
3677     sPhi    /= wtot ;
3678     sEtaPhi /= wtot ;
3679
3680     sEta    -= etaMean * etaMean ;
3681     sPhi    -= phiMean * phiMean ;
3682     sEtaPhi -= etaMean * phiMean ;
3683
3684     l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
3685     l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
3686   }
3687   else
3688   {
3689     l0   = 0. ;
3690     l1   = 0. ;
3691     dEta = 0. ; dPhi = 0. ; disp    = 0. ;
3692     sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
3693   }
3694
3695 }
3696
3697