//--------------------------------------------------------------------+
//								      |
//  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 ./FileSwitch -hbook to run the program with the HBook file   |
//			manager. This is also the default. The        |
//			output file is named FileSwitch.rz            |
//                                                                    |
//  Type ./FileSwitch -histo to run the program with the HistoScope   |
//			file manager. This is currently disabled.     |
//			The output file is named FileSwitch.hs        |
//                                                                    |
//  Type ./FileSwitch -root to run the program with the Root file     |
//			manager. The output file is named             |
//			FileSwitch.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/RanluxEngine.h"
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandBreitWigner.h"
#include "CLHEP/Random/RandGauss.h"
#include "CLHEP/Random/RandExponential.h"
#include "CLHEP/Units/PhysicalConstants.h"

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

#include "TObject.h"

#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 )

#include "HepTuple/HepFileSwitch.h"

typedef std::vector<void*> strayList;

class switchHandler : public HepFileSwitch {

public:

  void strayObjectHandler( strayList& strays ) {

//  Rather than do anything dangerous, just list the names of the stray objects

    if ( strays.empty() ) {
      cerr << " \nStray object list is empty" << endl;
    } else {
      cerr << "\nStray object list is NOT empty. The object names are:" << endl;
      for( strayList::const_iterator it(strays.begin());
				     it!=strays.end();
				     it++ ) {
        TObject* obj = (TObject*)(*it);
        cerr << "\t" << obj->GetName() << endl;
      }
    }
  }
};

// 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=15000;
  bool first;
  HepFileManager * man;
  void generate( HepFileManager*, Int4, bool, string );
  bool newFile( string );

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

  int compression = 1;
  string fileName = "zippo.xxx";
  string topDir = "FileSwitch";
  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" ) {
    fileName = "FileSwitch.rz";
    first = newFile( fileName );
    int blocksize = 4096;
    man = new HepHBookFileManager( fileName, blocksize, topDir );
#ifdef USE_HISTOSCOPE
  } else if ( mgr_type == "HISTO" ) {
    fileName = "FileSwitch.hs";
    first = newFile( fileName );
    man = new HepHistoFileManager( fileName );
#endif
#ifdef USE_ROOT
  } else if ( mgr_type == "ROOT") {
    fileName = "FileSwitch.root";
    first = newFile( fileName );
    man = new HepRootFileManager( fileName, topDir, compression );
#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( man, nevents, first, topDir );

  delete man;
  return 0;
}

