//--------------------------------------------------------------------+
//								      |
//  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.        |
//								      |
//  Since the histograms are identified internally by the title       |
//  string and not by the ID, the title strings must be unique.	      |
//								      |
//  Note that this program runs in several modes:                     |
//								      |
//  Type ./hepToyz -hbook to run the program with the HBook file      |
//			manager. This is also the default. The        |
//			output file is named hepToyz.rz               |
//                                                                    |
//  Type ./hepToyz -histo to run the program with the HistoScope file |
//            		manager. This is currently disabled. The      |
//			output file is named hepToyz.hs               |
//                                                                    |
//  Type ./hepToyz -root to run the program with the Root file        |
//			manager. The output file is named             |
//			hepToyz.root				      |
//                                                                    |
//	Note. The program determines if the appropriate output file   |
//            exists already. It so, it will be updated. If not, it   |
//            is created ab initio.                                   |
//                                                                    |
//--------------------------------------------------------------------+

#include "ZMutility/ZMenvironment.h"

#include "PhysicsVectors/PhysicsVectors.h"
ZM_USING_NAMESPACE( zmpv )

#ifdef USE_HBOOK
#include "HepTuple/HepHBookFileManager.h"
#endif
#ifdef USE_HISTOSCOPE
#include "HepTuple/HepHistoFileManager.h"
#endif
#ifdef USE_ROOT
#include "HepTuple/HepRootFileManager.h"
#endif

#include "HepTuple/HepHist1D.h"
#include "HepTuple/HepHist2D.h"
#include "HepTuple/HepHistProf.h"
#include "HepTuple/HepNtuple.h"
ZM_USING_NAMESPACE( zmht )

#include "CLHEP/Random/Randomize.h"
#include "CLHEP/Units/PhysicalConstants.h"

#include "ZMutility/FixedTypes.h"
#include "ZMutility/iostream"

#include "ZMutility/fstream"
ZM_USING_NAMESPACE( zmfxt )

#include "ZMutility/iostream"
USING( std::cerr )
USING( std::endl )

#include <cstdlib>
#include <cmath>
#ifndef ISOCXX__ISOCXX
  #include <STLUtility/fwd_vector.h>
#endif
#include <vector>
USING( std::vector )
#include <string>
USING( std::string )

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

HepFileManager *tManager;

// Define a handy function to compute the energy of a particle for a given
// momentum and mass.

Float8 energy( SpaceVector p, Float8 mass ) {
  return sqrt((p.dot(p)+mass*mass)); }

int main(Int4 argc, char** argv) {
  Int4 nevents=2500;
  bool first;
  void generate( Int4, bool, string, string );
  bool newFile( string );

//  Sort out which flavor of manager the user wants.
//  Default to HBook if nothing is said.

  string topdir = "zippo";
  string mgr_type = "HBOOK";
  if( argc > 1 )  {
    if( string(argv[argc-1]) == string("-histo")) {
      mgr_type = "HISTO";
      cerr << "\nSorry. HistoScope can't do ntuples just yet." << endl;
      exit(-1);
    }
    if( string(argv[argc-1]) == string("-root"))  mgr_type = "ROOT";
  }

//  Initialize the PhysicsVectors exception handler.

//ZMxPhysicsVectors::setHandler ( ZMexIgnoreAlways() );

//  Initialize a manager and open an output file.  Note that this is now the
//  only direct invocation of an manager-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 output 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.

  if( mgr_type == "HBOOK" ) {
    topdir = "hepToyz.rz";
    first = newFile( topdir );
    int blocksize = 4096;
    tManager = new HepHBookFileManager( topdir, blocksize );
#ifdef USE_HISTOSCOPE
  } else if ( mgr_type == "HISTO" ) {
    topdir = "hepToyz.hs";
    first = newFile( topdir );
    tManager = new HepHistoFileManager( topdir );
#endif
#ifdef USE_ROOT
  } else if ( mgr_type == "ROOT") {
    topdir = "hepToyz.root";
    first = newFile( topdir );
    tManager = new HepRootFileManager( topdir );
#endif
  } else {
    cerr << "\nCannot determine file manager type\n" << endl;
    exit(-1);
  }

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

  cerr << "Calling generate function" << endl;
  generate(nevents,first,topdir,mgr_type);
  delete tManager;
  return 0;
}

