//--------------------------------------------------------------------
//
//  Sample histogramming job that produces the same results as
//  the classic HBOOK and Fortran style (except for the histogram
//  ID's, which now use the default labelling).  
//
//  This approach uses the HepHistogram member functions, which
//  are the generic wrapper functions. Doing things this way is
//  more flexible and does not commit the user to HBOOK style.
//
//  There is one more ugly detail that must be understood.  Since
//  the histograms are identified internally by the title string
//  and not by the ID, the title strings must be unique.
//
//--------------------------------------------------------------------
#include "HepTuple/HepHBookFileManager.h"
#include "CLHEP/OldRandom/Random.h"
#include "CLHEP/Vector/ThreeVector.h"
#include <string>
#include <stdlib.h>
#include <iostream.h>
#include <math.h>
#include "HepTuple/HepHist1D.h"

//  Declare the pointer to the file manager at this level to
//  make it available at file scope.

HepHBookFileManager *tManager;

main()
{
  int lunit=10;
//int nsamp=100;
  int nsamp=50000;
  void generate(long);

//  Initialize HBOOK and open an RZ file.  Note that this is
//  now the only direct invocation of an HBOOK-specific method.
//  Using it here triggers everything else even though all
//  subsequent references below are to things like HepHistogram.

//  Note that the top directory on the RZ file will have the same name as
//  the file itself.  There's no choice in this but it isn't a big deal
//  since the name of the top directory isn't actually recorded anywhere.

  tManager = new HepHBookFileManager("newlib.rz");

//  Book and fill some histograms in a function named generate.

  generate(nsamp);

//  Shut down the corresponding RZ file, clean up and bail out.

  tManager->write();
  delete tManager;
  return 0;
}

void generate(long nsamples)
{
  int j;
  float w=1.0;
  long seed=9736821;
  double ptmean=2.0;
  double etamax=5.0;
  double TwoPi=2.0*M_PI;
  double pt, eta, eta1, eta2;
  double pmag1, pmag2, px1, py1, pz1, px2, py2, pz2, cost;
  void pxyz(double, double, HepRandom &, double *, double *, double *);

  cout << "\tFunction generate called to produce " << nsamples << " events." << endl;

//Instantiate a random number stream and set the initial seed.

  HepRandom(rnd);
  rnd.setSeed(seed);

//  Book a gang of histograms.  This time, we declare them as HepHistograms
//  and let the virtual functions sort out what kind of histograms we really
//  mean to produce.  In this case, we end up with HBookHistograms without
//  actually saying so.

  cout << "\n>>> Populate some histograms." << endl;

  HepHist1D  *ppt(tManager->newHist1D("Pt distribution   ",50,   0.0f, 10.0f));
  HepHist1D *peta(tManager->newHist1D("Eta distribution  ",55,  -5.5f,  5.5f));
  HepHist1D *pmom(tManager->newHist1D("P distribution    ",50,   0.0f,100.0f));
  HepHist1D  *ppx(tManager->newHist1D("Px distribution   ",50, -25.0f, 25.0f));
  HepHist1D  *ppy(tManager->newHist1D("Py distribution   ",50, -25.0f, 25.0f));
  HepHist1D  *ppz(tManager->newHist1D("Pz distribution   ",50,-100.0f,100.0f));
  HepHist1D *pcos(tManager->newHist1D("Cos(opening angle)",50,  -1.0f,  1.0f));

//  Make some entries. Make some pseudo minbias tracks using random numbers. 

  for(j=0;j<nsamples;j++) {

//  Generate a pseudo track.  Draw a value for pt from an
//  exponential distribution and histogram it

    pt = 0.0;
    while (pt < 0.5) {
      pt=rnd.exponential(ptmean);
    }
    ppt->accumulate((float)pt); 
  
//  Draw two values of |eta| from a flat distribution, keep the
//  larger, put a sign on it and histogram the result. The resulting
//  eta distribution is wierd but not boring.

    eta1=rnd.flat(etamax);
    eta2=rnd.flat(etamax);
    eta=eta1;
    if(eta2 > eta ) eta = eta2;
    if( rnd.flat(-1.0,1.0) < 0.0 ) eta=-eta;
    peta->accumulate((float)eta); 

//  Compute the three vector components, create a Hep3Vector, get
//  its magnitude and histogram all that.  The transformation from
//  (pt,eta) to (px,py,pz) is unique up to the decomposition of
//  pt into px and py.
// 

    pxyz( pt, eta, rnd, &px1, &py1, &pz1);
    Hep3Vector pone(px1, py1, pz1);
    pmag1=pone.mag();
    pmom->accumulate((float)pmag1);
    ppx->accumulate((float)px1);
    ppy->accumulate((float)py1);
    ppz->accumulate((float)pz1);

//  Repeat the exercise and generate another pseudotrack

    pt = 0.0;
    while (pt < 0.5) {
      pt=rnd.exponential(ptmean);
    }
    ppt->accumulate((float)pt); 
  
    eta1=rnd.flat(etamax);
    eta2=rnd.flat(etamax);
    eta=eta1;
    if( eta2 > eta ) eta=eta2;
    if( rnd.flat(-1.0,1.0) < 0.0 ) eta=-eta;
    peta->accumulate((float)eta); 

    pxyz( pt, eta, rnd, &px2, &py2, &pz2);
    Hep3Vector ptwo(px2, py2, pz2);
    pmag2=ptwo.mag();
    pmom->accumulate((float)pmag2);
    ppx->accumulate((float)px2);
    ppy->accumulate((float)py2);
    ppz->accumulate((float)pz2);

//  Compute the opening angle between these two "tracks" and
//  histogram it

    Hep3Vector unit1;
    Hep3Vector unit2;
    unit1 = pone.unit();
    unit2 = ptwo.unit();
    cost = unit1.dot(unit2);
    pcos->accumulate((float)cost);
  }

//  Finished. Write out the histograms

  cout << "\n>>> Histograms filled.  Write out the directory." << endl;
  tManager->writeDirectory();

//  Scrog everything in the current directory

  tManager->rm();

//  Make a directory listing.

  cout << "\n>>> Make a full directory listing." << endl;
  tManager->ls("//newlib.rz","T");
  return;
}	  

void pxyz(double pt, double eta, HepRandom &xrand, double *ppx, double *ppy, double *ppz)
{
//
// Given pt and eta, compute px, py and pz.  This transformation
// is unique up to the decomposition of pt into px and py.
// 
   double sign, abseta, rtheta, absrtheta, pmag, rphi;
   double pxx, pyy;
   double TwoPi=2.0*M_PI;
   long seed=3236521;

   sign = 1.0;
   if(eta < 0.0) sign = -1.0;
   abseta = sign*eta;
   absrtheta = exp(-abseta);
   rtheta = 2.0*sign*atan(absrtheta);
   pmag = pt/sin(rtheta);
   *ppz = sign*sqrt(pmag*pmag-pt*pt);
   rphi = xrand.flat(TwoPi);
   pxx = cos(rphi)*pt;
   *ppx = pxx;
   pyy = sqrt(pt*pt-pxx*pxx);
   if(xrand.flat() > 0.5)  pyy = -pyy;
   *ppy = pyy;
   return;
}