PWGGA/CaloTrackCorrelations/AliAnaInsideClusterInvariantMass.cxx

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