void generate(Int4 nevents, bool first, string topdir, string mgr_type)
{
  const Int4 maxTracks = 12;
  Int4   numTracks;
  Int4   sumTracks=0;
  Int4   sumPairs=0;
  Int4   StartSeed=5736321;
  Int4   LuxLevel=3;
  Int4   numDim=1;
  Float4 W=1.0;
  Float8 Etamax=5.0;
  Float8 piMass=0.13956995;
  Float8 KMass=0.493677;
  Float8 ResMass=50.00;
  Float8 ResWidth=8.0;
  Float8 Ptmean = 10.0;
  Float8 Pt, Eta, Phi, Eta1, Eta2;
  Float8 Pmag, Px, Py, Pz, Cost;
  Float8 PtArray[maxTracks];
  Float8 EtaArray[maxTracks];
  Float8 PhiArray[maxTracks];
  Float8 PxArray[maxTracks];
  Float8 PyArray[maxTracks];
  Float8 PzArray[maxTracks];
  Float8 defaultValue = -999.9;

  RanluxEngine    someEngine;
  RandFlat        Rflt(someEngine);
  RandGauss       Rgaus(someEngine);
  RandExponential Rexp(someEngine);
  RandBreitWigner RBW(someEngine);
  void pxyz(Float8, Float8, RandFlat &, Float8 *, Float8 *, Float8 *, Float8 *);

  cerr << "\n\tFunction generate called to produce " << nevents << " events." << endl;

//  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.
//
//  All the uproar with " p = first ? tManager-> ... : tManager-> ..."
//  is a trick to switch between a new run and a continuation run on the
//  basis of a single switch variable and still keep everything in scope.
//  A simple
//            if(first) {
//              HepHist1D & ppt = tManager-> hist1D(...
//                     :
//            } else {
//              HepHist1D & ppt = tManager-> retrieveHist1D(...
//                     :
//            }
//  doesn't do it because the { }'s define a scope for the declarations
//  inside, but it's the wrong scope. Consequently, we resort to this
//  monstrosity.

  cerr << "\nGenerate some fake events and populate some histograms." << endl;

  if( first )  tManager->mkdir("junk");
  tManager->cd("junk");

  HepHist1D &  ppt= first
    ? tManager->hist1D("P_{t} distribution",       50,   0.0f, 50.0f)
    : tManager->retrieveHist1D("P_{t} distribution" );
  HepHist1D & peta= first
    ? tManager->hist1D("#eta distribution",        60,  -6.0f,  6.0f)
    : tManager->retrieveHist1D("#eta distribution"  );
  HepHist1D & pphi= first
    ? tManager->hist1D("#Phi distribution",        64,  -3.2f,  3.2f)
    : tManager->retrieveHist1D("#Phi distribution"  );
  HepHist1D & pmom= first
    ? tManager->hist1D("P distribution",           50,   0.0f,100.0f)
    : tManager->retrieveHist1D("P distribution"     );
  HepHist1D &  ppx= first
    ? tManager->hist1D("P_{x} distribution",       50, -25.0f, 25.0f)
    : tManager->retrieveHist1D("P_{x} distribution" );
  HepHist1D &  ppy= first
    ? tManager->hist1D("P_{y} distribution",       50, -25.0f, 25.0f)
    : tManager->retrieveHist1D("P_{y} distribution" );
  HepHist1D &  ppz= first
    ? tManager->hist1D("P_{z} distribution",       50, -50.0f, 50.0f)
    : tManager->retrieveHist1D("P_{z} distribution" );
  HepHist1D & pcos= first
    ? tManager->hist1D("Cos( #theta )",    50,  -1.0f,  1.0f)
    : tManager->retrieveHist1D("Cos( #theta )" );

// Testing the new title setting member function
  if( first && mgr_type == "ROOT" )
    pcos.setTitle(string("Substitute Hist1D title string for Cos #theta"));

  HepHist1D & pmas= first
    ? tManager->hist1D("K pi Effective Mass",   50,   0.0f,100.0f)
    : tManager->retrieveHist1D("K pi Effective Mass"   );
  HepHistProf & prof= first
    ? tManager->histProf("Mass(Cos #theta) Profile",50,   0.0f,100.0f, -1.0f,  1.0f," ")
    : tManager->retrieveHistProf("Mass(Cos #theta) Profile" );

// Testing the new title setting member function
  if( first && mgr_type == "ROOT" )
    prof.setTitle(string("Substitute HistProf title string for Mass(Cos #theta )"));

  HepHist2D & psct= first
    ? tManager->hist2D("Mass vs cos #theta",  50,   0.0f,100.0f, 50,  -1.0f,  1.0f)
    : tManager->retrieveHist2D("Mass vs cos #theta");

// Testing the new title setting member function
  if( first && mgr_type == "ROOT" )
    psct.setTitle(string("Substitute Hist2D title string for Mass vs cos( #theta )"));

  HepHist1D & pntk= first
    ? tManager->hist1D("Number of tracks per event",20,   0.0f, 20.0f)
    : tManager->retrieveHist1D("Number of tracks per event");

//  Instantiate an ntuple, declare it to be columnwise (this is default
//  anyhow) and tell the manager where to find the data items when we tell
//  it to capture. This is a handy way to do this.  We establish a reference
//  to the data item here and give the item a value when the program logic
//  warrants it. After all items have been assigned a value, we tell the
//  ntuple machinery to grab everything. This works because the location
//  of the data doesn't change - just the value sitting in that location.
//  There's a lot of action here, so we walk through it carefully.

//  The whole business with the number and size of buffers is currently up
//  for discussion so we don't mess with any of that here. It's too much like
//  asking for trouble and we have enough without asking for more.

  HepNtuple & tupl = first
    ? tManager->ntuple("Kinematics CWN")
    : tManager->retrieveNtuple("Kinematics CWN");

  HepHist1D &  cmppx= first
    ? tManager->hist1D("CM pion Px distribution", 50, -50.0f, 50.0f)
    : tManager->retrieveHist1D("CM pion Px distribution" );
  HepHist1D &  cmppy= first
    ? tManager->hist1D("CM pion Py distribution", 50, -50.0f, 50.0f)
    : tManager->retrieveHist1D("CM pion Py distribution" );
  HepHist1D &  cmppz= first
    ? tManager->hist1D("CM pion Pz distribution", 50, -50.0f, 50.0f)
    : tManager->retrieveHist1D("CM pion Pz distribution" );

  tupl.setColumnWise();
  tupl.setDiskResident();

//  If this is the first of a series of runs, establish the random number
//  seed.  Otherwise, restore the state as of the end of the previous run.

  if( first ) {
    someEngine.setSeed( StartSeed, LuxLevel );
    someEngine.showStatus();
  } else {
    someEngine.restoreStatus("hepToyz.config");
    someEngine.showStatus();
  }

// Set some axis labels. For HBook and Histo, these functions are dummies.
// For Root, they insert text into the TAxis object.
  ppx.setXTitle("p_{x} (Gev)");
  ppx.setYTitle("dN/dp_{x}");
  ppy.setXTitle("p_{y} (Gev)");
  ppy.setYTitle("dN/dp_{y}");
  ppz.setXTitle("p_{z} (Gev)");
  ppz.setYTitle("dN/dp_{z}");
  ppt.setXTitle("p_{t} (Gev)");
  ppt.setYTitle("dN/dp_{t}");
  pphi.setXTitle("#Phi");
  pphi.setYTitle("dN/d#Phi");
  prof.setXTitle("K-#pi Mass (Gev)");
  prof.setYTitle("Cos #theta ");
  psct.setXTitle("K-#pi Mass (GeV)");
  psct.setYTitle("cos #theta ");
  psct.setZTitle("d^{2}N/dM d#theta ");
  peta.setXTitle("#eta");
  peta.setYTitle("dN/d#eta ");

  if( first ) {
    cerr << "\nHistograms and one CWNtuple instantiated." << endl;
  } else {
    cerr << "\nHistograms and one CWNtuple re-instantiated." << endl;
  }

//  N. B. In addition to defining the ntuple block names and column names,
//        these declaration also indirectly declare the data types of the
//        corresponding items, either through the type of the data item,
//        the type of the default value or the type of the referenced item.
//        This is fairly subtle and needs some care.

  if( first ) {
//  Define a block named Tracks and a column named Ntracks.  In addition,
//  declare this column to contain an index variable and specify it's range.

    tupl.columnAt("TRACKS::Ntracks",&numTracks,(Int4)1);
    tupl.span("TRACKS::Ntracks",0,maxTracks);	    // specify nametag

//  Declare a column array named Pt and say that it is a variable length
//  array indexed by Ntracks and has maximum dimension maxTracks.  Repeat
//  this for each column within the Tracks block. Note that we specifically
//  cast the default value as a Float8.  This fixes the type of the data.

//  There are a number of variations on how to declare indexed column arrays.
//  We use several of them here to demonstrate some of the possibilities.

    tupl.columnArray("TRACKS::Pt", numDim, defaultValue);
    tupl.dimension("TRACKS::Pt",maxTracks);	// specify nametag
    tupl.index("TRACKS::Pt","Ntracks");

    tupl.columnArray("TRACKS::Eta",numDim, defaultValue);
    tupl.dimension(maxTracks);			// default nametag to "current"
    tupl.index("TRACKS::Eta","Ntracks");

    tupl.columnArray("TRACKS::Phi",numDim, defaultValue);
    tupl.dimension(maxTracks);
    tupl.index("TRACKS::Phi","Ntracks");		// specify nametag

    tupl.columnArray("TRACKS::Px", numDim, defaultValue);
    tupl.dimension(maxTracks);
    tupl.index("Ntracks");			// default nametag to "current"

// Two-fold method cascade since each returns another reference

    tupl.columnArray("TRACKS::Py", numDim, defaultValue).dimension(maxTracks);
    tupl.index("TRACKS::Py","Ntracks");

// Three-fold method cascade since each returns another reference

    tupl.columnArray("TRACKS::Pz", numDim, defaultValue).dimension(maxTracks).index("Ntracks");
  } 

//  For a re-instantiation, most of the needed information is obtained from
//  the file itself.  However, we always need to tell the code where to find
//  the data items and what default values to use. That information is not
//  and, in the case of the pointers to the data, can not be recovered from
//  the file. 

  Float8* currentDefault = new Float8;
  tupl.columnAt("TRACKS::Pt", &PtArray[0] );
  if( tupl.columnDefault("TRACKS::Pt", currentDefault ) ) {
    cerr << "\nDefault value for TRACKS::Pt  has been previously set to " 
         << *currentDefault << endl;
  } else {
    cerr << "\nNo default value has been defined for TRACKS::Pt.  Set it now to " 
         << defaultValue << endl;
    tupl.setColumnDefault( "TRACKS::Pt", defaultValue );
  }

  tupl.columnAt("TRACKS::Eta",&EtaArray[0] );
  if( tupl.columnDefault("TRACKS::Eta", currentDefault ) ) {
    cerr << "Default value for TRACKS::Eta has been previously set to " 
         << *currentDefault << endl;
  } else {
    cerr << "No default value has been defined for TRACKS::Eta. Set it now to " 
         << defaultValue << endl;
    tupl.setColumnDefault( "TRACKS::Eta", defaultValue );
  }

  tupl.columnAt("TRACKS::Phi",&PhiArray[0] );
  if( tupl.columnDefault("TRACKS::Phi", currentDefault ) ) {
    cerr << "Default value for TRACKS::Phi has been previously set to " 
         << *currentDefault << endl;
  } else {
    cerr << "No default value has been defined for TRACKS::Phi. Set it now to " 
         << defaultValue << endl;
    tupl.setColumnDefault( "TRACKS::Phi", defaultValue );
  }

  tupl.columnAt("TRACKS::Px", &PxArray[0] );
  if( tupl.columnDefault("TRACKS::Px", currentDefault ) ) {
    cerr << "Default value for TRACKS::Px  has been previously set to " 
         << *currentDefault << endl;
  } else {
    cerr << "No default value has been defined for TRACKS::Px.  Set it now to " 
         << defaultValue << endl;
    tupl.setColumnDefault( "TRACKS::Px", defaultValue );
  }

  tupl.columnAt("TRACKS::Py", &PyArray[0] );
  if( tupl.columnDefault("TRACKS::Py", currentDefault ) ) {
    cerr << "Default value for TRACKS::Py  has been previously set to " 
         << *currentDefault << endl;
  } else {
    cerr << "No default value has been defined for TRACKS::Py.  Set it now to " 
         << defaultValue << endl;
    tupl.setColumnDefault( "TRACKS::Py", defaultValue );
  }

  tupl.columnAt("TRACKS::Pz", &PzArray[0] );
  if( tupl.columnDefault("TRACKS::Pz", currentDefault ) ) {
    cerr << "Default value for TRACKS::Pz  has been previously set to " 
         << *currentDefault << endl;
  } else {
    cerr << "No default value has been defined for TRACKS::Pz.  Set it now to " 
         << defaultValue << endl;
    tupl.setColumnDefault( "TRACKS::Pz", defaultValue );
  }

//  Finally, print out the block descriptor string that this produced.

  cerr << "\nPrint the chForm string for Block TRACKS to see what we got." << endl;
  cerr << tupl.blockFormat("TRACKS") << endl;

//  Release everything so we can define values for all quantities.

  tupl.clearDataBlock("TRACKS");

//  Put a useless histogram in another subdirectory and make some entries.
//  This is just to verify that all that machinery works as advertised.

  if( first )  tManager->mkdir("trash");
  tManager->cd("trash");
  HepHist1D &  xxx = first
    ? tManager->hist1D( "Gaussian distribution", 60, -6.0f, 6.0f, 10 )
    : tManager->retrieveHist1D("Gaussian distribution");
  for( Int4 l = 0; l < 10000; ++l ) {
    xxx.accumulate((Float4)Rgaus.fire(0.5,1.25),W);
  }
  xxx.release();
  tManager->cd( string("..") );

  string dirlist =  tManager->ls("//"+topdir,"r");
  cerr << "\nMake a full directory listing\n" << dirlist << endl;

//  Establish a list of tracks (more accurately, a list of pointers to tracks).

  vector<SpaceVector *> Tracks;

//  Make some "events". Generate pseudo tracks using random numbers.

  for( Int4 j=0; j<nevents; ++j ) {

//  Generate some pseudo tracks.

//  For the current "event" choose a number of tracks between 2 and maxTracks-2.
//  We're going to throw in the last two tracks later.

    numTracks = (Int4)Rflt.fireInt(2,maxTracks-2);
    pntk.accumulate((Float4)numTracks,W);
    sumTracks += numTracks;

//  For each track, generate a pt and an eta. Transform those into
//  (px,py,pz) and use them to make a three vector. Stick a pointer
//  to the final SpaceVector into our track list.

    for( Int4 k=0; k<numTracks; ++k )  {
      Pt = 0.0;
      while (Pt < 0.01) {
        Pt = Rexp.fire( Ptmean );
      }
      ppt.accumulate((Float4)Pt,W);

//  Use a flat eta distribution for the most part. Just for grins, once
//  in a while, throw in a variation:
//    Draw two values of |eta| from a flat distribution, keep the
//    larger, put a sign on it and histogram the result. The resulting
//    distribution is a bit less boring.

      if( j%8 == 0 )  {
        Eta1=Rflt.fire(Etamax);
        Eta2=Rflt.fire(Etamax);
        Eta=Eta1;
        if(Eta2 > Eta ) Eta = Eta2;
        if( Rflt.fire(-0.5,0.5) < 0.0 ) Eta=-Eta;
      } else {
        Eta=Rflt.fire(2.0*Etamax)-Etamax;
      }
      peta.accumulate((Float4)Eta,W);

//  Compute the three vector components, put them in a SpaceVector, add
//  it to the list, 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, Rflt, &Phi, &Px, &Py, &Pz);
      SpaceVector * pvec = new SpaceVector(Px, Py, Pz);
      Tracks.push_back(pvec);

//  Catch the data for the ntuple arrays and store it away.

      PtArray[k] = Pt;
      EtaArray[k] = Eta;
      PhiArray[k] = Phi;
      PxArray[k] = Px;
      PyArray[k] = Py;
      PzArray[k] = Pz;

//  For the last two tracks in this set, we assume it is an unstable particle
//  that decays into a pion and a kaon, work out the kinematics of the
//  decay and add the decay products to the list.

     if( k == numTracks-2 ) {
       double Mres = RBW.fire( ResMass, ResWidth, 3.0*ResWidth );
       SpaceVector Pres( Px, Py, Pz );
       double Eres = energy( Pres, Mres );
       LorentzVector CMres( Pres, Tcomponent( Eres ));
       SpaceVector beta = Pres/Eres;
       CMres.boost( beta );
       double Epi = ( Mres*Mres + piMass*piMass - KMass*KMass )/(2.0*Mres);
       double Ppi = sqrt( Epi*Epi - piMass*piMass );
       double cosTheta = Rflt.fire( -1.0, 1.0 );
       double Pzpi = Ppi*cosTheta;
       double Ptpi = sqrt( Ppi*Ppi - Pzpi*Pzpi );
       double CMPhi = Rflt.fire( -pi, pi );
       double Pxpi = Ptpi*cos( CMPhi );
       double Pypi = sqrt( Ptpi*Ptpi - Pxpi*Pxpi );
       if( Rflt.fire( -1.0, 1.0 ) < 0.0 ) Pypi = -Pypi;
       LorentzVector CMpi( Pxpi, Pypi, Pzpi, Tcomponent( Epi ));
       cmppx.accumulate((Float4)Pxpi,W);
       cmppy.accumulate((Float4)Pypi,W);
       cmppz.accumulate((Float4)Pzpi,W);
       double EK = ( Mres*Mres + KMass*KMass - piMass*piMass )/(2.0*Mres);
       LorentzVector CMK( -Pxpi, -Pypi, -Pzpi, Tcomponent( EK ));
       CMpi.boost( -beta );
       CMK.boost( -beta );
       SpaceVector * Labpi = new SpaceVector(CMpi.getV() );
       SpaceVector * LabK = new SpaceVector(CMK.getV() );
       Tracks.push_back( Labpi );
       Tracks.push_back( LabK );
     }
   }

//  Finished making tracks for this event. Make some histogram and ntuple
//  entries. Capture a bunch of ntuple entries before they go out of scope.
//  Note that we explicitly specify the event-by-event length for the
//  indexed arrays.

    tupl.capture("TRACKS::Ntracks",numTracks);
    tupl.capture("TRACKS::Pt" ,&(PtArray[0]),  numTracks);
    tupl.capture("TRACKS::Eta",&(EtaArray[0]), numTracks);
    tupl.capture("TRACKS::Phi",&(PhiArray[0]), numTracks);
    tupl.capture("TRACKS::Px", &(PxArray[0]),  numTracks);
    tupl.capture("TRACKS::Py", &(PyArray[0]),  numTracks);
    tupl.capture("TRACKS::Pz", &(PzArray[0]),  numTracks);
    tupl.storeCapturedData();
    tupl.clearData();
    if (j==0) {
      cerr << "\nPrint the (possibly) reordered chForm string for Block TRACKS." << endl;
      cerr << tupl.blockFormat("TRACKS") << endl;
    }

//  Run an iterator over the track list and plot things that depend
//  only on the properties of that one track.

    typedef vector<SpaceVector *> sp_vector;
    for( sp_vector::const_iterator iter(Tracks.begin());
         iter != Tracks.end();
         iter++) {
      SpaceVector *ptrk = *iter;
      Pmag=ptrk->mag();
      Px=ptrk->x();
      Py=ptrk->y();
      Pz=ptrk->z();
      pphi.accumulate((Float4)Phi,W);
      pmom.accumulate((Float4)Pmag,W);
      ppx.accumulate((Float4)Px,W);
      ppy.accumulate((Float4)Py,W);
      ppz.accumulate((Float4)Pz,W);
    }

//  Now do something that depends on more than one track. Compute the
//  opening angle between all distinct track pairs and histogram it.

//  Run two iterators over the same list of tracks and sift out the distinct
//  track pairs.  This is a little messy since the list is really a list of
//  pointers to tracks. Hence the iterator is a pointer to a pointer and
//  that's why we have to dereference it twice to get a track.

    for( sp_vector::const_iterator ptrk1(Tracks.begin());
         ptrk1 != Tracks.end();
         ptrk1++) {
      UnitVector p1Hat(**ptrk1);

//  Establish another iterator and loop over the same list, looking
//  for a distinct "track 2".

//  Run the second iterator only until the two pointers are equal. This
//  should prevent generating pairs that are permutations of each other.
//  Note that this works only because a vector container is ordered.

      for( sp_vector::const_iterator ptrk2(Tracks.begin());
           ptrk2 != ptrk1;
           ptrk2++) {
        Float8 Mass;
        sumPairs++;
        UnitVector p2Hat(**ptrk2);
        Cost = p1Hat.dot(p2Hat);
        pcos.accumulate((Float4)Cost,W);

//  Assuming one track is a pion and the other a K, compute the effective mass.

        LorentzVector v1( **ptrk1, Tcomponent( energy( **ptrk1, piMass )));
        LorentzVector v2( **ptrk2, Tcomponent( energy( **ptrk2,  KMass )));
        Mass = v1.invariantMass(v2);
        pmas.accumulate((Float4)Mass,W);
        psct.accumulate((Float4)Mass,(Float4)Cost,W);
        prof.accumulate((Float4)Mass,(Float4)Cost,W);

//  Switch the mass assignments, recompute Mass and enter it again

        v1.set( **ptrk1, Tcomponent( energy( **ptrk1,  KMass )));
        v2.set( **ptrk2, Tcomponent( energy( **ptrk2, piMass )));
        Mass = v1.invariantMass(v2);
        pmas.accumulate((Float4)Mass,W);
        psct.accumulate((Float4)Mass,(Float4)Cost,W);
        prof.accumulate((Float4)Mass,(Float4)Cost,W);
      }
    }

//  Clear out the track list (but don't delete it) and go do it all again.
//  Recall we instantiated the track vectors on the free store so the
//  SpaceVector destructor doesn't see them. It is our responsibility to make
//  them go away.

//  We don't need to delete the other SpaceVectors.  They were instantiated
//  on the stack and are going out of scope here. In that case, the destructor
//  takes care of them.

    for( sp_vector::const_iterator iter( Tracks.begin() );
                                   iter != Tracks.end();
                                   iter++ ) {
      delete *iter;
    }
    Tracks.clear( );
  }
  tupl.release();

//  Finished generating tracks. Summarize, clean up and return.

  cerr << "\nTotal number of tracks generated: " << sumTracks << endl;
  cerr << "\nTotal number of tracks pairs:     " << sumPairs  << endl;

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

  cerr << "\nHistograms filled.  Write the output file." << endl;

//  Save the state of the RandomEngine

  someEngine.saveStatus("hepToyz.config");
  tManager->cd( string("//")+topdir );

  return;
}

  void pxyz(Float8 Pt, Float8 Eta, RandFlat & Rand, Float8 *pphi, Float8 *ppx, Float8 *ppy, Float8 *ppz)
{

//  Given pt and eta, compute px, py and pz.  This transformation
//  is unique up to the decomposition of pt into px and py.

//  Note that we had to pass a reference to the random number stream
//  to this function. Otherwise, the function wouldn't know anything
//  about random numbers. The point is that the generator is NOT global.
//  It's just another object to be passed around.

    Float8 Sign, absEta, Theta, absTheta, Pmag, Phi;
    Float8 Px, Py;

    Sign = 1.0;
    if(Eta < 0.0) Sign = -1.0;
    absEta = Sign*Eta;
    absTheta = exp(-absEta);
    Theta = 2.0*Sign*atan(absTheta);
    Pmag = Pt/sin(Theta);
    *ppz = Sign*sqrt(Pmag*Pmag-Pt*Pt);
    Phi = Rand.fire(-pi,pi);
    *pphi = Phi;
    Px = cos(Phi)*Pt;
    *ppx = Px;
    Py = sqrt(Pt*Pt-Px*Px);
    if( Rand.fire(-0.5,0.5) > 0.0 )  Py = -Py;
    *ppy = Py;
    return;
}
bool newFile( string topdir ) {

//  Sort out if the file exists already. Try to open it
//  and see if the attempt fails.

  std::ifstream test(topdir.c_str(),std::ios::in);
  if( test ) {
    test.close();
    return false;
  } else {
    return true;
  }
}