void bhwAnalMakeClass(const Char_t *fileName) {
  // open ntuple file
  TNtuple *nt = makeBhwNtuple(fileName);
  // make class
  nt->MakeClass("bhwAnalClass");
}


void bhwAnal(const Char_t *fileName) {

  // open ntuple file
  TNtuple *nt = makeBhwNtuple(fileName);
  // create instance of bhwAnalClass
  bhwAnalClass *bhw = new bhwAnalClass(nt);
  // output file
  TFile *output = new TFile("bhwAnal.root","RECREATE");

  // histos
  TH1D *bspread = new TH1D("bspread","",100,0.0,1.1);
  TH1D *bstrahl = new TH1D("bstrahl","",100,0.0,1.1);
  TH1D *isradn  = new TH1D("isradn","",100,0.0,1.1);
  TH1D *total   = new TH1D("total","",100,0.0,1.1);
  TH1D *xreco   = new TH1D("xreco","",100,0.0,1.1);
  TH2D *scat    = new TH2D("scatter","",100,0.0,1.1,100,0.0,1.1);

  // loop over all events
  for(Int_t i=0;i<Int_t(nt->GetEntriesFast());i++) {
    bhw->GetEntry(i);

    TLorentzVector se = TLorentzVector(bhw->px1,bhw->py1,bhw->pz1,bhw->e1);
    TLorentzVector sp = TLorentzVector(bhw->px2,bhw->py2,bhw->pz2,bhw->e2);

    if(bhw->be != 1.0 || bhw->bp != 1.0) continue;

    bspread->Fill(bhw->se*bhw->sp);
    bstrahl->Fill(bhw->be*bhw->bp);
    isradn->Fill(bhw->isr);
    total->Fill(bhw->tot);

    double x = 0;
    if(se->Theta()+sp->Theta() >= TMath::Pi())
      x = sqrt((1/tan(se->Theta()/2))*(1/tan(sp->Theta()/2)));
    else
      x = sqrt((1/tan((TMath::Pi()-se->Theta())/2))*(1/tan((TMath::Pi()-sp->Theta())/2)));

    xreco->Fill(x);
    scat->Fill(bhw->tot,x);

  }

  // write output
  output->Write();

}