Jeng/scripts/qcd.C

#define qcd_cxx
#include "qcd.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TMultiGraph.h>
#include <TGraphErrors.h>
#include <TCanvas.h>
#include "Math/GSLIntegrator.h"
#include "Math/WrappedTF1.h"
#include <iostream>
#include <TLegend.h>
#include <TPaveStats.h>

int nbin = 50; // 200 for 1 and 2 jet bin
double xMin = 0.0;
double xMax = 2.0;

// Select fitting function:
TString fitfunc = "landau";  // landau, polX, gaus ...
// Select New or Old reliso:
TString isoname = "New";
// Jet bin:
Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5
double d0sigCut = 3.;

void qcd::Loop()
{
  //   In a ROOT session, you can do:
  //      Root > .L qcd.C++
  //      Root > qcd t
  //      Root > t.GetEntry(12); // Fill t data members with entry number 12
  //      Root > t.Show();       // Show values of entry 12
  //      Root > t.Show(16);     // Read and show values of entry 16
  //      Root > t.Loop();       // Loop on all entries
  
  TH1D *h_RelIso_ttbar = new TH1D("h_RelIso_ttbar","h_RelIso_ttbar",nbin,xMin,xMax);
  TH1D *h_RelIso_qcd = new TH1D("h_RelIso_qcd","h_RelIso_qcd",nbin,xMin,xMax);
  TH1D *h_RelIso_wjets = new TH1D("h_RelIso_wjets","h_RelIso_wjets",nbin,xMin,xMax);
  TH1D *h_RelIso_zjets = new TH1D("h_RelIso_zjets","h_RelIso_zjets",nbin,xMin,xMax);
  TH1D *h_RelIso_all = new TH1D("h_RelIso_all","h_RelIso_all",nbin,xMin,xMax);
    
  // define fit region:
  double AX0 = 0.;
  double AXf = 0.;
  double BX0 = 0.;
  double BXf = 0.;
  
  cout<<"Use "<<isoname<<" RelIso"<<endl;
  
  if ( isoname.Contains("Old") ) {
    // old
    AX0 = 0.95;
    AXf = 1.0;
    BX0 = 0.2;
    BXf = 0.8;
  }
  else if ( isoname.Contains("New") ) {
    // new
    AX0 = 0.0;
    AXf = 0.053;
    BX0 = 0.15;
    BXf = 0.6;
  }
  //   // old weights
  //   double wttbar = 0.0094;
  //   double wqcd = 0.3907;
  //   double wwjets = 0.0177;
  //   double wzjets = 0.041;
  
  // weights
  double wttbar = 0.0101;
  double wqcd   = 1.3161;
  double wwjets = 0.0977;
  double wzjets = 0.078;
  
  // Na = num of qcd in signal region
  double Na,Nb,Nfac,Ns, Ns2, NsErrP, NsErrM;
  Na = Nb = Nfac = Ns = Ns2 = NsErrP = NsErrM = 0;
  
  double Nttbar, NWjets, NZjets, Nqcd;
  Nttbar = NWjets = NZjets = Nqcd = 0;
  double reliso = 0.0;
  
  if (fChain == 0) return;
  
  Long64_t nentries = fChain->GetEntriesFast();
  
  Long64_t nbytes = 0, nb = 0;
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    // if (Cut(ientry) < 0) continue;
    bool signalRegion = false;
    bool fitRegion = false;
    
    TString filename(fChain->GetCurrentFile()->GetName());
    if (isoname.Contains("Old")) {
      reliso = top_muon_old_reliso[0];
      if ( reliso > AX0 ) signalRegion = true;
      else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
    }
    else if (isoname.Contains("New")) {
      reliso = top_muon_new_reliso[0];
      if ( reliso < AXf ) signalRegion = true;
      else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
    }
    else {
      cout<<"No RelIso defined!!!"<<endl;
      return;
    }
    
    double d0sig = TMath::Abs(top_muon_d0[0])/top_muon_d0Error[0];
    bool correctBin = false; // jet bin condition
    bool goodD0sig = false; // d0sig condition

    if ( jbin > 4 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
    if (top_njets == jbin && jbin != 4) correctBin = true;
    else if (top_njets >= jbin && jbin == 4) correctBin = true;
    if ( d0sig < d0sigCut ) goodD0sig = true;
    
    // Jet bin
    if( correctBin && goodD0sig ) {
      if ( filename.Contains("TTJets") ) {
        h_RelIso_ttbar->Fill(reliso);
        if ( signalRegion ) { Nttbar += wttbar; }
        else if ( fitRegion ) Nfac += wttbar;
      }
      if ( filename.Contains("MuPt15") ) {
        h_RelIso_qcd->Fill(reliso);
        if ( signalRegion ) { Na += wqcd; Nqcd+= wqcd; } // QCD in signal region
        else if ( fitRegion ) {
          Nb += wqcd; // QCD in fitted background region
          Nfac += wqcd;
        }
      }
      if ( filename.Contains("WJets") ) {
        h_RelIso_wjets->Fill(reliso);
        if ( signalRegion ) { NWjets += wwjets; }       
        else if ( fitRegion ) Nfac += wwjets;
      }
      if ( filename.Contains("ZJets") ) {
        h_RelIso_zjets->Fill(reliso);
        if ( signalRegion ) { NZjets += wzjets; }
        else if ( fitRegion ) Nfac += wzjets;
      }
    }
  }
  cout << "\n";
  cout << " Nb = " << Nb << "; Nfac = " << Nfac << endl << endl;
  cout << " N ttbar = " << Nttbar << endl;
  cout << " N qcd = " << Nqcd << endl;
  cout << " N Wjets = " << NWjets << endl;
  cout << " N Zjets = " << NZjets << endl << endl;
  
  gStyle->SetOptFit(01111);
  
  TCanvas *cv1 = new TCanvas("cv1","cv1",800,600);
  
  h_RelIso_qcd->Scale(wqcd);
  h_RelIso_wjets->Scale(wwjets);
  h_RelIso_zjets->Scale(wzjets);
  h_RelIso_ttbar->Scale(wttbar);

  h_RelIso_ttbar->SetXTitle("Combined Relative Isolation");
  h_RelIso_ttbar->SetMarkerColor(2);
  h_RelIso_ttbar->SetLineColor(2);   
  h_RelIso_ttbar->Draw();
  h_RelIso_all->Add(h_RelIso_ttbar,1.);
  h_RelIso_all->Add(h_RelIso_qcd,1.);
  h_RelIso_all->Add(h_RelIso_wjets,1.);
  h_RelIso_all->Add(h_RelIso_zjets,1.);
  h_RelIso_all->SetLineWidth(2);
  if( jbin == 4)
    h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 4 jets)");
  else if ( jbin < 4 ) {
    TString title = isoname+" RelIso distribution (";
    title += jbin;
    title += " jets)";
    h_RelIso_all->SetTitle(title);
  }
  else {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
  h_RelIso_all->SetMinimum(0);
  //   h_RelIso_all->Sumw2();
  //   h_RelIso_all->Draw("e1");
  h_RelIso_all->Draw("same");

  // Fitting here:
  FitBKG(fitfunc,h_RelIso_all,h_RelIso_qcd,cv1,AX0,AXf,BX0,BXf,Na,Nfac);
  
}

