Changed default weights and scaling between different detectors
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPhoton.cxx
CommitLineData
a3aebfff 1 /**************************************************************************
1c5acb87 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 hereby granted *
c8fe2783 9 * without fee, providGetMixedEvent()ed that the above copyright notice appears in all *
1c5acb87 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/* $Id: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
16
17//_________________________________________________________________________
18//
19// Class for the photon identification.
20// Clusters from calorimeters are identified as photons
21// and kept in the AOD. Few histograms produced.
22//
23// -- Author: Gustavo Conesa (LNF-INFN)
24//////////////////////////////////////////////////////////////////////////////
25
26
27// --- ROOT system ---
28#include <TH2F.h>
477d6cee 29#include <TClonesArray.h>
0c1383b5 30#include <TObjString.h>
477d6cee 31//#include <Riostream.h>
123fc3bd 32#include "TParticle.h"
1c5acb87 33
34// --- Analysis system ---
35#include "AliAnaPhoton.h"
36#include "AliCaloTrackReader.h"
123fc3bd 37#include "AliStack.h"
1c5acb87 38#include "AliCaloPID.h"
6639984f 39#include "AliMCAnalysisUtils.h"
ff45398a 40#include "AliFiducialCut.h"
477d6cee 41#include "AliAODCaloCluster.h"
591cc579 42#include "AliAODMCParticle.h"
c8fe2783 43#include "AliMixedEvent.h"
44
1c5acb87 45
46ClassImp(AliAnaPhoton)
47
48//____________________________________________________________________________
49 AliAnaPhoton::AliAnaPhoton() :
50 AliAnaPartCorrBaseClass(), fCalorimeter(""),
a3aebfff 51 fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
1e86c71e 52 fCheckConversion(kFALSE),fAddConvertedPairsToAOD(kFALSE), fMassCut(0),
2ac125bf 53 fTimeCutMin(-1), fTimeCutMax(9999999), fNCellsCut(0),
a3aebfff 54 fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
7175a03a 55 //MC
123fc3bd 56 fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
591cc579 57 fhPtMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0),
1c5acb87 58 fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
59 fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
6639984f 60 fhPtISR(0),fhPhiISR(0),fhEtaISR(0),
1c5acb87 61 fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
62 fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
63 fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
64 fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
65{
66 //default ctor
67
68 //Initialize parameters
69 InitParameters();
70
71}
78219bac 72/*
1c5acb87 73//____________________________________________________________________________
74AliAnaPhoton::AliAnaPhoton(const AliAnaPhoton & g) :
75 AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter),
1e86c71e 76 fMinDist(g.fMinDist),fMinDist2(g.fMinDist2), fMinDist3(g.fMinDist3),
77 fRejectTrackMatch(g.fRejectTrackMatch),
78 fCheckConversion(g.fCheckConversion),fAddConvertedPairsToAOD(g.fAddConvertedPairsToAOD),
2ac125bf 79 fMassCut(g.fMassCut), fTimeCutMin(g.fTimeCutMin), fTimeCutMax(g.fTimeCutMax), fNCellsCut(g.fNCellsCut),
1e86c71e 80 fhPtPhoton(g.fhPtPhoton),fhPhiPhoton(g.fhPhiPhoton),fhEtaPhoton(g.fhEtaPhoton),
1c5acb87 81 //MC
123fc3bd 82 fhDeltaE(g.fhDeltaE), fhDeltaPt(g.fhDeltaPt),
83 fhRatioE(g.fhRatioE), fhRatioPt(g.fhRatioPt),
84 fh2E(g.fh2E), fh2Pt(g.fh2Pt),
591cc579 85 fhPtMCPhoton(g.fhPtMCPhoton),fhPhiMCPhoton(g.fhPhiMCPhoton),fhEtaMCPhoton(g.fhEtaMCPhoton),
1c5acb87 86 fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
87 fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
6639984f 88 fhPtISR(g.fhPtISR),fhPhiISR(g.fhPhiISR),fhEtaISR(g.fhEtaISR),
1c5acb87 89 fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
90 fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
91 fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
6639984f 92 fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown)
93
1c5acb87 94{
95 // cpy ctor
96
97}
98
99//_________________________________________________________________________
100AliAnaPhoton & AliAnaPhoton::operator = (const AliAnaPhoton & g)
101{
102 // assignment operator
103
104 if(&g == this) return *this;
105
106 fCalorimeter = g.fCalorimeter ;
107 fMinDist = g.fMinDist;
108 fMinDist2 = g.fMinDist2;
109 fMinDist3 = g.fMinDist3;
4cf55759 110 fMassCut = g.fMassCut;
1e86c71e 111
4cf55759 112 fTimeCutMin = g.fTimeCutMin;
113 fTimeCutMax = g.fTimeCutMax;
2ac125bf 114 fNCellsCut = g.fNCellsCut;
115
1e86c71e 116 fRejectTrackMatch = g.fRejectTrackMatch;
117 fCheckConversion = g.fCheckConversion;
118 fAddConvertedPairsToAOD = g.fAddConvertedPairsToAOD;
119
6639984f 120 fhPtPhoton = g.fhPtPhoton ;
1c5acb87 121 fhPhiPhoton = g.fhPhiPhoton ;
122 fhEtaPhoton = g.fhEtaPhoton ;
123fc3bd 123
124 fhDeltaE = g.fhDeltaE;
125 fhDeltaPt = g.fhDeltaPt;
126 fhRatioE = g.fhRatioE;
127 fhRatioPt = g.fhRatioPt;
128 fh2E = g.fh2E;
129 fh2Pt = g.fh2Pt;
1c5acb87 130
591cc579 131 fhPtMCPhoton = g.fhPtMCPhoton ;
132 fhPhiMCPhoton = g.fhPhiMCPhoton ;
133 fhEtaMCPhoton = g.fhEtaMCPhoton ;
6639984f 134 fhPtPrompt = g.fhPtPrompt;
1c5acb87 135 fhPhiPrompt = g.fhPhiPrompt;
136 fhEtaPrompt = g.fhEtaPrompt;
6639984f 137 fhPtFragmentation = g.fhPtFragmentation;
1c5acb87 138 fhPhiFragmentation = g.fhPhiFragmentation;
139 fhEtaFragmentation = g.fhEtaFragmentation;
6639984f 140 fhPtISR = g.fhPtISR;
141 fhPhiISR = g.fhPhiISR;
142 fhEtaISR = g.fhEtaISR;
143 fhPtPi0Decay = g.fhPtPi0Decay;
1c5acb87 144 fhPhiPi0Decay = g.fhPhiPi0Decay;
145 fhEtaPi0Decay = g.fhEtaPi0Decay;
6639984f 146 fhPtOtherDecay = g.fhPtOtherDecay;
1c5acb87 147 fhPhiOtherDecay = g.fhPhiOtherDecay;
148 fhEtaOtherDecay = g.fhEtaOtherDecay;
6639984f 149 fhPtConversion = g. fhPtConversion;
1c5acb87 150 fhPhiConversion = g.fhPhiConversion;
151 fhEtaConversion = g.fhEtaConversion;
6639984f 152 fhPtUnknown = g.fhPtUnknown;
1c5acb87 153 fhPhiUnknown = g.fhPhiUnknown;
154 fhEtaUnknown = g.fhEtaUnknown;
155
156 return *this;
157
158}
78219bac 159*/
1c5acb87 160//____________________________________________________________________________
161AliAnaPhoton::~AliAnaPhoton()
162{
163 //dtor
164
165}
166
0c1383b5 167//________________________________________________________________________
168TObjString * AliAnaPhoton::GetAnalysisCuts()
169{
170 //Save parameters used for analysis
171 TString parList ; //this will be list of parameters used for this analysis.
172 char onePar[255] ;
173
174 sprintf(onePar,"--- AliAnaPhoton ---\n") ;
175 parList+=onePar ;
176 sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
177 parList+=onePar ;
178 sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
179 parList+=onePar ;
180 sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
181 parList+=onePar ;
182 sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
183 parList+=onePar ;
184 sprintf(onePar,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
185 parList+=onePar ;
186
187 //Get parameters set in base class.
188 parList += GetBaseParametersList() ;
189
190 //Get parameters set in PID class.
191 parList += GetCaloPID()->GetPIDParametersList() ;
192
193 //Get parameters set in FiducialCut class (not available yet)
194 //parlist += GetFidCut()->GetFidCutParametersList()
195
196 return new TObjString(parList) ;
197}
198
1c5acb87 199
200//________________________________________________________________________
201TList * AliAnaPhoton::GetCreateOutputObjects()
202{
477d6cee 203 // Create histograms to be saved in output file and
204 // store them in outputContainer
205 TList * outputContainer = new TList() ;
206 outputContainer->SetName("PhotonHistos") ;
4a745797 207
5a2dbc3c 208 Int_t nptbins = GetHistoPtBins();
209 Int_t nphibins = GetHistoPhiBins();
210 Int_t netabins = GetHistoEtaBins();
477d6cee 211 Float_t ptmax = GetHistoPtMax();
212 Float_t phimax = GetHistoPhiMax();
213 Float_t etamax = GetHistoEtaMax();
214 Float_t ptmin = GetHistoPtMin();
215 Float_t phimin = GetHistoPhiMin();
216 Float_t etamin = GetHistoEtaMin();
217
218 //Histograms of highest Photon identified in Event
219 fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
220 fhPtPhoton->SetYTitle("N");
221 fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
222 outputContainer->Add(fhPtPhoton) ;
223
224 fhPhiPhoton = new TH2F
225 ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
226 fhPhiPhoton->SetYTitle("#phi");
227 fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
228 outputContainer->Add(fhPhiPhoton) ;
229
230 fhEtaPhoton = new TH2F
231 ("hEtaPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
232 fhEtaPhoton->SetYTitle("#eta");
233 fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
234 outputContainer->Add(fhEtaPhoton) ;
235
236 if(IsDataMC()){
123fc3bd 237 fhDeltaE = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
238 fhDeltaE->SetXTitle("#Delta E (GeV)");
239 outputContainer->Add(fhDeltaE);
240
241 fhDeltaPt = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
242 fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
243 outputContainer->Add(fhDeltaPt);
244
245 fhRatioE = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
246 fhRatioE->SetXTitle("E_{reco}/E_{gen}");
247 outputContainer->Add(fhRatioE);
477d6cee 248
123fc3bd 249 fhRatioPt = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
250 fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
251 outputContainer->Add(fhRatioPt);
252
253 fh2E = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 254 fh2E->SetXTitle("E_{rec} (GeV)");
255 fh2E->SetYTitle("E_{gen} (GeV)");
123fc3bd 256 outputContainer->Add(fh2E);
257
258 fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
41e886c8 259 fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
260 fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
123fc3bd 261 outputContainer->Add(fh2Pt);
262
c8fe2783 263 fhPtMCPhoton = new TH1F("hPtMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
264 fhPtMCPhoton->SetYTitle("N");
265 fhPtMCPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
266 outputContainer->Add(fhPtMCPhoton) ;
267
268 fhPhiMCPhoton = new TH2F
269 ("hPhiMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
270 fhPhiMCPhoton->SetYTitle("#phi");
271 fhPhiMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
272 outputContainer->Add(fhPhiMCPhoton) ;
273
274 fhEtaMCPhoton = new TH2F
275 ("hEtaMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
276 fhEtaMCPhoton->SetYTitle("#eta");
277 fhEtaMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
278 outputContainer->Add(fhEtaMCPhoton) ;
279
591cc579 280 fhPtPrompt = new TH1F("hPtMCPrompt","Number of prompt #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 281 fhPtPrompt->SetYTitle("N");
282 fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
283 outputContainer->Add(fhPtPrompt) ;
284
285 fhPhiPrompt = new TH2F
a3aebfff 286 ("hPhiMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 287 fhPhiPrompt->SetYTitle("#phi");
288 fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
289 outputContainer->Add(fhPhiPrompt) ;
290
291 fhEtaPrompt = new TH2F
a3aebfff 292 ("hEtaMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 293 fhEtaPrompt->SetYTitle("#eta");
294 fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
295 outputContainer->Add(fhEtaPrompt) ;
296
591cc579 297 fhPtFragmentation = new TH1F("hPtMCFragmentation","Number of fragmentation #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 298 fhPtFragmentation->SetYTitle("N");
299 fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
300 outputContainer->Add(fhPtFragmentation) ;
301
302 fhPhiFragmentation = new TH2F
a3aebfff 303 ("hPhiMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 304 fhPhiFragmentation->SetYTitle("#phi");
305 fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
306 outputContainer->Add(fhPhiFragmentation) ;
307
308 fhEtaFragmentation = new TH2F
a3aebfff 309 ("hEtaMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 310 fhEtaFragmentation->SetYTitle("#eta");
311 fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
312 outputContainer->Add(fhEtaFragmentation) ;
313
a3aebfff 314 fhPtISR = new TH1F("hPtMCISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 315 fhPtISR->SetYTitle("N");
316 fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
317 outputContainer->Add(fhPtISR) ;
318
319 fhPhiISR = new TH2F
a3aebfff 320 ("hPhiMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 321 fhPhiISR->SetYTitle("#phi");
322 fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
323 outputContainer->Add(fhPhiISR) ;
324
325 fhEtaISR = new TH2F
a3aebfff 326 ("hEtaMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 327 fhEtaISR->SetYTitle("#eta");
328 fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
329 outputContainer->Add(fhEtaISR) ;
330
591cc579 331 fhPtPi0Decay = new TH1F("hPtMCPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 332 fhPtPi0Decay->SetYTitle("N");
333 fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
334 outputContainer->Add(fhPtPi0Decay) ;
335
336 fhPhiPi0Decay = new TH2F
a3aebfff 337 ("hPhiMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 338 fhPhiPi0Decay->SetYTitle("#phi");
339 fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
340 outputContainer->Add(fhPhiPi0Decay) ;
341
342 fhEtaPi0Decay = new TH2F
a3aebfff 343 ("hEtaMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 344 fhEtaPi0Decay->SetYTitle("#eta");
345 fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
346 outputContainer->Add(fhEtaPi0Decay) ;
347
a3aebfff 348 fhPtOtherDecay = new TH1F("hPtMCOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 349 fhPtOtherDecay->SetYTitle("N");
350 fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
351 outputContainer->Add(fhPtOtherDecay) ;
352
353 fhPhiOtherDecay = new TH2F
a3aebfff 354 ("hPhiMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 355 fhPhiOtherDecay->SetYTitle("#phi");
356 fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
357 outputContainer->Add(fhPhiOtherDecay) ;
358
359 fhEtaOtherDecay = new TH2F
a3aebfff 360 ("hEtaMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 361 fhEtaOtherDecay->SetYTitle("#eta");
362 fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
363 outputContainer->Add(fhEtaOtherDecay) ;
364
a3aebfff 365 fhPtConversion = new TH1F("hPtMCConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 366 fhPtConversion->SetYTitle("N");
367 fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
368 outputContainer->Add(fhPtConversion) ;
369
370 fhPhiConversion = new TH2F
a3aebfff 371 ("hPhiMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 372 fhPhiConversion->SetYTitle("#phi");
373 fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
374 outputContainer->Add(fhPhiConversion) ;
375
376 fhEtaConversion = new TH2F
a3aebfff 377 ("hEtaMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 378 fhEtaConversion->SetYTitle("#eta");
379 fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
380 outputContainer->Add(fhEtaConversion) ;
381
a3aebfff 382 fhPtUnknown = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
477d6cee 383 fhPtUnknown->SetYTitle("N");
384 fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
385 outputContainer->Add(fhPtUnknown) ;
386
387 fhPhiUnknown = new TH2F
a3aebfff 388 ("hPhiMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
477d6cee 389 fhPhiUnknown->SetYTitle("#phi");
390 fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
391 outputContainer->Add(fhPhiUnknown) ;
392
393 fhEtaUnknown = new TH2F
a3aebfff 394 ("hEtaMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
477d6cee 395 fhEtaUnknown->SetYTitle("#eta");
396 fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
397 outputContainer->Add(fhEtaUnknown) ;
a3aebfff 398
477d6cee 399 }//Histos with MC
0c1383b5 400
477d6cee 401 return outputContainer ;
402
1c5acb87 403}
404
405//____________________________________________________________________________
6639984f 406void AliAnaPhoton::Init()
407{
408
409 //Init
410 //Do some checks
1e86c71e 411 if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
591cc579 412 printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 413 abort();
414 }
1e86c71e 415 else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
591cc579 416 printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
6639984f 417 abort();
418 }
419
420}
421
422
423//____________________________________________________________________________
1c5acb87 424void AliAnaPhoton::InitParameters()
425{
426
427 //Initialize the parameters of the analysis.
a3aebfff 428 AddToHistogramsName("AnaPhoton_");
429
1c5acb87 430 fCalorimeter = "PHOS" ;
431 fMinDist = 2.;
432 fMinDist2 = 4.;
433 fMinDist3 = 5.;
1e86c71e 434 fMassCut = 0.03; //30 MeV
435
4cf55759 436 fTimeCutMin = -1;
437 fTimeCutMax = 9999999;
2ac125bf 438 fNCellsCut = 0;
439
1e86c71e 440 fRejectTrackMatch = kTRUE ;
441 fCheckConversion = kFALSE;
442 fAddConvertedPairsToAOD = kFALSE;
443
1c5acb87 444}
445
446//__________________________________________________________________
447void AliAnaPhoton::MakeAnalysisFillAOD()
448{
c8fe2783 449 //Do analysis and fill aods
450 //Search for photons in fCalorimeter
451
452 //Get vertex for photon momentum calculation
477d6cee 453
c8fe2783 454 for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) {
455 if (!GetMixedEvent())
456 GetReader()->GetVertex(GetVertex(iev));
457 else
458 GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev));
459 }
460
1e68a3f4 461// Double_t vertex2[] = {0,0,0} ; //vertex from second input aod
462//
463// if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
464// {
465// if(GetReader()->GetSecondInputAODTree())
466// GetReader()->GetSecondInputAODVertex(vertex2);
467// }
c8fe2783 468 //printf("Vertex 0: %f,%f,%f\n",vertex[0],vertex[1],vertex[2]);
469 //printf("Vertex 1: %f,%f,%f\n",vertex2[0],vertex2[1],vertex2[2]);
470 //Select the Calorimeter of the photon
471 TObjArray * pl = 0x0;
477d6cee 472 if(fCalorimeter == "PHOS")
473 pl = GetAODPHOS();
474 else if (fCalorimeter == "EMCAL")
475 pl = GetAODEMCAL();
476
c8fe2783 477 //Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods
1e86c71e 478 TLorentzVector mom, mom2 ;
479 Int_t nCaloClusters = pl->GetEntriesFast();
480 Bool_t * indexConverted = new Bool_t[nCaloClusters];
c8fe2783 481 for (Int_t i = 0; i < nCaloClusters; i++)
482 indexConverted[i] = kFALSE;
1e86c71e 483
484 for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
485
486 AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));
c8fe2783 487
488 Int_t evtIndex = 0 ;
489 if (GetMixedEvent()) {
490 evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
491 }
492 //Cluster selection, not charged, with photon id and in fiducial cut
233e0df8 493
c8fe2783 494 //Input from second AOD?
495 Int_t input = 0;
1e68a3f4 496// if (fCalorimeter == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= icalo)
497// input = 1 ;
498// else if(fCalorimeter == "PHOS" && GetReader()->GetAODPHOSNormalInputEntries() <= icalo)
499// input = 1;
233e0df8 500
c8fe2783 501 //Get Momentum vector,
502 if (input == 0)
503 calo->GetMomentum(mom,GetVertex(evtIndex)) ;//Assume that come from vertex in straight line
1e68a3f4 504// else if(input == 1)
505// calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
c8fe2783 506
507 //If too small or big pt, skip it
477d6cee 508 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
4cf55759 509 Double_t tof = calo->GetTOF()*1e9;
c8fe2783 510
511 if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
2ac125bf 512
c8fe2783 513 if(calo->GetNCells() <= fNCellsCut) continue;
2ac125bf 514
c8fe2783 515 //printf("AliAnaPhoton::Current Event %d; Current File Name : %s, E %f, pT %f, Ecl %f\n",GetReader()->GetEventNumber(),(GetReader()->GetCurrentFileName()).Data(), mom.E(), mom.Pt(),calo->E());
516
517 //Check acceptance selection
ff45398a 518 if(IsFiducialCutOn()){
519 Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
477d6cee 520 if(! in ) continue ;
521 }
c8fe2783 522
523 //Create AOD for analysis
477d6cee 524 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
c8fe2783 525 Int_t label = calo->GetLabel();
591cc579 526 aodph.SetLabel(label);
c8fe2783 527 aodph.SetInputFileIndex(input);
528
529 //printf("Index %d, Id %d\n",icalo, calo->GetID());
530 //Set the indeces of the original caloclusters
477d6cee 531 aodph.SetCaloLabel(calo->GetID(),-1);
532 aodph.SetDetector(fCalorimeter);
533 if(GetDebug() > 1)
ff45398a 534 printf("AliAnaPhoton::MakeAnalysisFillAOD() - Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());
477d6cee 535
c8fe2783 536 //Check Distance to Bad channel, set bit.
537 Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
477d6cee 538 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
539 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
540 continue ;
541
a3aebfff 542 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
477d6cee 543
544 if(distBad > fMinDist3) aodph.SetDistToBad(2) ;
545 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
546 else aodph.SetDistToBad(0) ;
547
c8fe2783 548 //Skip matched clusters with tracks
549 if(fRejectTrackMatch && calo->GetNTracksMatched() > 0) continue ;
550
551 //Check PID
552 //PID selection or bit setting
477d6cee 553 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
c8fe2783 554 //Get most probable PID, check PID weights (in MC this option is mandatory)
555 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
a3aebfff 556 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
c8fe2783 557 //If primary is not photon, skip it.
477d6cee 558 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
559 }
560 else if(IsCaloPIDOn()){
c8fe2783 561
562 //Get most probable PID, 2 options check PID weights
563 //or redo PID, recommended option for EMCal.
477d6cee 564 if(!IsCaloPIDRecalculationOn())
c8fe2783 565 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
477d6cee 566 else
c8fe2783 567 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
477d6cee 568
a3aebfff 569 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
477d6cee 570
c8fe2783 571 //If cluster does not pass pid, not photon, skip it.
477d6cee 572 if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
573
574 }
575 else{
c8fe2783 576 //Set PID bits for later selection (AliAnaPi0 for example)
577 //GetPDG already called in SetPIDBits.
477d6cee 578 GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);
a3aebfff 579 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
477d6cee 580 }
581
a3aebfff 582 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
477d6cee 583
c8fe2783 584 //Play with the MC stack if available
585 //Check origin of the candidates
477d6cee 586 if(IsDataMC()){
c8fe2783 587
588 aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
591cc579 589 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
477d6cee 590 }//Work with stack also
591
c8fe2783 592
593 // Check if cluster comes from a conversion in the material in front of the calorimeter
594 // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
595
596 if(fCheckConversion && nCaloClusters > 1){
597 Bool_t bConverted = kFALSE;
598 Int_t id2 = -1;
1e86c71e 599
c8fe2783 600 //Check if set previously as converted couple, if so skip its use.
601 if (indexConverted[icalo]) continue;
1e86c71e 602
c8fe2783 603 for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
604 //Check if set previously as converted couple, if so skip its use.
605 if (indexConverted[jcalo]) continue;
606 //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
607 AliAODCaloCluster * calo2 = (AliAODCaloCluster*) (pl->At(jcalo)); //Get cluster kinematics
608 Int_t evtIndex2 = 0 ;
609 if (GetMixedEvent()) {
610 evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
611 }
612 calo2->GetMomentum(mom2,GetVertex(evtIndex2));
613 //Check only certain regions
614 Bool_t in2 = kTRUE;
615 if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
616 if(!in2) continue;
617
618 //Get mass of pair, if small, take this pair.
619 //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut);
620 if((mom+mom2).M() < fMassCut){
621 bConverted = kTRUE;
622 id2 = calo2->GetID();
623 indexConverted[jcalo]=kTRUE;
624 break;
625 }
1e86c71e 626
c8fe2783 627 }//Mass loop
1e86c71e 628
c8fe2783 629 if(bConverted){
630 if(fAddConvertedPairsToAOD){
631 //Create AOD of pair analysis
632 TLorentzVector mpair = mom+mom2;
633 AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
634 aodpair.SetLabel(aodph.GetLabel());
635 aodpair.SetInputFileIndex(input);
636
637 //printf("Index %d, Id %d\n",icalo, calo->GetID());
638 //Set the indeces of the original caloclusters
639 aodpair.SetCaloLabel(calo->GetID(),id2);
640 aodpair.SetDetector(fCalorimeter);
641 aodpair.SetPdg(aodph.GetPdg());
642 aodpair.SetTag(aodph.GetTag());
643
644 //Add AOD with pair object to aod branch
645 AddAODParticle(aodpair);
646 //printf("\t \t both added pair\n");
647 }
648
649 //Do not add the current calocluster
650 continue;
651 }//converted pair
652 }//check conversion
653 //printf("\t \t added single cluster %d\n",icalo);
1e86c71e 654
c8fe2783 655 //Add AOD with photon object to aod branch
477d6cee 656 AddAODParticle(aodph);
657
658 }//loop
659
4a745797 660 delete [] indexConverted;
661
a3aebfff 662 if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() End fill AODs \n");
477d6cee 663
1c5acb87 664}
665
666//__________________________________________________________________
667void AliAnaPhoton::MakeAnalysisFillHistograms()
668{
123fc3bd 669 //Do analysis and fill histograms
670
591cc579 671 // Access MC information in stack if requested, check that it exists.
672 AliStack * stack = 0x0;
673 TParticle * primary = 0x0;
674 TClonesArray * mcparticles0 = 0x0;
675 TClonesArray * mcparticles1 = 0x0;
676 AliAODMCParticle * aodprimary = 0x0;
677 if(IsDataMC()){
678
679 if(GetReader()->ReadStack()){
680 stack = GetMCStack() ;
681 if(!stack) {
682 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
683 abort();
684 }
685
686 }
687 else if(GetReader()->ReadAODMCParticles()){
688
689 //Get the list of MC particles
690 mcparticles0 = GetReader()->GetAODMCParticles(0);
691 if(!mcparticles0 && GetDebug() > 0) {
692 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
693 }
1e68a3f4 694// if(GetReader()->GetSecondInputAODTree()){
695// mcparticles1 = GetReader()->GetAODMCParticles(1);
696// if(!mcparticles1 && GetDebug() > 0) {
697// printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
698// }
699// }
591cc579 700
701 }
702 }// is data and MC
703
123fc3bd 704 //Loop on stored AOD photons
705 Int_t naod = GetOutputAODBranch()->GetEntriesFast();
706 if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
707
708 for(Int_t iaod = 0; iaod < naod ; iaod++){
c8fe2783 709 AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
710 Int_t pdg = ph->GetPdg();
711
712 if(GetDebug() > 3)
713 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
714
715 //If PID used, fill histos with photons in Calorimeter fCalorimeter
716 if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
717 if(ph->GetDetector() != fCalorimeter) continue;
718
719 if(GetDebug() > 2)
720 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
721
722 //Fill photon histograms
723 Float_t ptcluster = ph->Pt();
724 Float_t phicluster = ph->Phi();
725 Float_t etacluster = ph->Eta();
726 Float_t ecluster = ph->E();
727
728 fhPtPhoton ->Fill(ptcluster);
729 fhPhiPhoton ->Fill(ptcluster,phicluster);
730 fhEtaPhoton ->Fill(ptcluster,etacluster);
731
732 //Play with the MC data if available
733 if(IsDataMC()){
734
735 Int_t tag =ph->GetTag();
736
737 if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
738 {
739 fhPtMCPhoton ->Fill(ptcluster);
740 fhPhiMCPhoton ->Fill(ptcluster,phicluster);
741 fhEtaMCPhoton ->Fill(ptcluster,etacluster);
123fc3bd 742
c8fe2783 743 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
744 {
745 fhPtConversion ->Fill(ptcluster);
746 fhPhiConversion ->Fill(ptcluster,phicluster);
747 fhEtaConversion ->Fill(ptcluster,etacluster);
748 }
591cc579 749
c8fe2783 750 if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
751 fhPtPrompt ->Fill(ptcluster);
752 fhPhiPrompt ->Fill(ptcluster,phicluster);
753 fhEtaPrompt ->Fill(ptcluster,etacluster);
591cc579 754 }
c8fe2783 755 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
756 {
757 fhPtFragmentation ->Fill(ptcluster);
758 fhPhiFragmentation ->Fill(ptcluster,phicluster);
759 fhEtaFragmentation ->Fill(ptcluster,etacluster);
760 }
761 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
762 {
763 fhPtISR ->Fill(ptcluster);
764 fhPhiISR ->Fill(ptcluster,phicluster);
765 fhEtaISR ->Fill(ptcluster,etacluster);
766 }
767 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
768 {
769 fhPtPi0Decay ->Fill(ptcluster);
770 fhPhiPi0Decay ->Fill(ptcluster,phicluster);
771 fhEtaPi0Decay ->Fill(ptcluster,etacluster);
772 }
773 else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
774 {
775 fhPtOtherDecay ->Fill(ptcluster);
776 fhPhiOtherDecay ->Fill(ptcluster,phicluster);
777 fhEtaOtherDecay ->Fill(ptcluster,etacluster);
778 }
779 }
780 else{
781 fhPtUnknown ->Fill(ptcluster);
782 fhPhiUnknown ->Fill(ptcluster,phicluster);
783 fhEtaUnknown ->Fill(ptcluster,etacluster);
784
591cc579 785// printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
786// ph->GetLabel(),ph->Pt());
787// for(Int_t i = 0; i < 20; i++) {
788// if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
789// }
790// printf("\n");
791
c8fe2783 792 }
793
794
795 // Access MC information in stack if requested, check that it exists.
796 Int_t label =ph->GetLabel();
797 if(label < 0) {
798 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);
799 continue;
800 }
801
802 Float_t eprim = 0;
803 Float_t ptprim = 0;
804 if(GetReader()->ReadStack()){
805
806 if(label >= stack->GetNtrack()) {
807 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack());
808 continue ;
809 }
810
811 primary = stack->Particle(label);
812 if(!primary){
813 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
123fc3bd 814 continue;
c8fe2783 815 }
816 eprim = primary->Energy();
817 ptprim = primary->Pt();
818
819 }
820 else if(GetReader()->ReadAODMCParticles()){
821 //Check which is the input
822 if(ph->GetInputFileIndex() == 0){
823 if(!mcparticles0) continue;
824 if(label >= mcparticles0->GetEntriesFast()) {
825 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
826 label, mcparticles0->GetEntriesFast());
827 continue ;
591cc579 828 }
c8fe2783 829 //Get the particle
830 aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
123fc3bd 831
c8fe2783 832 }
833 else {//Second input
834 if(!mcparticles1) continue;
835 if(label >= mcparticles1->GetEntriesFast()) {
836 if(GetDebug() > 2) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",
837 label, mcparticles1->GetEntriesFast());
838 continue ;
591cc579 839 }
c8fe2783 840 //Get the particle
841 aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
591cc579 842
c8fe2783 843 }//second input
844
845 if(!aodprimary){
846 printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label);
847 continue;
848 }
849
850 eprim = aodprimary->E();
851 ptprim = aodprimary->Pt();
852
853 }
854
855 fh2E ->Fill(ecluster, eprim);
856 fh2Pt ->Fill(ptcluster, ptprim);
857 fhDeltaE ->Fill(eprim-ecluster);
858 fhDeltaPt->Fill(ptprim-ptcluster);
859 if(eprim > 0) fhRatioE ->Fill(ecluster/eprim);
860 if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
861
862 }//Histograms with MC
863
864 }// aod loop
591cc579 865
1c5acb87 866}
867
868
869//__________________________________________________________________
870void AliAnaPhoton::Print(const Option_t * opt) const
871{
477d6cee 872 //Print some relevant parameters set for the analysis
873
874 if(! opt)
875 return;
876
877 printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
878 AliAnaPartCorrBaseClass::Print(" ");
a3aebfff 879
477d6cee 880 printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
881 printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
882 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
883 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
a3aebfff 884 printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1e86c71e 885 printf("Check Pair Conversion = %d\n",fCheckConversion);
886 printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD);
887 printf("Conversion pair mass cut = %f\n",fMassCut);
4cf55759 888 printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax);
2ac125bf 889 printf("Number of cells in cluster is > %d \n", fNCellsCut);
477d6cee 890 printf(" \n") ;
1c5acb87 891
892}