Upload New File

CODE/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),
+	initx(0.*m),
+	inity(0.*m),
+	initz(-200.*m),
+	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);
+	particleDirection
+		= G4ThreeVector(0., 0.,1.);
+	
+	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);
+
+		}
+	}
+}
+