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