void  qcd::FitBKG(TString s0,TH1D *h0,TH1D *h1,TCanvas *cv0,
                  double ax0,double axf,double bx0,double bxf,double na,double nfac) {
  
  double ns, ns2;
  ns = ns2 = 0;
  
  // Fit QCD in background region
  cout<<"Fitting with "<<s0<<" function!"<<endl;
  TH1D *h_all = (TH1D*)h0->Clone();
  TH1D *h_qcd = (TH1D*)h1->Clone();
  
  h_all->Fit(s0,"","",bx0,bxf);
  
  //Get a pointer to the fitted function
  TF1 *f1 = (TF1*)h_all->GetFunction(s0);
  Int_t npar = f1->GetNpar();
  //   Double_t chi2 = f1->GetChisquare();
  double * par0 = new double [npar];
  double * parErr0 = new double [npar];
  printf("Function has %i parameters. Chisquare = %g\n",
         npar,
         f1->GetChisquare());
  for (Int_t i=0;i<npar;i++) {
    printf("%s = %g +- %g\n",
           f1->GetParName(i),
           f1->GetParameter(i),
           f1->GetParError(i)
           );
    par0[i] = f1->GetParameter(i);
    parErr0[i] = f1->GetParError(i);
    //       cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
  }
  f1->SetLineColor(2);
  h_qcd->SetLineColor(40);
  h_qcd->SetFillStyle(3018);
  h_qcd->SetFillColor(38);
  h_qcd->Draw("same");
  f1->Draw("same");
    
  TF1 *f2 = (TF1*)f1->Clone();
  f2->SetRange(ax0,axf);
  f2->SetLineColor(4);
  f2->Draw("same");
  int np = 100;
  double *x=new double[np];
  double *w=new double[np];
  f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
  ns = f2->IntegralFast(np,x,w,ax0,axf);
  ns/=(f2->Integral(bx0,bxf));
  ns*=nfac;
  cout<<"Ns by usual integral = "<<ns<<endl;
  delete [] x;
  delete [] w;
  
  // Another way to do integration
  TF1 *g = (TF1*)f1->Clone();
  ROOT::Math::GSLIntegrator ig(1.E-8,1.E-8,1000); 
  ROOT::Math::WrappedTF1 wf(*g);
  ig.SetFunction(wf);
  ns2 = ig.Integral(ax0,axf);
  ns2/=(g->Integral(bx0,bxf));
  ns2*=nfac;
  cout<<"Ns by MathMore integral = "<<ns2<<endl;
    
  // Yet another estimation
  double sum=0.,limitLow=ax0,limitHigh=axf,vx=0.,dx=0.;
  dx=(limitHigh-limitLow)/nbin;
  //   cout<<"limitHigh = "<<limitHigh<<"; dx = "<<dx<<endl;
  for(int i=1;i<nbin;i++){
    vx=h_all->GetXaxis()->GetBinCenter(i);
    if(vx <= limitHigh){
      //       cout<<"f2("<<vx<<")= "<<f2->Eval(vx)<<endl;
      sum+=f2->Eval(vx);
    }
  }
  cout<<"sum = "<<sum<<endl;
  cv0->Update();
  
  TPaveStats *p1 = (TPaveStats*)h_all->GetListOfFunctions()->FindObject("stats");
  p1->SetTextColor(kBlue);
  p1->SetX1NDC(0.6);
  p1->SetX2NDC(0.88);
  p1->SetY1NDC(0.62);
  p1->SetY2NDC(0.88);
  p1->Draw();
  
  // Label histo
  TLegend * leg1 = new TLegend(0.55,0.35,0.88,0.60);
  leg1->SetHeader("Fit with "+fitfunc);
  leg1->AddEntry(h_all,"All events (S+B)");
  leg1->AddEntry(h_qcd,"QCD events");
  leg1->AddEntry(f1,"Fit to all events in control region");
  leg1->AddEntry(f2,"Extrapolation of QCD events in signal region");
  leg1->Draw();
  // Print results
  cout << " Observed Ns = " << ns;
  cout <<"\n";
  cout << " Expected Na = " << na << endl;
}
Revision:	1.1
Committed:	Tue Apr 14 21:46:09 2009 UTC (16 years ago) by jengbou
Content type:	text/plain
Branch:	MAIN
Log Message:	QCD extrapolation fitting
#	User	Rev	Content
1	jengbou	1.1	#define qcd_cxx
2			#include "qcd.h"
3			#include <TH2.h>
4			#include <TStyle.h>
5			#include <TCanvas.h>
6			#include <TMath.h>
7			#include <TMultiGraph.h>
8			#include <TGraphErrors.h>
9			#include <TCanvas.h>
10			#include "Math/GSLIntegrator.h"
11			#include "Math/WrappedTF1.h"
12			#include <iostream>
13			#include <TLegend.h>
14			#include <TPaveStats.h>
15
16			int nbin = 50; // 200 for 1 and 2 jet bin
17			double xMin = 0.0;
18			double xMax = 2.0;
19
20			// Select fitting function:
21			TString fitfunc = "landau"; // landau, polX, gaus ...
22			// Select New or Old reliso:
23			TString isoname = "New";
24			// Jet bin:
25			Int_t jbin = 4; // 1-3, 4 for >=4. Put in number less than 5
26			double d0sigCut = 3.;
27
28			void qcd::Loop()
29			{
30			// In a ROOT session, you can do:
31			// Root > .L qcd.C++
32			// Root > qcd t
33			// Root > t.GetEntry(12); // Fill t data members with entry number 12
34			// Root > t.Show(); // Show values of entry 12
35			// Root > t.Show(16); // Read and show values of entry 16
36			// Root > t.Loop(); // Loop on all entries
37
38			TH1D *h_RelIso_ttbar = new TH1D("h_RelIso_ttbar","h_RelIso_ttbar",nbin,xMin,xMax);
39			TH1D *h_RelIso_qcd = new TH1D("h_RelIso_qcd","h_RelIso_qcd",nbin,xMin,xMax);
40			TH1D *h_RelIso_wjets = new TH1D("h_RelIso_wjets","h_RelIso_wjets",nbin,xMin,xMax);
41			TH1D *h_RelIso_zjets = new TH1D("h_RelIso_zjets","h_RelIso_zjets",nbin,xMin,xMax);
42			TH1D *h_RelIso_all = new TH1D("h_RelIso_all","h_RelIso_all",nbin,xMin,xMax);
43
44			// define fit region:
45			double AX0 = 0.;
46			double AXf = 0.;
47			double BX0 = 0.;
48			double BXf = 0.;
49
50			cout<<"Use "<<isoname<<" RelIso"<<endl;
51
52			if ( isoname.Contains("Old") ) {
53			// old
54			AX0 = 0.95;
55			AXf = 1.0;
56			BX0 = 0.2;
57			BXf = 0.8;
58			}
59			else if ( isoname.Contains("New") ) {
60			// new
61			AX0 = 0.0;
62			AXf = 0.053;
63			BX0 = 0.15;
64			BXf = 0.6;
65			}
66			// // old weights
67			// double wttbar = 0.0094;
68			// double wqcd = 0.3907;
69			// double wwjets = 0.0177;
70			// double wzjets = 0.041;
71
72			// weights
73			double wttbar = 0.0101;
74			double wqcd = 1.3161;
75			double wwjets = 0.0977;
76			double wzjets = 0.078;
77
78			// Na = num of qcd in signal region
79			double Na,Nb,Nfac,Ns, Ns2, NsErrP, NsErrM;
80			Na = Nb = Nfac = Ns = Ns2 = NsErrP = NsErrM = 0;
81
82			double Nttbar, NWjets, NZjets, Nqcd;
83			Nttbar = NWjets = NZjets = Nqcd = 0;
84			double reliso = 0.0;
85
86			if (fChain == 0) return;
87
88			Long64_t nentries = fChain->GetEntriesFast();
89
90			Long64_t nbytes = 0, nb = 0;
91			for (Long64_t jentry=0; jentry<nentries;jentry++) {
92			Long64_t ientry = LoadTree(jentry);
93			if (ientry < 0) break;
94			nb = fChain->GetEntry(jentry); nbytes += nb;
95			// if (Cut(ientry) < 0) continue;
96			bool signalRegion = false;
97			bool fitRegion = false;
98
99			TString filename(fChain->GetCurrentFile()->GetName());
100			if (isoname.Contains("Old")) {
101			reliso = top_muon_old_reliso[0];
102			if ( reliso > AX0 ) signalRegion = true;
103			else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
104			}
105			else if (isoname.Contains("New")) {
106			reliso = top_muon_new_reliso[0];
107			if ( reliso < AXf ) signalRegion = true;
108			else if ( reliso > BX0 && reliso < BXf ) fitRegion = true;
109			}
110			else {
111			cout<<"No RelIso defined!!!"<<endl;
112			return;
113			}
114
115			double d0sig = TMath::Abs(top_muon_d0[0])/top_muon_d0Error[0];
116			bool correctBin = false; // jet bin condition
117			bool goodD0sig = false; // d0sig condition
118
119			if ( jbin > 4 ) {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
120			if (top_njets == jbin && jbin != 4) correctBin = true;
121			else if (top_njets >= jbin && jbin == 4) correctBin = true;
122			if ( d0sig < d0sigCut ) goodD0sig = true;
123
124			// Jet bin
125			if( correctBin && goodD0sig ) {
126			if ( filename.Contains("TTJets") ) {
127			h_RelIso_ttbar->Fill(reliso);
128			if ( signalRegion ) { Nttbar += wttbar; }
129			else if ( fitRegion ) Nfac += wttbar;
130			}
131			if ( filename.Contains("MuPt15") ) {
132			h_RelIso_qcd->Fill(reliso);
133			if ( signalRegion ) { Na += wqcd; Nqcd+= wqcd; } // QCD in signal region
134			else if ( fitRegion ) {
135			Nb += wqcd; // QCD in fitted background region
136			Nfac += wqcd;
137			}
138			}
139			if ( filename.Contains("WJets") ) {
140			h_RelIso_wjets->Fill(reliso);
141			if ( signalRegion ) { NWjets += wwjets; }
142			else if ( fitRegion ) Nfac += wwjets;
143			}
144			if ( filename.Contains("ZJets") ) {
145			h_RelIso_zjets->Fill(reliso);
146			if ( signalRegion ) { NZjets += wzjets; }
147			else if ( fitRegion ) Nfac += wzjets;
148			}
149			}
150			}
151			cout << "\n";
152			cout << " Nb = " << Nb << "; Nfac = " << Nfac << endl << endl;
153			cout << " N ttbar = " << Nttbar << endl;
154			cout << " N qcd = " << Nqcd << endl;
155			cout << " N Wjets = " << NWjets << endl;
156			cout << " N Zjets = " << NZjets << endl << endl;
157
158			gStyle->SetOptFit(01111);
159
160			TCanvas *cv1 = new TCanvas("cv1","cv1",800,600);
161
162			h_RelIso_qcd->Scale(wqcd);
163			h_RelIso_wjets->Scale(wwjets);
164			h_RelIso_zjets->Scale(wzjets);
165			h_RelIso_ttbar->Scale(wttbar);
166
167			h_RelIso_ttbar->SetXTitle("Combined Relative Isolation");
168			h_RelIso_ttbar->SetMarkerColor(2);
169			h_RelIso_ttbar->SetLineColor(2);
170			h_RelIso_ttbar->Draw();
171			h_RelIso_all->Add(h_RelIso_ttbar,1.);
172			h_RelIso_all->Add(h_RelIso_qcd,1.);
173			h_RelIso_all->Add(h_RelIso_wjets,1.);
174			h_RelIso_all->Add(h_RelIso_zjets,1.);
175			h_RelIso_all->SetLineWidth(2);
176			if( jbin == 4)
177			h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 4 jets)");
178			else if ( jbin < 4 ) {
179			TString title = isoname+" RelIso distribution (";
180			title += jbin;
181			title += " jets)";
182			h_RelIso_all->SetTitle(title);
183			}
184			else {cout<<"Don't mess around!!! jbin should be less than 5!"<<endl;return;}
185			h_RelIso_all->SetMinimum(0);
186			// h_RelIso_all->Sumw2();
187			// h_RelIso_all->Draw("e1");
188			h_RelIso_all->Draw("same");
189
190			// Fitting here:
191			FitBKG(fitfunc,h_RelIso_all,h_RelIso_qcd,cv1,AX0,AXf,BX0,BXf,Na,Nfac);
192
193			}
194
195			void qcd::FitBKG(TString s0,TH1D h0,TH1D h1,TCanvas *cv0,
196			double ax0,double axf,double bx0,double bxf,double na,double nfac) {
197
198			double ns, ns2;
199			ns = ns2 = 0;
200
201			// Fit QCD in background region
202			cout<<"Fitting with "<<s0<<" function!"<<endl;
203			TH1D h_all = (TH1D)h0->Clone();
204			TH1D h_qcd = (TH1D)h1->Clone();
205
206			h_all->Fit(s0,"","",bx0,bxf);
207
208			//Get a pointer to the fitted function
209			TF1 f1 = (TF1)h_all->GetFunction(s0);
210			Int_t npar = f1->GetNpar();
211			// Double_t chi2 = f1->GetChisquare();
212			double * par0 = new double [npar];
213			double * parErr0 = new double [npar];
214			printf("Function has %i parameters. Chisquare = %g\n",
215			npar,
216			f1->GetChisquare());
217			for (Int_t i=0;i<npar;i++) {
218			printf("%s = %g +- %g\n",
219			f1->GetParName(i),
220			f1->GetParameter(i),
221			f1->GetParError(i)
222			);
223			par0[i] = f1->GetParameter(i);
224			parErr0[i] = f1->GetParError(i);
225			// cout<<"Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
226			}
227			f1->SetLineColor(2);
228			h_qcd->SetLineColor(40);
229			h_qcd->SetFillStyle(3018);
230			h_qcd->SetFillColor(38);
231			h_qcd->Draw("same");
232			f1->Draw("same");
233
234			TF1 f2 = (TF1)f1->Clone();
235			f2->SetRange(ax0,axf);
236			f2->SetLineColor(4);
237			f2->Draw("same");
238			int np = 100;
239			double *x=new double[np];
240			double *w=new double[np];
241			f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
242			ns = f2->IntegralFast(np,x,w,ax0,axf);
243			ns/=(f2->Integral(bx0,bxf));
244			ns*=nfac;
245			cout<<"Ns by usual integral = "<<ns<<endl;
246			delete [] x;
247			delete [] w;
248
249			// Another way to do integration
250			TF1 g = (TF1)f1->Clone();
251			ROOT::Math::GSLIntegrator ig(1.E-8,1.E-8,1000);
252			ROOT::Math::WrappedTF1 wf(*g);
253			ig.SetFunction(wf);
254			ns2 = ig.Integral(ax0,axf);
255			ns2/=(g->Integral(bx0,bxf));
256			ns2*=nfac;
257			cout<<"Ns by MathMore integral = "<<ns2<<endl;
258
259			// Yet another estimation
260			double sum=0.,limitLow=ax0,limitHigh=axf,vx=0.,dx=0.;
261			dx=(limitHigh-limitLow)/nbin;
262			// cout<<"limitHigh = "<<limitHigh<<"; dx = "<<dx<<endl;
263			for(int i=1;i<nbin;i++){
264			vx=h_all->GetXaxis()->GetBinCenter(i);
265			if(vx <= limitHigh){
266			// cout<<"f2("<<vx<<")= "<<f2->Eval(vx)<<endl;
267			sum+=f2->Eval(vx);
268			}
269			}
270			cout<<"sum = "<<sum<<endl;
271			cv0->Update();
272
273			TPaveStats p1 = (TPaveStats)h_all->GetListOfFunctions()->FindObject("stats");
274			p1->SetTextColor(kBlue);
275			p1->SetX1NDC(0.6);
276			p1->SetX2NDC(0.88);
277			p1->SetY1NDC(0.62);
278			p1->SetY2NDC(0.88);
279			p1->Draw();
280
281			// Label histo
282			TLegend * leg1 = new TLegend(0.55,0.35,0.88,0.60);
283			leg1->SetHeader("Fit with "+fitfunc);
284			leg1->AddEntry(h_all,"All events (S+B)");
285			leg1->AddEntry(h_qcd,"QCD events");
286			leg1->AddEntry(f1,"Fit to all events in control region");
287			leg1->AddEntry(f2,"Extrapolation of QCD events in signal region");
288			leg1->Draw();
289			// Print results
290			cout << " Observed Ns = " << ns;
291			cout <<"\n";
292			cout << " Expected Na = " << na << endl;
293			}