Replace construction.cc

GEANT4/construction.cc

@@ -2,6 +2,8 @@
 
 MyDetectorConstruction::MyDetectorConstruction()
 {
+//************************//
+// Acorde al numero de capas que se desea simular, en este caso estamos simulando 5 capas de atmosfera//
     nCols = 5;
     nRows = 5;
 
@@ -14,6 +16,8 @@ MyDetectorConstruction::MyDetectorConstruction()
     DefineMaterials();
 
     isAtmosphere = true;
+//************************//
+//Dimensiones del Mundo y de la Atmosfera // 
 
     xWorld = 1000.*m;
     yWorld = 1000.*m;
@@ -23,73 +27,135 @@ MyDetectorConstruction::MyDetectorConstruction()
 MyDetectorConstruction::~MyDetectorConstruction()
 {}
 
+//***********************************//
+// Definimos los Materiales y el modelo de densidad atmosferica dependiente de la altura //
+
+
 void MyDetectorConstruction::DefineMaterials()
 {
-    G4NistManager *nist = G4NistManager::Instance();
+    G4NistManager* nist = G4NistManager::Instance();
 
     worldMat = nist->FindOrBuildMaterial("G4_AIR");
-
+    
+    //**********************//
+    //Masa Atomica de los Elementos (He agregado Argon) //
+    
     G4double density0 = 1.29*kg/m3;
     G4double aN = 14.01*g/mole;
     G4double aO = 16.00*g/mole;
-
+    G4double aAr = 39.95*g/mole;
+   
+   //************************//
+   //Se crean los elementos presentes en la atmosfera con el formato G4Element("Nombre del elemento", "Abreviacion", numero atomico, masa atomica) //
+    
     G4Element *elN = new G4Element("Nitrogen", "N", 7, aN);
     G4Element *elO = new G4Element("Oxygen", "O", 8, aO);
-
-    G4double f = 3;
-    G4double R = 8.3144626181532;
-    G4double g0 = 9.81;
-    G4double kappa = (f+2)/f;
-    G4double T = 293.15;
-    G4double M = (0.3*14.01 + 0.7*16.0)/1000.;
-
-    for(G4int i = 0; i < 5; i++)
+    G4Element* elAr = new G4Element("Argon", "Ar", 18, aAr);
+
+	//**********************//
+	//Constantes que necesita el modelo de densidad //
+	
+    G4double f = 3; //grados de libertad
+    G4double R = 8.3144626181532; // Constante de los gases
+    G4double g0 = 9.81; // Gravedad
+    G4double kappa = (f + 2) / f; // No se que es pero es importante
+    G4double T = 293.15; // Temperatura estandar al nivel del mar
+    G4double M = (0.23 * 14.01 + 0.76 * 16.0 + 0.01 * 39.95 ) / 1000.; // la masa molar promedio de la mezcla de elementos.
+    
+    
+    //*****************************//
+
+	// El for crea las diferentes densidades del aire, i = 5 es el numero de capas de aire (con diferentes densidades) que deseo crear
+	
+    for (G4int i = 0; i < 5; i++)
     {
         std::stringstream stri;
         stri << i;
-        G4double h = 4.e3/5.*i;
-        G4double density = density0*pow((1-(kappa-1)/kappa*M*g0*h/(R*T)),(1/(kappa-1)));
-        Air[i] = new G4Material("G4_AIR_"+stri.str(), density, 2);
-        Air[i]->AddElement(elN, 70*perCent);
-        Air[i]->AddElement(elO, 30*perCent);
+        G4double h = 4.e3 / 5. * i;
+        //h es la altura de cada capa de la atmosfera (Si se quiere cambiar el largo, ZWorld, de la atmosfera, 4e3 se debe reemplazar por el largo total deseado)
+        //Esta divido entre 5 porque es el numero de capas que queremos si queremos mas o menos capas hay que reemplazar el 5, por el numero de capas deseado.
+        
+        G4double density = density0 * pow((1 - (kappa - 1) / kappa * M * g0 * h / (R * T)), (1 / (kappa - 1)));
+        //density0 * pow((1 - (kappa - 1) / kappa * M * g0 * h / (R * T)), (1 / (kappa - 1)) es la formula del modelo de presion que usamos.
+        
+        Air[i] = new G4Material("G4_AIR_" + stri.str(), density, 3); // G4Material(Nombre del material, densidad, numero de elementos que tiene el material) // ("G4_AIR_" + stri.str() is the name of the material. It uses a string that combines "G4_AIR_" with the value of stri.str() to create unique names for each material.)
+        
+        
+        // La composición del material. Especificamos los elementos y sus proporciones.
+        Air[i]->AddElement(elN, 76 * perCent);
+        Air[i]->AddElement(elO, 23 * perCent);
+        Air[i]->AddElement(elAr, 1 * perCent);
     }
 }
+	//FIN DE LOS MATERIALES
+
+//********************************//
+
+//COnstruccion de la Atmosfera 
 
 void MyDetectorConstruction::ConstructAtmosphere()
 {
     solidAir = new G4Box("solidAir", xWorld, yWorld, zWorld/5.);
 
+//*********************************//
+// Se crean las capas y se les agrega el material Air[0],Air[1]..., Air[4]
+
     for(G4int i = 0; i < 5; i++)
     {
         logicAir[i] = new G4LogicalVolume(solidAir, Air[i], "logicAir");
 
-        physAir[i] = new G4PVPlacement(0, G4ThreeVector(0, 0, zWorld/5.*2.*i - zWorld + zWorld/5.), logicAir[i], "physAir", logicWorld, false, i, true);
+        physAir[i] = new G4PVPlacement(0, G4ThreeVector(0, 0, zWorld/5.*2.*i - zWorld + zWorld/5.), logicAir[i], "physAir", logicWorld, false, i, true); //Aqui tambien se debe reemplazar el 5 que divide zWorld, si se desea cambiar el numero de capas de la atmosfera.
     }
 }
+//FIN DE LA CONSTRUCCION DE LA ATMOSFERA
+
+
+//************************************//
+
+
+// AQUI MEZCLAMOS TODO
 
 G4VPhysicalVolume *MyDetectorConstruction::Construct()
 {
+
+// Se crea el Mundo
+
     solidWorld = new G4Box("solidWorld", xWorld, yWorld, zWorld);
 
     logicWorld = new G4LogicalVolume(solidWorld, worldMat, "logicWorld");
 
     physWorld = new G4PVPlacement(0, G4ThreeVector(0., 0., 0.), logicWorld, "physWorld", 0, false, 0, true);
+    
+// Se anexa la atmosfera
 
     if(isAtmosphere)
         ConstructAtmosphere();
-//////////////////////////////
-    G4Box *solidDetector = new G4Box("solidDetector", 1000*m,1000*m,400*m);
-    logicDetector = new G4LogicalVolume(solidDetector, worldMat, "logicDetector");
+        
+        
+//*****************************//
+// Se Anexa el Detector //
+
+    G4Box *solidDetector = new G4Box("solidDetector", xWorld, yWorld, 400*m); //Dimension del detector, z debe ser la altura de cada capa entre 2, es decir, h/2. No pude hacerlo menos alto, me lanza errores de Overlap, pero si el detector era del largo de la capa entonces funcionaba sin problemas asi que lo deje asi. EL Problema en este caso Es que el suelo no sera nunca 0, el suelo en nuestro caso esta en 800m porque ahi es donde se encuentra la primera capa del detector. Asi que al final las particulas se estan propagando por 3200m de atmosfera y son medidos a una altura de 800msnm.
+    
+    logicDetector = new G4LogicalVolume(solidDetector, worldMat, "logicDetector"); // El material dentro del detector es aire generico. No el que creamos. Pero no importa, igual no afecta en nada las mediciones.
+
+
+    G4VPhysicalVolume *physDetector = new G4PVPlacement(0, G4ThreeVector(0.*m, 0.*m,1600.*m), logicDetector, "physDetector", logicWorld, false, 0, true); // Ubicacion del detector en el mundo. 1600 indica que se esta midiendo en el primer nivel de observacion (cerca del suelo). SI se quiere medir en el segundo nivel z es 800, en el 3er nivel z = 0, en el 4to nivel z= -800, y en el 5to (la que esta pegada a la fuente) z = -1600. Pendiente del signo en Z, es importante. Estos numero son asi porq el largo de la atmosfera es 4000m si se cambia el largo de la atmosfera estos numeros deben cambiar, CUIDADO!!
+
+//*********************************// 
+//FIN DEL DETECTOR
 
-    G4VPhysicalVolume *physDetector = new G4PVPlacement(0, G4ThreeVector(0.*m, 0.*m,800.*m), logicDetector, "physDetector", logicWorld, false, 0, true);
 
 	    return physWorld;
 
 }
 
+//*****************************//
+//Le agregamos sensibilidad al detector
+
 void MyDetectorConstruction::ConstructSDandField()
 {
 	MySensitiveDetector *sensDet = new MySensitiveDetector("SensitiveDetector");
-	logicDetector->SetSensitiveDetector(sensDet);
+	logicDetector->SetSensitiveDetector(sensDet); // LOGICDETECTOR tiene ahora sensibilidad y detecta particulas.
 }