void plot_beam_slope(){
  double x[400],dx[400],y[400],dy[400];
  double sigx[400],sigy[400],dsigx[400],dsigy[400],z[400];
 
  TF1 *f1=new TF1("f1","gaus(0)+gaus(3)",-10,10);
  f1->SetParameters(1000,0,0.2,100,0,1.);
  
  unsigned int first_bin=35,last_bin=110,bin_range=last_bin-first_bin;
  for (int i=first_bin;i<last_bin+1;i++){
    TH1D *projx=TwoTrackXYZ->ProjectionX ("px",0,-1,i+1,i+1);
    f1->SetParameter(0,projx->GetMaximum());
    f1->SetParameter(1,0.);
    f1->SetParameter(2,0.25);
    f1->SetParameter(3,0.1*projx->GetMaximum());
    f1->SetParameter(4,0);
    f1->SetParameter(5,0.5);
    projx->Fit("f1","q"); 
    x[i]=f1->GetParameter(1);
    dx[i]=f1->GetParError(1); 
    sigx[i]=f1->GetParameter(2);
    dsigx[i]=f1->GetParError(2);
    TH1D *projy=TwoTrackXYZ->ProjectionY("py",0,-1,i+1,i+1);
    f1->SetParameter(0,projy->GetMaximum());
    f1->SetParameter(1,0);
    f1->SetParameter(2,0.25);
    f1->SetParameter(3,0.1*projy->GetMaximum());
    f1->SetParameter(4,0);
    f1->SetParameter(5,0.5);
    projy->Fit("f1","q"); 
    y[i]=f1->GetParameter(1);
    dy[i]=f1->GetParError(1);
    sigy[i]=f1->GetParameter(2);
    dsigy[i]=f1->GetParError(2);
    z[i]=0.5*double(i)+30.25;

    //cout <<" x " <<  x[i] << " y " << y[i] << "z " << z[i] << endl;
    
  }
  //  return;
  TCanvas *c2=new TCanvas();
  c2->cd();
  TH2F *h1=new TH2F("h1","Beam position as function of z",140,30,100,100,-0.5,0.5);
  h1->SetXTitle("z [cm]");
  h1->SetYTitle("vertex x or y [cm]");
  h1->Draw();

  TGraphErrors *g1=new TGraphErrors(bin_range,&z[first_bin],&x[first_bin],0,&dx[first_bin]);
  g1->SetMarkerStyle(20);
  g1->Draw("p");
  g1->Fit("pol1","q");
  TF1 *fx=g1->GetFunction("pol1");
  
  TGraphErrors *g2=new TGraphErrors(bin_range,&z[first_bin],&y[first_bin],0,&dy[first_bin]);
  g2->SetMarkerStyle(20);
  g2->SetMarkerColor(2);
  g2->SetLineColor(2);
  g2->Draw("p");
  g2->Fit("pol1","q");
  TF1 *fy=g2->GetFunction("pol1");
  fy->SetLineColor(2);
  
  TLegend *legend=new TLegend(0.15,0.7,0.25,0.8);
  legend->AddEntry(g1,"x","p");
  legend->AddEntry(g2,"y","p");
  legend->Draw();

  TCanvas *c3=new TCanvas();
  c3->cd();
  
  TH2F *h2=new TH2F("h2","Beam spot size as function of z",140,30,100,100,-0.5,0.5);
  h2->SetXTitle("z [cm]");
  h2->SetYTitle("vertex #sigma_{x} or #sigma_{y} [cm]");
  h2->Draw();
  TGraphErrors *g3=new TGraphErrors(bin_range,&z[first_bin],&sigx[first_bin],0,&dsigx[first_bin]);
  g3->SetMarkerStyle(20);
  g3->Draw("p");
  g3->Fit("pol1","q");
  
  TGraphErrors *g4=new TGraphErrors(bin_range,&z[first_bin],&sigy[first_bin],0,&dsigy[first_bin]);
  g4->SetMarkerStyle(20);
  g4->SetMarkerColor(2);
  g4->SetLineColor(2);
  g4->Draw("p");
  g4->Fit("pol1","q");
  g4->GetFunction("pol1")->SetLineColor(2);

  TLegend *legend2=new TLegend(0.15,0.7,0.25,0.8);
  legend2->AddEntry(g3,"#sigma_{x}","p");
  legend2->AddEntry(g4,"#sigma_{y}","p");
  legend2->Draw();

  double mysigx=g3->GetFunction("pol1")->Eval(65.);
  double mysigy=g4->GetFunction("pol1")->Eval(65.);
  double myx=fx->Eval(65.);
  double myy=fy->Eval(65.);
  double myxslope=fx->GetParameter(1);
  double myyslope=fy->GetParameter(1);
  cout<<"Beam position at z=65 cm: ("<<setprecision(4)<<myx<<","<<myy<<")"<<endl;
  cout<<"Beam spot size at z=65 cm: ("<<mysigx<< ","<< mysigy << ")"<< endl;
  cout << "Beam slopes: (" <<myxslope<< ","<<myyslope<<")"<<endl;
  
  ofstream myfile("beamspot.dat");
  myfile << myx << " " << myy << " 65. " << mysigx*mysigx << " 0. " 
	 << mysigy*mysigy << " 0. 0. 0. " << myxslope << " 0. " << myyslope 
	 <<endl;
  myfile.close();
}