#include "G4VUserPrimaryGeneratorAction.hh"

#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "Randomize.hh"

#include <fstream>
using namespace std;

extern float P_GENERATED; // GeV
extern float THETA_GENERATED; // radians
extern float PHI_GENERATED;  // radians
extern float MASS_GENERATED; // GeV

class PiMuG4PrimaryGeneratorAction:public G4VUserPrimaryGeneratorAction{
	public:

		PiMuG4PrimaryGeneratorAction(const char *ptype="pi+"){
		
			double E = G4RandFlat::shoot(0.5,5.0)*GeV;
			E = 2.0*GeV;
		
			// z=-220 is just upstream of FCAL
			// z=-843 is at nominal CPP target position
			G4ThreeVector position(0.0, 0.0, -220.0*cm);
			//G4ThreeVector position(0.0, 0.0, -843.0*cm);
			G4ThreeVector momentumDirection(0.0, 0.0, 1.0);
			G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
			G4ParticleDefinition* particle = particleTable->FindParticle(ptype);
			fParticleGun  = new G4ParticleGun(1);
			fParticleGun->SetParticleDefinition(particle);
			fParticleGun->SetParticleEnergy(E);
			fParticleGun->SetParticlePosition(position);
			fParticleGun->SetParticleMomentumDirection(momentumDirection);
			
			MASS_GENERATED = particle->GetPDGMass()/GeV;
			
			Nevents=0;
			
			ifs = new ifstream("pb_pol70_10days_pi+.ascii");
			if(!ifs->is_open()){
				G4cout<<"Unable to open input file \"pb_pol70_10days_pi+.ascii\"!"<<G4endl;
//				exit(-1);
			}
		}

		void GeneratePrimaries(G4Event* anEvent){
		
//			// Read event from input file
//			//1 pi+ 0.139570 +1 0.045039 -0.060762 1.403137 1.412089
//			int id, charge;
//			string pname;
//			double m, px, py, pz, E;
//			*ifs >> id >> pname >> m >> charge >> px >> py >> pz >> E;
//			if(ifs->eof()){
//				G4cout<<"NO MORE EVENTS IN INPUT FILE!"<<G4endl;
//				return;
//			}
//			
//			double p = sqrt(px*px + py*py + pz*pz);
//			P_GENERATED = p;
//			
//			// Overwrite energy by recalculating using particle mass.
//			// This is to allow muons or pions to be the generated particle,
//			// but have the momentum come from the input file.
//			E = sqrt(P_GENERATED*P_GENERATED + MASS_GENERATED*MASS_GENERATED);
//			
//			// Directions
//			G4ThreeVector dir(px/p, py/p, pz/p);
//			THETA_GENERATED = dir.theta();
//			PHI_GENERATED = dir.phi();
//			
//			fParticleGun->SetParticleEnergy(E*GeV);
//			fParticleGun->SetParticleMomentumDirection(dir);
		
			fParticleGun->GeneratePrimaryVertex(anEvent);
			
			if(++Nevents%100 == 0){
				G4cout<<" "<<Nevents<<" events generated \r";
				G4cout.flush();
			}
		}

	private:

		G4ParticleGun*  fParticleGun;
		int Nevents;
		ifstream *ifs;
};