void generate( HepFileManager* man, Int4 nevents,
               bool first, string topDir)
{
  const Int4 maxTracks = 12;
  Int4   numTracks;
  Int4   sumTracks=0;
  Int4   sumPairs=0;
  Int4   StartSeed=5736329;
  Int4   LuxLevel=3;
  Int4   numDim=1;
  bool	 reset = false;
  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 ? man-> ... : man-> ..."
//  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 = man-> hist1D(...
//                     :
//            } else {
//              HepHist1D & ppt = man-> 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 )  man->mkdir("junk");
  man->cd("junk");

  HepHist1D &  ppt= first 
    ? man->hist1D("Pt distribution",       50,   0.0f, 50.0f)
    : man->retrieveHist1D("Pt distribution"       );
  HepHist1D & peta= first 
    ? man->hist1D("Eta distribution",      60,  -6.0f,  6.0f)
    : man->retrieveHist1D("Eta distribution"      );
  HepHist1D & pphi= first 
    ? man->hist1D("Phi distribution",      64,  -3.2f,  3.2f)
    : man->retrieveHist1D("Phi distribution"      );
  HepHist1D & pmom= first 
    ? man->hist1D("P distribution",        50,   0.0f,100.0f)
    : man->retrieveHist1D("P distribution"        );
  HepHist1D &  ppx= first 
    ? man->hist1D("Px distribution",       50, -25.0f, 25.0f)
    : man->retrieveHist1D("Px distribution"       );
  HepHist1D &  ppy= first 
    ? man->hist1D("Py distribution",       50, -25.0f, 25.0f)
    : man->retrieveHist1D("Py distribution"       );
  HepHist1D &  ppz= first 
    ? man->hist1D("Pz distribution",       50, -50.0f, 50.0f)
    : man->retrieveHist1D("Pz distribution"       );
  HepHist1D & pcos= first 
    ? man->hist1D("Cos(opening angle)",    50,  -1.0f,  1.0f)
    : man->retrieveHist1D("Cos(opening angle)"    );
  HepHist1D & pmas= first 
    ? man->hist1D("K pi Effective Mass",   50,   0.0f,100.0f)
    : man->retrieveHist1D("K pi Effective Mass"   );
  HepHistProf & prof= first 
    ? man->histProf("Mass(Cost) Profile",50,   0.0f,100.0f, -1.0f,  1.0f," ")
    : man->retrieveHistProf("Mass(Cost) Profile"  );
  HepHist2D & psct= first 
    ? man->hist2D("Mass vs cos(Theta)",  50,   0.0f,100.0f, 50,  -1.0f,  1.0f)
    : man->retrieveHist2D("Mass vs cos(Theta)"    );
  HepHist1D & pntk= first 
    ? man->hist1D("Number of tracks per event",20,   0.0f, 20.0f)
    : man->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 
    ? man->ntuple("Kinematics CWN")
    : man->retrieveNtuple("Kinematics CWN");

  HepHist1D &  cmppx= first 
    ? man->hist1D("CM pion Px distribution", 50, -50.0f, 50.0f)
    : man->retrieveHist1D("CM pion Px distribution" );
  HepHist1D &  cmppy= first 
    ? man->hist1D("CM pion Py distribution", 50, -50.0f, 50.0f)
    : man->retrieveHist1D("CM pion Py distribution" );
  HepHist1D &  cmppz= first 
    ? man->hist1D("CM pion Pz distribution", 50, -50.0f, 50.0f)
    : man->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("FileSwitch.state");
//  someEngine.showStatus();
  }

  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 )  man->mkdir("trash");
  man->cd("trash");
  HepHist1D &  xxx = first 
    ? man->hist1D( "Gaussian distribution", 60, -6.0f, 6.0f )
    : man->retrieveHist1D("Gaussian distribution");
  for( Int4 l = 0; l < 9000; ++l ) {
    xxx.accumulate((Float4)Rgaus.fire(0.5,1.25),W);
  }
  if( first ) man->mkdir("bilge");
  man->cd("bilge");
  HepHist1D &  xxx2 = first 
    ? man->hist1D( "Gaussian distribution2", 60, -6.0f, 6.0f )
    : man->retrieveHist1D("Gaussian distribution2");
  for( Int4 l = 0; l < 9000; ++l ) {
    xxx2.accumulate((Float4)Rgaus.fire(-0.5,1.00),W);
  }
  man->cd( string("//"+topDir) );
  HepHist1D &  xxx3 = first 
    ? man->hist1D( "Gaussian distribution3", 60, -6.0f, 6.0f )
    : man->retrieveHist1D("Gaussian distribution3");
  for( Int4 l = 0; l < 9000; ++l ) {
    xxx3.accumulate((Float4)Rgaus.fire(0.0,0.75),W);
  }
  man->cd("junk");

  string dirlist =  man->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=1; 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( );

// Instantiate a spillover file and keep going

    if( j == 6000 ) {
      string nextFileName = topDir + "_2.root";
      cerr << "\nAfter 6000 events, switch to output file "
           << nextFileName << endl;
      bool status = man->switchFile( nextFileName, reset );
    }

// That was so much fun, let's do it again

    if( j == 11000 ) {
      string nextFileName = topDir + "_3.root";
      cerr << "\nAfter 11000 events, switch to output file "
           << nextFileName << endl;
      switchHandler cleanup;
      bool status = man->switchFile( nextFileName, reset, cleanup );
    }
  }
  man->cd("trash");
  for( Int4 l = 0; l < 11000; ++l ) {
    xxx.accumulate((Float4)Rgaus.fire(0.5,1.25),W);
  }
  man->cd("bilge");
  for( Int4 l = 0; l < 11000; ++l ) {
    xxx2.accumulate((Float4)Rgaus.fire(-0.5,1.00),W);
  }
  man->cd( string("//"+topDir) );
  for( Int4 l = 0; l < 11000; ++l ) {
    xxx3.accumulate((Float4)Rgaus.fire(0.0,0.75),W);
  }
  man->cd("junk");

//  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("FileSwitch.state");
//man->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 fileName ) {

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

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