//  R. Shanks      February (?) 1999
//  J. Bogart      31 Mar 1999   Changed all float to double

#include "MCDiag.h"
//#ifdef WIN32
//#define hepevt_ HEPEVT
//#endif

struct hepevt HEPEVT;
  //////////////////////////////////////////////////////////////////////////
  //                                                                      //
  // MCDiag                                                               //
  //                                                                      //
  // The MCDiag class produces a stream of N particles with their         //
  // momentum spread out in magnitude and direction. The resulting        //
  // particles are output to a file in StdHep format.                     //
  //                                                                      //
  //////////////////////////////////////////////////////////////////////////

Int_t MCDiag::m_known_parts[] = {11,13,111,211,2212,2112,321,310,130};
Double_t MCDiag::m_known_masses[] = {0.000511,0.105658,0.1349739,0.1395675,
				    0.938272,0.939566,0.493646,0.497671,
				    0.497671};
Int_t MCDiag::m_num_known_parts = 9;

ClassImp(MCDiag)

///_________________________________________________________________________
 MCDiag::MCDiag(char* ofile):m_MagSpread(0.0),m_CosThSpread(0.0),m_PhiSpread(0.0),m_NumParticles(1), m_xpoint(0),m_ypoint(0),m_zpoint(0),m_prod_time(0)
{
  // Default constructor
  m_mag = 5.00; // GeV
  m_phi = TMath::Pi()/4;
  m_costheta = 0.707;
  m_type = 13; // Muon by default
  hepevt_.nevhep = 0;
  factor = new TRandom();

  ilbl = 1;
  istr = 0; 

  if(StdHepXdrWriteInit(ofile,"Robs file", 1, istr)){
    printf("MCDiag: Failed to initialize output file "%s".n", ofile);
    exit(1);
  }// All StdHep routines return non zero if an error has occured	
  
}
///_________________________________________________________________________
 void MCDiag::Set_Seed(UInt_t seed){
  // Set the seed of the random number generator. If zero then seed
  // is set to current machine clock
  factor->SetSeed(seed);
}
///_________________________________________________________________________
 void MCDiag::Set_Spread(Double_t magspread, Double_t costhspread, Double_t phispread){
  // The user is prompted for the momentum spread parameters
  m_MagSpread = magspread;
  m_CosThSpread = costhspread;
  m_PhiSpread = phispread;
}
///_________________________________________________________________________
 void MCDiag::Set_Momentum(Double_t mag, Double_t costheta, Double_t phi){
  // Sets the mean components of the momentum vectors
  m_mag = mag;
  m_phi = phi;
  m_costheta = costheta;
}
///_________________________________________________________________________

 void MCDiag::Set_Start_Point(Double_t xpoint, Double_t ypoint, Double_t zpoint){
  // Sets the starting point of the particles(s).
  // Default is (0,0,0) aka the IP.
  m_xpoint = xpoint;
  m_ypoint = ypoint;
  m_zpoint = zpoint;
}
///_________________________________________________________________________
 void MCDiag::Set_Production_Time(Double_t  time){
  // Sets the time of production of the particles(s).
  // Default is 0.0
	m_prod_time = time;
}
///_________________________________________________________________________
 void MCDiag::Set_Particles(Int_t num){
  // Sets the number of particles
  m_NumParticles = num;
}
///_________________________________________________________________________
 void MCDiag::Set_Type(Int_t type){
  // Sets the particle type. NOTE: Must use StdHep numbering scheme.Particle 
  // mass is set from type
  m_type = type;
  for(Int_t i =0;i<m_num_known_parts;i++){
    if((m_known_parts[i] == type)||(m_known_parts[i] == -type)){
      m_mass = m_known_masses[i]; // mass
      return;
    }
    else continue;
  }
  printf("Particle type %i unknown.n Supported types include:n",type);
  printf("11 electronn13 muonn111 pi-zeron211 pi-plusn2112 neutronn2212 protonn321 k-plusn310 k0-shortn130 k0-long");
  exit(1); 
      
}
///_________________________________________________________________________
 Int_t MCDiag::getEvent(Event* event){

  Int_t ierr;
  if(ierr = doit()){
    printf("MCDiag::Error in filling common block.");
    return(1);
  }
  readStdFile* rsf = new readStdFile();
  if(ierr = rsf->makeMCs(event)){
    printf("MCDiag::Error in filling the event structure.");
    return(1);
  }
  return(0);
}
///_________________________________________________________________________
 Int_t MCDiag::doit(){
  // Create the stream of Particles
  hepevt_.nevhep += 1;
  hepevt_.nhep = m_NumParticles;
  for(Int_t i =0;i<m_NumParticles; i++){
    Double_t newmag = m_mag + (2*factor->Rndm()-1)*(m_MagSpread);
    Double_t newcth = m_costheta +(2*factor->Rndm()-1)*(m_CosThSpread);
    Double_t sinth = sqrt(1.-newcth*newcth);
    Double_t newphi = m_phi + (2*factor->Rndm()-1)*(m_PhiSpread);
    if(newphi < 0) newphi += 2*TMath::Pi();
    if(newphi > 2*TMath::Pi()) newphi -= 2*TMath::Pi();

    hepevt_.isthep[i] = 1;                // Status code
    hepevt_.idhep[i] = m_type;            // Particle ID
    hepevt_.jmohep[i][0] = 0;             // index of mother
    hepevt_.phep[i][0] = newmag*sinth*TMath::Cos(newphi);// x momentum
    hepevt_.phep[i][1] = newmag*sinth*TMath::Sin(newphi);// y momentum
    hepevt_.phep[i][2] = newmag*newcth;// z momentum
    hepevt_.phep[i][4] = m_mass; 
    hepevt_.phep[i][3] = sqrt(newmag*newmag + 
			      hepevt_.phep[i][4]*hepevt_.phep[i][4]);
    hepevt_.vhep[i][0] = m_xpoint; // x position
    hepevt_.vhep[i][1] = m_ypoint; // y position
    hepevt_.vhep[i][2] = m_zpoint; // z position
    hepevt_.vhep[i][3] = m_prod_time; // production time
  }
  return(0);
}
///_________________________________________________________________________
 Int_t MCDiag::spew(char* ofile)const {
    // Save the particle data into the hepevt common block
    // Write the common block out ot file using StdHepXdrWrite()

  if(StdHepXdrWrite(ilbl,istr)){
    printf("Write failed.n");
    exit(1);
  }
  StdHepZero();
    return(0);
}
///________________________________________________________________________
 void MCDiag::cleanup(){
  // Closes the output stream.
    StdHepXdrEnd(istr);
}