diff --git a/data/CORSIKAValidation.bin b/data/CORSIKAValidation.bin new file mode 100644 index 0000000000..4dde254d3a Binary files /dev/null and b/data/CORSIKAValidation.bin differ diff --git a/macros/validation/CORSIKA.config.mac b/macros/validation/CORSIKA.config.mac new file mode 100644 index 0000000000..8f3d1ae612 --- /dev/null +++ b/macros/validation/CORSIKA.config.mac @@ -0,0 +1,28 @@ +/Generator/CORSIKAGenerator/corsika_file data/CORSIKAValidation.bin +/Generator/CORSIKAGenerator/height 4 m +/Generator/CORSIKAGenerator/plane_x_size 5 m +/Generator/CORSIKAGenerator/plane_y_size 5 m + +/run/verbose 0 +/event/verbose 0 +/tracking/verbose 0 +/process/em/verbose 0 + +/PhysicsList/RegisterPhysics G4EmStandardPhysics_option4 +/PhysicsList/RegisterPhysics G4DecayPhysics +/PhysicsList/RegisterPhysics G4RadioactiveDecayPhysics +/PhysicsList/RegisterPhysics G4HadronElasticPhysicsHP +/PhysicsList/RegisterPhysics G4HadronPhysicsQGSP_BERT_HP +/PhysicsList/RegisterPhysics G4StoppingPhysics +/PhysicsList/RegisterPhysics G4IonPhysics +/PhysicsList/RegisterPhysics NexusPhysics +/PhysicsList/RegisterPhysics G4StepLimiterPhysics + +/Geometry/XeSphere/LXe false +/Geometry/XeSphere/pressure 15. bar +/Geometry/XeSphere/radius 2. m + +### JOB CONTROL +/nexus/random_seed 1 +/nexus/persistency/start_id 0 +/nexus/persistency/outputFile CORSIKA.next diff --git a/macros/validation/CORSIKA.init.mac b/macros/validation/CORSIKA.init.mac new file mode 100644 index 0000000000..adaf3c11cb --- /dev/null +++ b/macros/validation/CORSIKA.init.mac @@ -0,0 +1,15 @@ +# This macro generates cosmic rays taken from a CORSIKA binary +# 4 m above a xenon sphere with radius 2 m, in a plane 5 m x 5 m + +/control/execute macros/physics/DefaultPhysicsList.mac +/nexus/RegisterGeometry XeSphere +/nexus/RegisterGenerator CORSIKAGenerator +/nexus/RegisterPersistencyManager PersistencyManager + +/nexus/RegisterMacro macros/validation/CORSIKA.config.mac + +/control/verbose 0 +/run/verbose 0 +/event/verbose 0 +/tracking/verbose 0 +/process/em/verbose 0 \ No newline at end of file diff --git a/source/generators/CORSIKAGenerator.cc b/source/generators/CORSIKAGenerator.cc new file mode 100644 index 0000000000..3a7031e911 --- /dev/null +++ b/source/generators/CORSIKAGenerator.cc @@ -0,0 +1,278 @@ +// ---------------------------------------------------------------------------- +// nexus | CORSIKAGenerator.cc +// +// This class generates particles taking as input a CORSIKA binary file +// +// The NEXT Collaboration +// ---------------------------------------------------------------------------- + +#include "CORSIKAGenerator.h" +#include "DetectorConstruction.h" +#include "GeometryBase.h" +#include "FactoryBase.h" + +#include +#include +#include +#include +#include +#include +#include + + +using namespace nexus; + +REGISTER_CLASS(CORSIKAGenerator, G4VPrimaryGenerator) + +CORSIKAGenerator::CORSIKAGenerator(): + G4VPrimaryGenerator(), msg_(0), height_(4 * m), plane_x_size_(5 * m), plane_y_size_(5 * m), geom_(0) +{ + msg_ = new G4GenericMessenger(this, "/Generator/CORSIKAGenerator/", + "Control commands of CORSIKAGenerator."); + + msg_->DeclareProperty("corsika_file", corsika_file_, + "Name of the CORSIKA binary file."); + + G4GenericMessenger::Command& height = msg_->DeclareProperty("height", height_, + "Height of the cosmic generation plane."); + height.SetUnitCategory("Length"); + height.SetParameterName("height", false); + height.SetRange("height>0."); + + G4GenericMessenger::Command& plane_x_size = msg_->DeclareProperty("plane_x_size", plane_x_size_, + "Size of the cosmic generation plane on the x axis."); + plane_x_size.SetUnitCategory("Length"); + plane_x_size.SetParameterName("plane_x_size", false); + plane_x_size.SetRange("plane_x_size > 0."); + + G4GenericMessenger::Command& plane_y_size = msg_->DeclareProperty("plane_y_size", plane_y_size_, + "Size of the cosmic generation plane on the y axis."); + plane_y_size.SetUnitCategory("Length"); + plane_y_size.SetParameterName("plane_y_size", false); + plane_y_size.SetRange("plane_y_size > 0."); + + msg_->DeclareProperty("region", region_, + "Set the region of the geometry where the vertex will be generated."); + + DetectorConstruction* detconst = (DetectorConstruction*) G4RunManager::GetRunManager()->GetUserDetectorConstruction(); + geom_ = detconst->GetGeometry(); +} + + +CORSIKAGenerator::~CORSIKAGenerator() +{ + delete msg_; +} + +void CORSIKAGenerator::OpenFile(G4String corsika_binary) +{ + input_ = new std::ifstream(corsika_binary); + + while(input_->read(buf_.ch, 4)) { + unsigned reclen = buf_.in[0]; + // CORSIKA records are in units of 4 bytes + if(reclen % 4) { + throw std::runtime_error("[CORSIKAGenerator] Error: record size not a multiple of 4"); + } + + // We will be looking at at least 8 bytes to determine the + // input file format, and all real CORSIKA records are longer + // than that. + if(reclen < 2*4) { + throw std::runtime_error("[CORSIKAGenerator] Error: reclen too small"); + } + + // Read the full record + if(!input_->read(buf_.ch, reclen)) { + break; + } + + // Determine the format and and store the decision for future blocks. + // We are starting file read, so should see the RUNH marker + // In COMPACT format each block is preceded by 4 bytes + // giving the size of the block in words. + + if(!strncmp(buf_.ch+0, "RUNH", 4)) { + G4cout << "[CORSIKAGenerator] Reading NORMAL format" << G4endl; + infmt_ = Format::NORMAL; + } + else if(!strncmp(buf_.ch+4, "RUNH", 4)) { + G4cout << "[CORSIKAGenerator] Reading COMPACT format" << G4endl; + infmt_ = Format::COMPACT; + } + else { + throw std::runtime_error("[CORSIKAGenerator] Error: did not find the RUNH record to determine COMPACT flag"); + } + + unsigned iword = 0; + if(infmt_ == Format::COMPACT) { + // Move to the beginning of the actual block + ++iword; + } + + break; + + } + input_->clear(); + input_->seekg(0, std::ios::beg); +} + +double CORSIKAGenerator::WrapVarBoxNo(const float var, const float low, const float high, int &boxno) +{ + //wrap variable so that it's always between low and high + boxno = int(floor(var / (high - low))); + return (var - (high - low) * floor(var / (high - low))) + low; +} + +void CORSIKAGenerator::GeneratePrimaryVertex(G4Event* event) +{ + // FORTRAN sequential records are prefixed with their length + // in a 4-byte word + if (!is_open_) { + OpenFile(corsika_file_); + is_open_ = true; + } + + unsigned n_part = 0; + std::map, std::vector< std::pair>> particles_map; + + while(input_->read(buf_.ch, 4)) { + + unsigned reclen = buf_.in[0]; + // CORSIKA records are in units of 4 bytes + if(reclen % 4) { + throw std::runtime_error("Error: record size not a multiple of 4"); + } + + // We will be looking at at least 8 bytes to determine the + // input file format, and all real CORSIKA records are longer + // than that. + if(reclen < 2*4) { + throw std::runtime_error("Error: reclen too small"); + } + + if(reclen > 4*fbsize_words_) { + throw std::runtime_error("Error: reclen too big"); + } + + // Read the full record + if(!input_->read(buf_.ch, reclen)) { + break; + } + + //================================================================ + // Go over blocks in the record + for(unsigned iword = 0; iword < reclen/4; ) { + + unsigned block_words = (infmt_ == Format::COMPACT) ? buf_.in[iword] : 273; + + if(!block_words) { + throw std::runtime_error("Got block_words = 0\n"); + } + + if(infmt_ == Format::COMPACT) { + // Move to the beginning of the actual block + ++iword; + } + + std::string event_marker = (infmt_ == Format::NORMAL || !event_count_) ? "EVTH" : "EVHW"; + + // Determine the type of the data block + if (!strncmp(buf_.ch+4*iword, "RUNH", 4)) { + run_number_ = lrint(buf_.fl[1+iword]); + } + else if(!strncmp(buf_.ch+4*iword, "RUNE", 4)) { + unsigned end_run_number = lrint(buf_.fl[1+iword]); + unsigned end_event_count = lrint(buf_.fl[2+iword]); + if(end_run_number != run_number_) { + throw std::runtime_error("Error: run number mismatch in end of run record\n"); + } + if(event_count_ != end_event_count) { + std::cerr<<"RUNE: _event_count = " << event_count_ << " end record = " << end_event_count << G4endl; + throw std::runtime_error("Error: event count mismatch in end of run record\n"); + } + // Exit the read loop at the end of run + // _primaries = 0; + G4cout << "[CORSIKAGenerator] End of CORSIKA file" << G4endl; + break; + } + else if(!strncmp(buf_.ch+4*iword, event_marker.data(), 4)) { + ++event_count_; + current_event_number_ = lrint(buf_.fl[1+iword]); + G4cout << "[CORSIKAGenerator] Event " << current_event_number_ << G4endl; + particles_map.erase(particles_map.begin(), particles_map.end()); + // ++_primaries; + } else if(!strncmp(buf_.ch+4*iword, "EVTE", 4)) { + unsigned end_event_number = lrint(buf_.fl[1+iword]); + if (end_event_number != current_event_number_) { + throw std::runtime_error("Error: event number mismatch in end of event record\n"); + } + } else { + for (unsigned i_part = 0; i_part < block_words; i_part+=7) { + unsigned id = buf_.fl[iword + i_part] / 1000; + if (id == 0){ + continue; + } + n_part++; + const int pdgId = corsikaToPdgId.at(id); + + // The footprint of a CORSIKA shower is very large. We tile the event + // with user-defined dimensions, then we pick only the first tile, discarding the rest. + int box_x = 0; + const float x = WrapVarBoxNo(buf_.fl[iword + i_part + 4] * cm, -plane_x_size_/2, plane_x_size_/2, box_x); + int box_y = 0; + const float y = WrapVarBoxNo(buf_.fl[iword + i_part + 5] * cm, -plane_y_size_/2, plane_y_size_/2, box_y); + std::pair xy(box_x, box_y); + + const float px = buf_.fl[iword + i_part + 1] * GeV; + const float py = buf_.fl[iword + i_part + 2] * GeV; + const float pz = -buf_.fl[iword + i_part + 3] * GeV; // In CORSIKA the momentum is in the -z direction + const float t = buf_.fl[iword + i_part + 6]; + G4PrimaryVertex* vertex = new G4PrimaryVertex(G4ThreeVector(x,y,height_), t); + + G4ParticleDefinition* particle_definition_ = + G4ParticleTable::GetParticleTable()->FindParticle(pdgId); + G4PrimaryParticle* particle = new G4PrimaryParticle(particle_definition_, px, py, pz); + + // Here we store the particles in a map where the key is the tile ID and + // the value is a pair of vertex, particle + std::pair vp(vertex, particle); + if (particles_map.count(xy) == 0) { + std::vector> vps; + vps.push_back(vp); + particles_map.insert({xy, vps}); + } else { + particles_map[xy].push_back(vp); + } + + } + + } + + // Move to the next block + iword += block_words; + + } // loop over blocks in a record + + // Here we expect the FORTRAN end of record padding, + // read and verify its value. + if(!input_->read(buf_.ch, 4)) { + break; + } + if(buf_.in[0] != reclen) { + throw std::runtime_error("Error: unexpected FORTRAN record end padding"); + } + + if (n_part > 0) { + // Pick the first tile and insert the particle into the event + for (auto vp: particles_map.begin()->second) { + vp.first->SetPrimary(vp.second); + event->AddPrimaryVertex(vp.first); + } + G4cout << "[CORSIKAGenerator] End event " << event_count_ << " with " << n_part << " particles" << G4endl; + n_part = 0; + break; + } + + } // loop over records +} diff --git a/source/generators/CORSIKAGenerator.h b/source/generators/CORSIKAGenerator.h new file mode 100644 index 0000000000..86463aaabd --- /dev/null +++ b/source/generators/CORSIKAGenerator.h @@ -0,0 +1,205 @@ +// ---------------------------------------------------------------------------- +// nexus | CORSIKAGenerator.h +// +// This class generates particles taking as input a CORSIKA binary file +// +// The NEXT Collaboration +// ---------------------------------------------------------------------------- + +#ifndef CORSIKA_GENERATOR_H +#define CORSIKA_GENERATOR_H + +#include +#include + +class G4GenericMessenger; +class G4Event; +class G4ParticleDefinition; + + +namespace nexus { + + enum class Format { UNDEFINED, NORMAL, COMPACT }; + const std::map corsikaToPdgId = { + {1, 22}, // gamma + {2, -11}, // e+ + {3, 11}, // e- + {5, -13}, // mu+ + {6, 13}, // mu- + {7, 111}, // pi0 + {8, 211}, // pi+ + {9, -211}, // pi- + {10, 130}, // K0_L + {11, 321}, // K+ + {12, -321}, // K- + {13, 2112}, // n + {14, 2212}, // p + {15, -2212}, // pbar + {16, 310}, // K0_S + {17, 221}, // eta + {18, 3122}, // Lambda + {19, 3222}, // Sigma+ + {20, 3212}, // Sigma0 + {21, 3112}, // Sigma- + {22, 3322}, // Cascade0 + {23, 3312}, // Cascade- + {24, 3334}, // Omega- + {25, -2112}, // nbar + {26, -3122}, // Lambdabar + {27, -3112}, // Sigma-bar + {28, -3212}, // Sigma0bar + {29, -3222}, // Sigma+bar + {30, -3322}, // Cascade0bar + {31, -3312}, // Cascade+bar + {32, -3334}, // Omega+bar + + {50, 223}, // omega + {51, 113}, // rho0 + {52, 213}, // rho+ + {53, -213}, // rho- + {54, 2224}, // Delta++ + {55, 2214}, // Delta+ + {56, 2114}, // Delta0 + {57, 1114}, // Delta- + {58, -2224}, // Delta--bar + {59, -2214}, // Delta-bar + {60, -2114}, // Delta0bar + {61, -1114}, // Delta+bar + {62, 10311}, // K*0 + {63, 10321}, // K*+ + {64, -10321}, // K*- + {65, -10311}, // K*0bar + {66, 12}, // nu_e + {67, -12}, // nu_ebar + {68, 14}, // nu_mu + {69, -14}, // nu_mubar + + {116, 421}, // D0 + {117, 411}, // D+ + {118, -411}, // D-bar + {119, -421}, // D0bar + {120, 431}, // D+_s + {121, -431}, // D-_sbar + {122, 441}, // eta_c + {123, 423}, // D*0 + {124, 413}, // D*+ + {125, -413}, // D*-bar + {126, -423}, // D*0bar + {127, 433}, // D*+_s + {128, -433}, // D*-_s + + {130, 443}, // J/Psi + {131, -15}, // tau+ + {132, 15}, // tau- + {133, 16}, // nu_tau + {134, -16}, // nu_taubar + + {137, 4122}, // Lambda+_c + {138, 4232}, // Cascade+_c + {139, 4132}, // Cascade0_c + {140, 4222}, // Sigma++_c + {141, 4212}, // Sigma+_c + {142, 4112}, // Sigma0_c + {143, 4322}, // Cascade'+_c + {144, 4312}, // Cascade'0_c + {145, 4332}, // Omega0_c + {149, -4122}, // Lambda-_cbar + {150, -4232}, // Cascade-_cbar + {151, -4132}, // Cascade0_cbar + {152, -4222}, // Sigma--_cbar + {153, -4212}, // Sigma-_cbar + {154, -4112}, // Sigma0_cbar + {155, -4322}, // Cascade'-_cbar + {156, -4312}, // Cascade'0_cbar + {157, -4332}, // Omega0_cbar + {161, 4224}, // Sigma*++_c + {162, 1214}, // Sigma*+_c + {163, 4114}, // Sigma*0_c + + {171, -4224}, // Sigma*--_cbar + {172, -1214}, // Sigma*-_cbar + {173, -4114}, // Sigma*0_cbar + {176, 511}, // B0 + {177, 521}, // B+ + {178, -521}, // B-bar + {179, -511}, // B0bar + {180, 531}, // B0_s + {181, -531}, // B0_sbar + {182, 541}, // B+_c + {183, -541}, // B-_cbar + {184, 5122}, // Lambda0_b + {185, 5112}, // Sigma-_b + {186, 5222}, // Sigma+_b + {187, 5232}, // Cascade0_b + {188, 5132}, // Cascade-_b + {189, 5332}, // Omega-_b + {190, -5112}, // Lambda0_bbar + {191, -5222}, // Sigma+_bbar + {192, -5112}, // Sigma-_bbar + {193, -5232}, // Cascade0_bbar + {194, -5132}, // Cascade+_bbar + {195, -5332}, // Omega+_bbar + + {201, 1000010020}, // Deuteron + {301, 1000010030}, // Tritium + {402, 1000020040}, // alpha + {5626, 1000260560}, // Iron + {1206, 1000080120}, // Carbon + {1407, 1000070140}, // Nitrogen + {1608, 1000080160}, // Oxygen + {2713, 1000130270}, // Aluminum + {3216, 1000160320}, // Sulfur + {2814, 1000140280}, // Silicon + {9900, 22} // Cherenkov gamma + }; + + class GeometryBase; + + class CORSIKAGenerator: public G4VPrimaryGenerator + { + public: + /// Constructor + CORSIKAGenerator(); + /// Destructor + ~CORSIKAGenerator(); + + /// This method is invoked at the beginning of the event. It sets + /// a primary vertex (that is, a particle in a given position and time) + /// in the event. + void GeneratePrimaryVertex(G4Event*); + void OpenFile(G4String); + double WrapVarBoxNo(const float var, const float low, const float high, int &boxno); + + private: + + private: + static constexpr unsigned fbsize_words_ = 5733 + 2; + + union { + float fl[fbsize_words_]; + unsigned in[fbsize_words_]; + char ch[sizeof(fl[0])*fbsize_words_]; + } buf_; + + Format infmt_ = Format::UNDEFINED; + G4GenericMessenger* msg_; + + G4String region_; + G4String corsika_file_; + G4double height_; + G4double plane_x_size_; + G4double plane_y_size_; + G4bool is_open_ = false; + unsigned run_number_ = 0; + unsigned current_event_number_ = 0; + unsigned event_count_ = 0; + + std::ifstream *input_; + + const GeometryBase* geom_; ///< Pointer to the detector geometry + + }; + +} // end namespace nexus + +#endif