AliAnaPi0: Add possibility to cut on opening angle of pairs, to reject background...
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPi0.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 /* $Id: $ */
16
17 //_________________________________________________________________________
18 // Class to collect two-photon invariant mass distributions for
19 // extractin raw pi0 yield.
20 //
21 //-- Author: Dmitri Peressounko (RRC "KI") 
22 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
23 //-- and Gustavo Conesa (INFN-Frascati)
24 //_________________________________________________________________________
25
26
27 // --- ROOT system ---
28 #include "TH3.h"
29 #include "TH2D.h"
30 //#include "Riostream.h"
31 #include "TCanvas.h"
32 #include "TPad.h"
33 #include "TROOT.h"
34 #include "TClonesArray.h"
35 #include "TObjString.h"
36
37 //---- AliRoot system ----
38 #include "AliAnaPi0.h"
39 #include "AliCaloTrackReader.h"
40 #include "AliCaloPID.h"
41 #include "AliStack.h"
42 #include "AliFiducialCut.h"
43 #include "TParticle.h"
44 #include "AliAODCaloCluster.h"
45 #include "AliVEvent.h"
46 #include "AliESDCaloCluster.h"
47 #include "AliESDEvent.h"
48 #include "AliAODEvent.h"
49 #include "AliNeutralMesonSelection.h"
50
51 ClassImp(AliAnaPi0)
52
53 //________________________________________________________________________________________________________________________________________________  
54 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
55 fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
56 fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
57 fNModules(12), fUseAngleCut(kFALSE), fEventsList(0x0), //fhEtalon(0x0), 
58 fhReMod(0x0), fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0),
59 fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0),
60 fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0),
61 fhPrimOpeningAngle(0x0),fhPrimCosOpeningAngle(0x0)
62 {
63 //Default Ctor
64  InitParameters();
65  
66 }
67
68 //________________________________________________________________________________________________________________________________________________
69 AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnaPartCorrBaseClass(ex),  
70 fNCentrBin(ex.fNCentrBin),fNZvertBin(ex.fNZvertBin),fNrpBin(ex.fNrpBin),
71 fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv),fZvtxCut(ex.fZvtxCut), fCalorimeter(ex.fCalorimeter),
72 fNModules(ex.fNModules), fUseAngleCut(ex.fUseAngleCut), fEventsList(ex.fEventsList), //fhEtalon(ex.fhEtalon), 
73 fhReMod(ex.fhReMod), fhRe1(ex.fhRe1),fhMi1(ex.fhMi1),fhRe2(ex.fhRe2),fhMi2(ex.fhMi2),
74 fhRe3(ex.fhRe3),fhMi3(ex.fhMi3),fhEvents(ex.fhEvents),
75 fhRealOpeningAngle(ex.fhRealOpeningAngle),fhRealCosOpeningAngle(ex.fhRealCosOpeningAngle),
76 fhPrimPt(ex.fhPrimPt), fhPrimAccPt(ex.fhPrimAccPt), fhPrimY(ex.fhPrimY), 
77 fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi),
78 fhPrimOpeningAngle(ex.fhPrimOpeningAngle),fhPrimCosOpeningAngle(ex.fhPrimCosOpeningAngle)
79 {
80   // cpy ctor
81   //Do not need it
82 }
83
84 //________________________________________________________________________________________________________________________________________________
85 AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & ex)
86 {
87   // assignment operator
88   
89   if(this == &ex)return *this;
90   ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
91   
92   fNCentrBin = ex.fNCentrBin  ; fNZvertBin = ex.fNZvertBin  ; fNrpBin = ex.fNrpBin  ; 
93   fNPID = ex.fNPID  ; fNmaxMixEv = ex.fNmaxMixEv  ; fZvtxCut = ex.fZvtxCut  ;  fCalorimeter = ex.fCalorimeter  ;  
94   fNModules = ex.fNModules; fEventsList = ex.fEventsList  ;  //fhEtalon = ex.fhEtalon  ; 
95   fhRe1 = ex.fhRe1  ; fhMi1 = ex.fhMi1  ; fhRe2 = ex.fhRe2  ; fhMi2 = ex.fhMi2  ; fhReMod = ex.fhReMod; 
96   fhRe3 = ex.fhRe3  ; fhMi3 = ex.fhMi3  ; fhEvents = ex.fhEvents  ; fUseAngleCut = ex.fUseAngleCut;
97   fhPrimPt = ex.fhPrimPt  ;  fhPrimAccPt = ex.fhPrimAccPt  ;  fhPrimY = ex.fhPrimY  ;  
98   fhPrimAccY = ex.fhPrimAccY  ;  fhPrimPhi = ex.fhPrimPhi  ;  fhPrimAccPhi = ex.fhPrimAccPhi ;
99   fhRealOpeningAngle = ex.fhRealOpeningAngle; fhRealCosOpeningAngle = ex.fhRealCosOpeningAngle;
100   fhPrimOpeningAngle = ex.fhPrimOpeningAngle; fhPrimCosOpeningAngle = ex.fhPrimCosOpeningAngle;
101
102   return *this;
103   
104 }
105
106 //________________________________________________________________________________________________________________________________________________
107 AliAnaPi0::~AliAnaPi0() {
108   // Remove event containers
109   if(fEventsList){
110     for(Int_t ic=0; ic<fNCentrBin; ic++){
111       for(Int_t iz=0; iz<fNZvertBin; iz++){
112         for(Int_t irp=0; irp<fNrpBin; irp++){
113           fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp]->Delete() ;
114           delete fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] ;
115         }
116       }
117     }
118     delete[] fEventsList; 
119     fEventsList=0 ;
120   }
121         
122 }
123
124 //________________________________________________________________________________________________________________________________________________
125 void AliAnaPi0::InitParameters()
126 {
127 //Init parameters when first called the analysis
128 //Set default parameters
129   SetInputAODName("PWG4Particle");
130   
131   AddToHistogramsName("AnaPi0_");
132   fNModules = 12; // set maximum to maximum number of EMCAL modules
133   fNCentrBin = 1;
134   fNZvertBin = 1;
135   fNrpBin    = 1;
136   fNPID      = 9;
137   fNmaxMixEv = 10;
138   fZvtxCut   = 40;
139   fCalorimeter  = "PHOS";
140   fUseAngleCut = kFALSE;
141         
142 }
143 //________________________________________________________________________________________________________________________________________________
144 //void AliAnaPi0::Init()
145 //{  
146   //Init some data members needed in analysis
147   
148   //Histograms binning and range
149  // if(!fhEtalon){                                                   //  p_T      alpha   d m_gg    
150  //   fhEtalon = new TH3D("hEtalon","Histo with binning parameters",50,0.,25.,10,0.,1.,200,0.,1.) ; 
151  //   fhEtalon->SetXTitle("P_{T} (GeV)") ;
152  //   fhEtalon->SetYTitle("#alpha") ;
153  //   fhEtalon->SetZTitle("m_{#gamma#gamma} (GeV)") ;
154  // }
155   
156 //}
157
158 //________________________________________________________________________________________________________________________________________________
159 TList * AliAnaPi0::GetCreateOutputObjects()
160 {  
161   // Create histograms to be saved in output file and 
162   // store them in fOutputContainer
163   
164   //create event containers
165   fEventsList = new TList*[fNCentrBin*fNZvertBin*fNrpBin] ;
166         
167   for(Int_t ic=0; ic<fNCentrBin; ic++){
168     for(Int_t iz=0; iz<fNZvertBin; iz++){
169       for(Int_t irp=0; irp<fNrpBin; irp++){
170         fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] = new TList() ;
171       }
172     }
173   }
174   
175   //If Geometry library loaded, do geometry selection during analysis.
176   if(fCalorimeter=="PHOS"){
177     if(!GetReader()->GetPHOSGeometry()) printf("AliAnaPi0::GetCreateOutputObjects() - Initialize PHOS geometry!\n");
178     GetReader()->InitPHOSGeometry();
179     
180   }
181   else if(fCalorimeter=="EMCAL"){
182     if(!GetReader()->GetEMCALGeometry()) printf("AliAnaPi0::GetCreateOutputObjects() - Initialize EMCAL geometry!\n");
183     GetReader()->InitEMCALGeometry();
184   }
185         
186   TList * outputContainer = new TList() ; 
187   outputContainer->SetName(GetName()); 
188         
189   fhReMod = new TH3D*[fNModules] ;
190   fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
191   fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
192   fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
193   fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
194   fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
195   fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
196   
197   char key[255] ;
198   char title[255] ;
199   Int_t nptbins   = GetHistoPtBins();
200   Int_t nphibins  = GetHistoPhiBins();
201   Int_t netabins  = GetHistoEtaBins();
202   Float_t ptmax   = GetHistoPtMax();
203   Float_t phimax  = GetHistoPhiMax();
204   Float_t etamax  = GetHistoEtaMax();
205   Float_t ptmin   = GetHistoPtMin();
206   Float_t phimin  = GetHistoPhiMin();
207   Float_t etamin  = GetHistoEtaMin();   
208         
209   Int_t nmassbins = GetHistoMassBins();
210   Int_t nasymbins = GetHistoAsymmetryBins();
211   Float_t massmax = GetHistoMassMax();
212   Float_t asymmax = GetHistoAsymmetryMax();
213   Float_t massmin = GetHistoMassMin();
214   Float_t asymmin = GetHistoAsymmetryMin();
215         
216   for(Int_t ic=0; ic<fNCentrBin; ic++){
217     for(Int_t ipid=0; ipid<fNPID; ipid++){
218                 
219       //Distance to bad module 1
220       sprintf(key,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
221       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
222       
223       //fhEtalon->Clone(key);
224       //fhRe1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
225       //fhRe1[ic*fNPID+ipid]->SetName(key) ;
226       //fhRe1[ic*fNPID+ipid]->SetTitle(title) ;
227           fhRe1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
228       outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
229       
230       sprintf(key,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
231       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
232       //fhMi1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
233       //fhMi1[ic*fNPID+ipid]->SetName(key) ;
234       //fhMi1[ic*fNPID+ipid]->SetTitle(title) ;
235           fhMi1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
236       outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
237       
238       //Distance to bad module 2
239       sprintf(key,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
240       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
241       //fhRe2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
242       //fhRe2[ic*fNPID+ipid]->SetName(key) ;
243       //fhRe2[ic*fNPID+ipid]->SetTitle(title) ;
244           fhRe2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
245       outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
246       
247       sprintf(key,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
248       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
249       //fhMi2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
250       //fhMi2[ic*fNPID+ipid]->SetName(key) ;
251       //fhMi2[ic*fNPID+ipid]->SetTitle(title) ;
252           fhMi2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
253       outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
254       
255       //Distance to bad module 3
256       sprintf(key,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
257       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
258       //fhRe3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
259       //fhRe3[ic*fNPID+ipid]->SetName(key) ; 
260       //fhRe3[ic*fNPID+ipid]->SetTitle(title) ;
261           fhRe3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
262       outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
263       
264       sprintf(key,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
265       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
266       //fhMi3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
267       //fhMi3[ic*fNPID+ipid]->SetName(key) ;
268       //fhMi3[ic*fNPID+ipid]->SetTitle(title) ;
269           fhMi3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
270       outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
271     }
272   }
273   
274   
275   fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
276                     fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
277   outputContainer->Add(fhEvents) ;
278   
279         
280   fhRealOpeningAngle  = new TH2D
281   ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5); 
282   fhRealOpeningAngle->SetYTitle("#theta(rad)");
283   fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
284   outputContainer->Add(fhRealOpeningAngle) ;
285
286   fhRealCosOpeningAngle  = new TH2D
287   ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1); 
288   fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
289   fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
290   outputContainer->Add(fhRealCosOpeningAngle) ;
291         
292         
293   //Histograms filled only if MC data is requested      
294   if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
295    // if(fhEtalon->GetXaxis()->GetXbins() && fhEtalon->GetXaxis()->GetXbins()->GetSize()){ //Variable bin size
296    //   fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
297    //   fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",fhEtalon->GetXaxis()->GetNbins(),
298         //                   fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
299    // }
300    // else{
301    //   fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
302    //   fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",
303         //                   fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
304    // }
305
306         fhPrimPt     = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
307     fhPrimAccPt  = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
308     outputContainer->Add(fhPrimPt) ;
309     outputContainer->Add(fhPrimAccPt) ;
310     
311     fhPrimY      = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ; 
312     outputContainer->Add(fhPrimY) ;
313     
314     fhPrimAccY   = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ; 
315     outputContainer->Add(fhPrimAccY) ;
316     
317         fhPrimPhi    = new TH1D("hPrimaryPhi","Azimithal of primary pi0",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; 
318     outputContainer->Add(fhPrimPhi) ;
319     
320     fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; 
321     outputContainer->Add(fhPrimAccPhi) ;
322     
323     
324     fhPrimOpeningAngle  = new TH2D
325       ("hPrimOpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5); 
326     fhPrimOpeningAngle->SetYTitle("#theta(rad)");
327     fhPrimOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
328     outputContainer->Add(fhPrimOpeningAngle) ;
329     
330     fhPrimCosOpeningAngle  = new TH2D
331       ("hPrimCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1); 
332     fhPrimCosOpeningAngle->SetYTitle("cos (#theta) ");
333     fhPrimCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
334     outputContainer->Add(fhPrimCosOpeningAngle) ;
335     
336   }
337   
338   for(Int_t imod=0; imod<fNModules; imod++){
339     //Module dependent invariant mass
340     sprintf(key,"hReMod_%d",imod) ;
341     sprintf(title,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
342     //fhEtalon->Clone(key);
343     //fhReMod[imod]=(TH3D*)fhEtalon->Clone(key) ;
344     //fhReMod[imod]->SetName(key) ;
345     //fhReMod[imod]->SetTitle(title) ;
346     fhReMod[imod]  = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
347     outputContainer->Add(fhReMod[imod]) ;
348   }
349   
350   
351   //Save parameters used for analysis
352   TString parList ; //this will be list of parameters used for this analysis.
353   char onePar[255] ;
354   sprintf(onePar,"--- AliAnaPi0 ---\n") ;
355   parList+=onePar ;     
356   sprintf(onePar,"Number of bins in Centrality:  %d \n",fNCentrBin) ;
357   parList+=onePar ;
358   sprintf(onePar,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
359   parList+=onePar ;
360   sprintf(onePar,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
361   parList+=onePar ;
362   sprintf(onePar,"Depth of event buffer: %d \n",fNmaxMixEv) ;
363   parList+=onePar ;
364   sprintf(onePar,"Number of different PID used:  %d \n",fNPID) ;
365   parList+=onePar ;
366   sprintf(onePar,"Cuts: \n") ;
367   parList+=onePar ;
368   sprintf(onePar,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
369   parList+=onePar ;
370   sprintf(onePar,"Calorimeter: %s \n",fCalorimeter.Data()) ;
371   parList+=onePar ;
372   sprintf(onePar,"Number of modules: %d \n",fNModules) ;
373   parList+=onePar ;
374   
375   TObjString *oString= new TObjString(parList) ;
376   outputContainer->Add(oString);
377   
378   return outputContainer;
379 }
380
381 //_________________________________________________________________________________________________________________________________________________
382 void AliAnaPi0::Print(const Option_t * /*opt*/) const
383 {
384   //Print some relevant parameters set for the analysis
385   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
386   AliAnaPartCorrBaseClass::Print(" ");
387
388   printf("Number of bins in Centrality:  %d \n",fNCentrBin) ;
389   printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
390   printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
391   printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
392   printf("Number of different PID used:  %d \n",fNPID) ;
393   printf("Cuts: \n") ;
394   printf("Z vertex position: -%2.3f < z < %2.3f \n",fZvtxCut,fZvtxCut) ;
395   printf("Number of modules:             %d \n",fNModules) ;
396   printf("Select pairs with their angle: %d \n",fUseAngleCut) ;
397   printf("------------------------------------------------------\n") ;
398
399
400 //____________________________________________________________________________________________________________________________________________________
401 Int_t AliAnaPi0::GetModuleNumber(AliAODPWG4Particle * particle) 
402 {
403   //Get the EMCAL/PHOS module number that corresponds to this particle
404   
405   Int_t absId = -1;
406   if(fCalorimeter=="EMCAL"){
407     GetReader()->GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
408     if(GetDebug() > 2) 
409       printf("AliAnaPi0::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
410              particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId));
411     return GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId) ;
412   }//EMCAL
413   else{//PHOS
414     Int_t    relId[4];
415     if(!strcmp((GetReader()->GetInputEvent())->GetName(),"AliESDEvent"))   {
416       AliESDCaloCluster *cluster = ((AliESDEvent*)GetReader()->GetInputEvent())->GetCaloCluster(particle->GetCaloLabel(0));
417       if ( cluster->GetNCells() > 0) {
418         absId = cluster->GetCellAbsId(0);
419         if(GetDebug() > 2) 
420           printf("AliAnaPi0::GetModuleNumber(ESD) - PHOS: cluster eta %f, phi %f, e %f, e cluster %f, absId %d\n",
421                  particle->Eta(), particle->Phi()*TMath::RadToDeg(), particle->E(), cluster->E(), absId);
422       }
423       else return -1;
424     }//ESDs
425     else{
426       AliAODCaloCluster *cluster = ((AliAODEvent*)GetReader()->GetInputEvent())->GetCaloCluster(particle->GetCaloLabel(0));
427       if ( cluster->GetNCells() > 0) {
428         absId = cluster->GetCellAbsId(0);
429         if(GetDebug() > 2) 
430           printf("AliAnaPi0::GetModuleNumber(AOD) - PHOS: cluster eta %f, phi %f, e %f, e cluster %f, absId %d\n",
431                  particle->Eta(), particle->Phi()*TMath::RadToDeg(), particle->E(), cluster->E(), absId);
432       }
433       else return -1;
434     }//AODs
435     
436     if ( absId >= 0) {
437       GetReader()->GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
438       if(GetDebug() > 2) 
439         printf("PHOS: Module %d\n",relId[0]-1);
440       return relId[0]-1;
441     }
442     else return -1;
443   }//PHOS
444   
445   return -1;
446 }
447
448 //____________________________________________________________________________________________________________________________________________________
449 void AliAnaPi0::MakeAnalysisFillHistograms() 
450 {
451   //Process one event and extract photons from AOD branch 
452   // filled with AliAnaPhoton and fill histos with invariant mass
453   
454   //Apply some cuts on event: vertex position and centrality range  
455   Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
456   if(IsBadRun(iRun)) return ;   
457   
458   Double_t vert[]={0,0,0} ; //vertex ;
459   GetReader()->GetVertex(vert);
460   if(vert[2]<-fZvtxCut || vert[2]> fZvtxCut) return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
461   
462   //Get Centrality and calculate centrality bin
463   //Does not exist in ESD yet???????
464   Int_t curCentrBin=0 ;
465   
466   //Get Reaction Plain position and calculate RP bin
467   //does not exist in ESD yet????
468   Int_t curRPBin=0 ;    
469   
470   Int_t curZvertBin=(Int_t)(0.5*fNZvertBin*(vert[2]+fZvtxCut)/fZvtxCut) ;
471   
472   fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
473   
474   Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
475   if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
476   Int_t module1 = -1;
477   Int_t module2 = -1;
478   for(Int_t i1=0; i1<nPhot-1; i1++){
479     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
480     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
481     //Get Module number
482     module1 = GetModuleNumber(p1);
483     for(Int_t i2=i1+1; i2<nPhot; i2++){
484       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
485       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
486       //Get module number
487       module2 = GetModuleNumber(p2);
488       Double_t m  = (photon1 + photon2).M() ;
489       Double_t pt = (photon1 + photon2).Pt();
490       Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
491       if(GetDebug() > 2)
492         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
493                p1->Pt(), p2->Pt(), pt,m,a);
494                                 
495       //Check if opening angle is too large or too small compared to what is expected   
496       Double_t angle   = photon1.Angle(photon2.Vect());
497       //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
498       //printf("angle %f\n",angle);
499       if(fUseAngleCut && angle < 0.1) continue;
500       fhRealOpeningAngle   ->Fill(pt,angle);
501       fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
502       
503       //Fill module dependent histograms
504       //if(module1==module2) printf("mod1 %d\n",module1);
505       if(module1==module2 && module1 >=0 && module1<fNModules)
506         fhReMod[module1]->Fill(pt,a,m) ;
507       
508       for(Int_t ipid=0; ipid<fNPID; ipid++)
509         {
510           if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ 
511             fhRe1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
512             if(p1->DistToBad()>0 && p2->DistToBad()>0){
513               fhRe2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
514               if(p1->DistToBad()>1 && p2->DistToBad()>1){
515                 fhRe3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
516               }
517             }
518           }
519         } 
520     }
521   }
522   
523   //Fill mixed
524   
525   TList * evMixList=fEventsList[curCentrBin*fNZvertBin*fNrpBin+curZvertBin*fNrpBin+curRPBin] ;
526   Int_t nMixed = evMixList->GetSize() ;
527   for(Int_t ii=0; ii<nMixed; ii++){  
528     TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
529     Int_t nPhot2=ev2->GetEntriesFast() ;
530     Double_t m = -999;
531     if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d\n", ii, nPhot);
532     
533     for(Int_t i1=0; i1<nPhot; i1++){
534       AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
535       TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
536       for(Int_t i2=0; i2<nPhot2; i2++){
537         AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
538         
539         TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
540         m =           (photon1+photon2).M() ; 
541         Double_t pt = (photon1 + photon2).Pt();
542         Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
543         
544         //Check if opening angle is too large or too small compared to what is expected
545         Double_t angle   = photon1.Angle(photon2.Vect());
546         //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
547         if(fUseAngleCut && angle < 0.1) continue;  
548         
549         if(GetDebug() > 2)
550           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
551                  p1->Pt(), p2->Pt(), pt,m,a);                   
552         for(Int_t ipid=0; ipid<fNPID; ipid++){ 
553           if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ 
554             fhMi1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
555             if(p1->DistToBad()>0 && p2->DistToBad()>0){
556               fhMi2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
557               if(p1->DistToBad()>1 && p2->DistToBad()>1){
558                 fhMi3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
559               }
560               
561             }
562           }
563         }
564       }
565     }
566   }
567   
568   TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
569   //Add current event to buffer and Remove redandant events 
570   if(currentEvent->GetEntriesFast()>0){
571     evMixList->AddFirst(currentEvent) ;
572     currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
573     if(evMixList->GetSize()>=fNmaxMixEv)
574       {
575         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
576         evMixList->RemoveLast() ;
577         delete tmp ;
578       }
579   } 
580   else{ //empty event
581     delete currentEvent ;
582     currentEvent=0 ; 
583   }
584   
585   //Acceptance
586   if(IsDataMC() && GetReader()->ReadStack()){   
587     AliStack * stack = GetMCStack();
588     if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
589       for(Int_t i=0 ; i<stack->GetNprimary(); i++){
590         TParticle * prim = stack->Particle(i) ;
591         if(prim->GetPdgCode() == 111){
592           Double_t pi0Pt = prim->Pt() ;
593           //printf("pi0, pt %2.2f\n",pi0Pt);
594           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception          
595           Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
596           Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
597           if(TMath::Abs(pi0Y) < 0.5){
598             fhPrimPt->Fill(pi0Pt) ;
599           }
600           fhPrimY  ->Fill(pi0Y) ;
601           fhPrimPhi->Fill(phi) ;
602           
603           //Check if both photons hit Calorimeter
604           Int_t iphot1=prim->GetFirstDaughter() ;
605           Int_t iphot2=prim->GetLastDaughter() ;
606           if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
607             TParticle * phot1 = stack->Particle(iphot1) ;
608             TParticle * phot2 = stack->Particle(iphot2) ;
609             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
610               //printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
611               //        phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
612               
613               TLorentzVector lv1, lv2;
614               phot1->Momentum(lv1);
615               phot2->Momentum(lv2);
616               
617               Bool_t inacceptance = kFALSE;
618               if(fCalorimeter == "PHOS"){
619                 if(GetReader()->GetPHOSGeometry()){
620                   Int_t mod ;
621                   Double_t x,z ;
622                   if(GetReader()->GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetReader()->GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x)) 
623                     inacceptance = kTRUE;
624                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
625                 }
626                 else{
627                   
628                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter)) 
629                     inacceptance = kTRUE ;
630                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
631                 }
632                 
633               }    
634               else if(fCalorimeter == "EMCAL"){
635                 if(GetReader()->GetEMCALGeometry()){
636                   if(GetReader()->GetEMCALGeometry()->Impact(phot1) && GetReader()->GetEMCALGeometry()->Impact(phot2)) 
637                     inacceptance = kTRUE;
638                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
639                 }
640                 else{
641                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter)) 
642                     inacceptance = kTRUE ;
643                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
644                 }
645               }   
646               
647               if(inacceptance){
648                 
649                 fhPrimAccPt->Fill(pi0Pt) ;
650                 fhPrimAccPhi->Fill(phi) ;
651                 fhPrimAccY->Fill(pi0Y) ;
652                 Double_t angle  = lv1.Angle(lv2.Vect());
653                 fhPrimOpeningAngle   ->Fill(pi0Pt,angle);
654                 fhPrimCosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
655                 
656               }//Accepted
657             }// 2 photons      
658           }//Check daughters exist
659         }// Primary pi0
660       }//loop on primaries      
661     }//stack exists and data is MC
662   }//read stack
663   else if(GetReader()->ReadAODMCParticles()){
664     if(GetDebug() >= 0)  printf("AliAnaPi0::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
665   }     
666   
667 }       
668
669 //________________________________________________________________________
670 void AliAnaPi0::ReadHistograms(TList* outputList)
671 {
672   // Needed when Terminate is executed in distributed environment
673   // Refill analysis histograms of this class with corresponding histograms in output list. 
674   
675   // Histograms of this analsys are kept in the same list as other analysis, recover the position of
676   // the first one and then add the next.
677   Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
678   
679   if(!fhRe1) fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
680   if(!fhRe2) fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
681   if(!fhRe3) fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
682   if(!fhMi1) fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
683   if(!fhMi2) fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
684   if(!fhMi3) fhMi3 = new TH3D*[fNCentrBin*fNPID] ;      
685   if(!fhReMod) fhReMod = new TH3D*[fNModules] ; 
686   
687   for(Int_t ic=0; ic<fNCentrBin; ic++){
688     for(Int_t ipid=0; ipid<fNPID; ipid++){
689       fhRe1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
690       fhMi1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
691       fhRe2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
692       fhMi2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
693       fhRe3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
694       fhMi3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
695     }
696   }
697   
698   fhEvents = (TH3D *) outputList->At(index++); 
699   
700   //Histograms filled only if MC data is requested      
701   if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
702     fhPrimPt     = (TH1D*)  outputList->At(index++);
703     fhPrimAccPt  = (TH1D*)  outputList->At(index++);
704     fhPrimY      = (TH1D*)  outputList->At(index++);
705     fhPrimAccY   = (TH1D*)  outputList->At(index++);
706     fhPrimPhi    = (TH1D*)  outputList->At(index++);
707     fhPrimAccPhi = (TH1D*)  outputList->At(index++);
708   }
709   
710   for(Int_t imod=0; imod < fNModules; imod++)
711     fhReMod[imod] = (TH3D*) outputList->At(index++);
712   
713 }
714
715
716 //____________________________________________________________________________________________________________________________________________________
717 void AliAnaPi0::Terminate(TList* outputList) 
718 {
719   //Do some calculations and plots from the final histograms.
720   
721   printf(" *** %s Terminate:\n", GetName()) ; 
722   
723   //Recover histograms from output histograms list, needed for distributed analysis.    
724   ReadHistograms(outputList);
725   
726   if (!fhRe1) {
727     printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
728     return;
729   }
730   
731   printf("AliAnaPi0::Terminate()         Mgg Real        : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(),   fhRe1[0]->GetRMS() ) ;
732   
733   char nameIM[128];
734   sprintf(nameIM,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
735   TCanvas  * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
736   cIM->Divide(2, 2);
737   
738   cIM->cd(1) ; 
739   //gPad->SetLogy();
740   TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ(Form("IMPtAll_%s",fCalorimeter.Data()));
741   hIMAllPt->SetLineColor(2);
742   hIMAllPt->SetTitle("No cut on  p_{T, #gamma#gamma} ");
743   hIMAllPt->Draw();
744
745   cIM->cd(2) ; 
746   TH3F * hRe1Pt5 = (TH3F*)fhRe1[0]->Clone(Form("IMPt5_%s",fCalorimeter.Data()));
747   hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
748   TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
749   hIMPt5->SetLineColor(2);  
750   hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
751   hIMPt5->Draw();
752   
753   cIM->cd(3) ; 
754   TH3F * hRe1Pt10 =  (TH3F*)fhRe1[0]->Clone(Form("IMPt10_%s",fCalorimeter.Data()));
755   hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
756   TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
757   hIMPt10->SetLineColor(2);  
758   hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
759   hIMPt10->Draw();
760   
761   cIM->cd(4) ; 
762   TH3F * hRe1Pt20 =  (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
763   hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
764   TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
765   hIMPt20->SetLineColor(2);  
766   hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
767   hIMPt20->Draw();
768    
769   char nameIMF[128];
770   sprintf(nameIMF,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
771   cIM->Print(nameIMF);
772
773   char namePt[128];
774   sprintf(namePt,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
775   TCanvas  * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
776   cPt->Divide(2, 2);
777
778   cPt->cd(1) ; 
779   //gPad->SetLogy();
780   TH1D * hPt = (TH1D*) fhRe1[0]->Project3D("x");
781   hPt->SetLineColor(2);
782   hPt->SetTitle("No cut on  M_{#gamma#gamma} ");
783   hPt->Draw();
784
785   cPt->cd(2) ; 
786   TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
787   hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
788   TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
789   hPtIM1->SetLineColor(2);  
790   hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
791   hPtIM1->Draw();
792   
793   cPt->cd(3) ; 
794   TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
795   hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
796   TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
797   hPtIM2->SetLineColor(2);  
798   hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
799   hPtIM2->Draw();
800
801   cPt->cd(4) ; 
802   TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
803   hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
804   TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
805   hPtIM3->SetLineColor(2);  
806   hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
807   hPtIM3->Draw();
808    
809   char namePtF[128];
810   sprintf(namePtF,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
811   cPt->Print(namePtF);
812
813  
814   char line[1024] ; 
815   sprintf(line, ".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ; 
816   gROOT->ProcessLine(line);
817   sprintf(line, ".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data()); 
818   gROOT->ProcessLine(line);
819  
820   printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
821
822 }
823
824
825
826
827