PrimarySpectrum.cc

// Geant libraries
//
#include "G4SystemOfUnits.hh"
#include "globals.hh"
#include "g4root.hh"
#include "Randomize.hh"
#include "G4ios.hh"


// Local libraries
//
#include "PrimarySpectrum.hh"
#include "EventAction.hh"


// c++ libraries
//
#include <math.h>
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <vector>
#include <string>
#include <iomanip>
#include <stdio.h>


PrimarySpectrum::PrimarySpectrum():
	pi(3.14159265358979323846),
	Ro(500.*m), // Radio de la fuente
	initx(0.*m),
	inity(0.*m),
	initz(-200.*m), // Altura a la que se encuentra la fuente 
	r(0.*m),
	phi(0.*m),
	rphi(0.),
	signo1(0),
	signo2(0),
	//line(),
	px(0.),
	py(0.),
	pz(0.),
	condicion(0.)
//	inputFile(),
{
	particlePosition 
		= G4ThreeVector(0.*m,0.*m,-200.*m); // Cambiar aqui tambien la posicion de la fuente
	particleDirection
		= G4ThreeVector(0., 0.,1.); // la particula va en la direccion positiva del eje 'z'
	
	char inputFile[] ="salidasimu.shw.bz2"; //name of the file

	openFile(inputFile);


}


PrimarySpectrum::~PrimarySpectrum()
{}


int PrimarySpectrum::openFile(char *nfi)
{
	char tmpc[256]; 

	snprintf(tmpc,256,"bzip2 -d -c %s",nfi); //lee y descomprime 
	inFile = popen(tmpc,"r");

	return 1;
}


void PrimarySpectrum::primaryMomento()
{
	G4int search = 1;
	char line[256];

	while(search)
	{
	
		if(feof(inFile))
		{
			char inputFile[] = "salidasimu.shw.bz2"; //name of file
			openFile(inputFile);
		}

		if(fgets(line,256,inFile))
		{
			if(line[0] != '#')
			{
				sscanf(
					line,
					"%d %lf %lf %lf %lf %lf %lf %d %d %lf %lf %lf\n",
					&crkId,
					&px,
					&py, 
					&pz,
					&x,
					&y,
					&z,
					&shwId,
					&prmId,
					&prmEner,
					&prmThe,
					&prmPhi
				      ); 
	
				px = px*GeV;
				py = py*GeV;
				pz = pz*GeV;
			}
		}

		search = 0;

		switch(crkId) //CONECTION OF CORSIKA AND GEANT4
		{
			case 1:
				parId = "gamma";
				break;
			case 2:
				parId = "e+";
				break;
			case 3:
				parId = "e-";
				break;
			case 5:
				parId = "mu+";
				break;
			case 6:
				parId = "mu-";
				break;
			case 7:
				parId = "pi(1300)0";
				break;
			case 8:
				parId = "pi(1300)+";
				break;
			case 9:
				parId = "pi(1300)-";
				break;
			case 10:
				parId = "kaon0L";
				break;
			case 11:
				parId = "kaon+";
				break;
			case 12:
				parId = "kaon-";
				break;
			case 13:
				parId = "neutron";
				break;
			case 14:
				parId = "proton";
				break;
			case 15:
				parId = "anti_proton";
				break;
			case 25:
				parId = "anti_neutron";
				break;
			default:
				search = 1;
		}
	}
	
	particleDirection
		= G4ThreeVector(px, py, pz);

	//G4cout <<"ID="<< crkId<<G4endl;
	//G4cout<<"Momento del archivo: "<<"px="<<px<<" "<<"py="<<py<<" "<<"pz="<<pz<<G4endl;

}

void PrimarySpectrum::primaryPosition()
{
	G4int ok = 1;

	signo1 = (int) (10*G4UniformRand());

	signo2 = (int) (10*G4UniformRand());

	while(ok)
	{
		r = Ro*G4UniformRand()*pow((-1),signo1); 
		phi = Ro*G4UniformRand()*pow((-1),signo2); 
	
		rphi = sqrt(r*r + phi*phi);

		if(rphi < Ro)
		{
			ok = 0;
			initx = r;
			inity = phi;
		
			particlePosition 
				= G4ThreeVector(initx, inity, initz);

		}
	}
}