diff --git a/projects/dataclasses/public/SIREN/dataclasses/ParticleTypes.def b/projects/dataclasses/public/SIREN/dataclasses/ParticleTypes.def index 73b578f9..af36148c 100644 --- a/projects/dataclasses/public/SIREN/dataclasses/ParticleTypes.def +++ b/projects/dataclasses/public/SIREN/dataclasses/ParticleTypes.def @@ -194,6 +194,7 @@ X(NuclInt, -2000001004) X(MuPair, -2000001005) X(Hadrons, -2000001006) X(Decay, -2000001007) +X(Hadronization, 99) X(ContinuousEnergyLoss, -2000001111) X(FiberLaser, -2000002100) X(N2Laser, -2000002101) diff --git a/projects/detector/private/DetectorModel.cxx b/projects/detector/private/DetectorModel.cxx index 75a86dc7..37cce521 100644 --- a/projects/detector/private/DetectorModel.cxx +++ b/projects/detector/private/DetectorModel.cxx @@ -683,7 +683,7 @@ double DetectorModel::GetInteractionDensity(Geometry::IntersectionList const & i std::vector const & total_cross_sections, double const & total_decay_length) const { Vector3D direction = p0 - intersections.position; - if(direction.magnitude() == 0 || direction.magnitude() <= 1e-6) { + if(direction.magnitude() == 0 || direction.magnitude() <= 1e-5) { direction = intersections.direction; } else { direction.normalize(); @@ -987,8 +987,9 @@ double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const if(distance == 0.0) { return 0.0; } - if(direction.magnitude() <= 1e-6) { + if(direction.magnitude() <= 1e-5) { direction = intersections.direction; + // return 1e-05; // I have to do this to ensure that it works } direction.normalize(); diff --git a/projects/distributions/private/secondary/vertex/SecondaryBoundedVertexDistribution.cxx b/projects/distributions/private/secondary/vertex/SecondaryBoundedVertexDistribution.cxx index a86fffb1..43705a3e 100644 --- a/projects/distributions/private/secondary/vertex/SecondaryBoundedVertexDistribution.cxx +++ b/projects/distributions/private/secondary/vertex/SecondaryBoundedVertexDistribution.cxx @@ -53,8 +53,6 @@ double log_one_minus_exp_of_negative(double x) { void SecondaryBoundedVertexDistribution::SampleVertex(std::shared_ptr rand, std::shared_ptr detector_model, std::shared_ptr interactions, siren::dataclasses::SecondaryDistributionRecord & record) const { - // std::cout << "in sample bounded vertex" << std::endl; - siren::math::Vector3D pos = record.initial_position; siren::math::Vector3D dir = record.direction; @@ -123,8 +121,6 @@ void SecondaryBoundedVertexDistribution::SampleVertex(std::shared_ptr detector_model, std::shared_ptr interactions, siren::dataclasses::InteractionRecord const & record) const { - // std::cout << "in sample bounded vertex gen prob" << std::endl; - siren::math::Vector3D dir(record.primary_momentum[1], record.primary_momentum[2], record.primary_momentum[3]); dir.normalize(); siren::math::Vector3D vertex(record.interaction_vertex); diff --git a/projects/distributions/private/secondary/vertex/SecondaryPhysicalVertexDistribution.cxx b/projects/distributions/private/secondary/vertex/SecondaryPhysicalVertexDistribution.cxx index d49a42f0..f5d327b7 100644 --- a/projects/distributions/private/secondary/vertex/SecondaryPhysicalVertexDistribution.cxx +++ b/projects/distributions/private/secondary/vertex/SecondaryPhysicalVertexDistribution.cxx @@ -51,22 +51,15 @@ double log_one_minus_exp_of_negative(double x) { void SecondaryPhysicalVertexDistribution::SampleVertex(std::shared_ptr rand, std::shared_ptr detector_model, std::shared_ptr interactions, siren::dataclasses::SecondaryDistributionRecord & record) const { - // // std::cout << "in sample physical vertex" << std::endl; siren::math::Vector3D pos = record.initial_position; siren::math::Vector3D dir = record.direction; - // // std::cout << "in sample physical vertex-1" << std::endl; siren::math::Vector3D endcap_0 = pos; // treat hadronizations differntely if (interactions->HasHadronizations()) { - // std::cout << "in sample physical vertex-hadron" << std::endl; - record.SetLength(0); return; } - // // std::cout << "in sample physical vertex-shouldnm't be here" << std::endl; - - siren::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), std::numeric_limits::infinity()); path.ClipToOuterBounds(); @@ -86,7 +79,6 @@ void SecondaryPhysicalVertexDistribution::SampleVertex(std::shared_ptr detector_model, std::shared_ptr interactions, siren::dataclasses::InteractionRecord const & record) const { - // std::cout << "in sample physical vertex gen prob" << std::endl; - siren::math::Vector3D dir(record.primary_momentum[1], record.primary_momentum[2], record.primary_momentum[3]); dir.normalize(); siren::math::Vector3D vertex(record.interaction_vertex); diff --git a/projects/distributions/private/secondary/vertex/SecondaryVertexPositionDistribution.cxx b/projects/distributions/private/secondary/vertex/SecondaryVertexPositionDistribution.cxx index 70af0664..87b44e63 100644 --- a/projects/distributions/private/secondary/vertex/SecondaryVertexPositionDistribution.cxx +++ b/projects/distributions/private/secondary/vertex/SecondaryVertexPositionDistribution.cxx @@ -15,7 +15,6 @@ namespace distributions { // class SecondaryVertexPositionDistribution : InjectionDistribution //--------------- void SecondaryVertexPositionDistribution::Sample(std::shared_ptr rand, std::shared_ptr detector_model, std::shared_ptr interactions, siren::dataclasses::SecondaryDistributionRecord & record) const { - // // // // // std::cout << "sampling vertex" << std::endl; SampleVertex(rand, detector_model, interactions, record); } diff --git a/projects/injection/private/Injector.cxx b/projects/injection/private/Injector.cxx index 9bf04c1b..a77338cf 100644 --- a/projects/injection/private/Injector.cxx +++ b/projects/injection/private/Injector.cxx @@ -188,7 +188,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record double fake_prob; // if contains hadronization, then perform only hadronization if (interactions->HasHadronizations()) { - // std::cout << "saw hadronization" << std::endl; double total_frag_prob = 0; std::vector frag_probs; for(auto const & hadronization : interactions->GetHadronizations() ) { @@ -209,9 +208,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record } } - // std::cout << "Hadronization finished signatures" << std::endl; - - // now choose the specific charmed hadron to fragment into double r = random->Uniform(0, total_frag_prob); unsigned int index = 0; @@ -224,34 +220,22 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record matching_hadronizations[index]->SampleFinalState(xsec_record, random); xsec_record.Finalize(record); - - // std::cout << "hadronization done!" << std::endl; - } else { if (interactions->HasCrossSections()) { - // std::cout << "saw xsec" << std::endl; for(auto const target : available_targets) { if(possible_targets.find(target) != possible_targets.end()) { - // std::cout << "saw xsec: in first for loop" << std::endl; // Get target density double target_density = detector_model->GetParticleDensity(intersections, DetectorPosition(interaction_vertex), target); // Loop over cross sections that have this target std::vector> const & target_cross_sections = interactions->GetCrossSectionsForTarget(target); for(auto const & cross_section : target_cross_sections) { - // std::cout << "saw xsec: in second for loop" << std::endl; // Loop over cross section signatures with the same target std::vector signatures = cross_section->GetPossibleSignaturesFromParents(record.signature.primary_type, target); for(auto const & signature : signatures) { - // std::cout << "saw xsec: in third for loop" << std::endl; - fake_record.signature = signature; fake_record.target_mass = detector_model->GetTargetMass(target); // Add total cross section times density to the total prob - // std::cout << "about to sample total xsec" << std::endl; - fake_prob = target_density * cross_section->TotalCrossSection(fake_record); - // std::cout << "finished sampling total xsec" << std::endl; - total_prob += fake_prob; xsec_prob += fake_prob; // Add total prob to probs @@ -266,7 +250,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record } } if (interactions->HasDecays()) { - // std::cout << "saw decay" << std::endl; for(auto const & decay : interactions->GetDecays() ) { for(auto const & signature : decay->GetPossibleSignaturesFromParent(record.signature.primary_type)) { fake_record.signature = signature; @@ -282,7 +265,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record } } } - // std::cout << "continuing...." << std::endl; if(total_prob == 0) throw(siren::utilities::InjectionFailure("No valid interactions for this event!")); // Throw a random number @@ -293,7 +275,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record record.signature.target_type = matching_targets[index]; record.signature = matching_signatures[index]; double selected_prob = 0.0; - // std::cout << "finished initializing stuff" << std::endl; for(unsigned int i=0; i 0 ? probs[i] - probs[i - 1] : probs[i]); @@ -303,14 +284,9 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record throw(siren::utilities::InjectionFailure("No valid interactions for this event!")); record.target_mass = detector_model->GetTargetMass(record.signature.target_type); siren::dataclasses::CrossSectionDistributionRecord xsec_record(record); - // std::cout << "finished sampling from primary process" << std::endl; if(r <= xsec_prob) { - // std::cout << "about to sample primary process final state" << std::endl; - matching_cross_sections[index]->SampleFinalState(xsec_record, random); } else { - // std::cout << "about to sample primary process final state" << std::endl; - matching_decays[index - matching_cross_sections.size()]->SampleFinalState(xsec_record, random); } xsec_record.Finalize(record); @@ -321,8 +297,6 @@ void Injector::SampleCrossSection(siren::dataclasses::InteractionRecord & record // // Modifies an InteractionRecord with the new event siren::dataclasses::InteractionRecord Injector::SampleSecondaryProcess(siren::dataclasses::SecondaryDistributionRecord & secondary_record) const { - // std::cout << "sampling secondary" << std::endl; - // std::cout << "secondary record type is " << secondary_record.type << " " << secondary_record.id << std::endl; std::shared_ptr secondary_process = secondary_process_map.at(secondary_record.type); std::shared_ptr secondary_interactions = secondary_process->GetInteractions(); std::vector> secondary_distributions = secondary_process->GetSecondaryInjectionDistributions(); @@ -332,32 +306,15 @@ siren::dataclasses::InteractionRecord Injector::SampleSecondaryProcess(siren::da size_t failed_tries = 0; while(true) { tries += 1; - // // std::cout << "gotcha" << std::endl; - // for(auto & distribution : secondary_distributions) { - // // std::cout << "sample distribution" << std::endl; - // distribution->Sample(random, detector_model, secondary_process->GetInteractions(), secondary_record); - // } - // siren::dataclasses::InteractionRecord record; - // secondary_record.Finalize(record); - // // // std::cout << "sample distribution" << std::endl; - - // SampleCrossSection(record, secondary_interactions); - // return record; - try { for(auto & distribution : secondary_distributions) { - // std::cout << "sample distribution" << std::endl; distribution->Sample(random, detector_model, secondary_process->GetInteractions(), secondary_record); } siren::dataclasses::InteractionRecord record; secondary_record.Finalize(record); - // // std::cout << "sample distribution" << std::endl; - SampleCrossSection(record, secondary_interactions); return record; } catch(siren::utilities::InjectionFailure const & e) { - // std::cout << "caught error" << std::endl; - failed_tries += 1; if(failed_tries > max_tries) { throw(siren::utilities::InjectionFailure("Failed to generate secondary process!")); @@ -379,7 +336,6 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() { size_t tries = 0; size_t failed_tries = 0; // Initial Process - // std::cout << "Sampling primary interactions" << std::endl; while(true) { tries += 1; try { @@ -388,7 +344,6 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() { distribution->Sample(random, detector_model, primary_process->GetInteractions(), primary_record); } primary_record.Finalize(record); - // std::cout << "primary record fixed" << std:: endl; SampleCrossSection(record); break; } catch(siren::utilities::InjectionFailure const & e) { @@ -408,12 +363,9 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() { std::shared_ptr parent = tree.add_entry(record); // Secondary Processes - // std::cout << "Sampling primary interactions 2" << std::endl; - std::deque, std::shared_ptr>> secondaries; std::function)> add_secondaries = [&](std::shared_ptr parent) { for(size_t i=0; irecord.signature.secondary_types.size(); ++i) { - // // std::cout << "for loop 1" << std::endl; siren::dataclasses::ParticleType const & type = parent->record.signature.secondary_types[i]; std::map>::iterator it = secondary_process_map.find(type); if(it == secondary_process_map.end()) { @@ -430,70 +382,18 @@ siren::dataclasses::InteractionTree Injector::GenerateEvent() { }; add_secondaries(parent); - // std::cout << "num secondaries: " << secondaries.size() << std::endl; while(secondaries.size() > 0) { - // // std::cout << "while loop 1" << std::endl; - for(int i = secondaries.size() - 1; i >= 0; --i) { - // // std::cout << "for loop 2" << std::endl; - std::shared_ptr parent = std::get<0>(secondaries[i]); std::shared_ptr secondary_dist = std::get<1>(secondaries[i]); - // // std::cout << "for loop 2-1" << std::endl; secondaries.erase(secondaries.begin() + i); - - // std::cout << "Primary Type: " << secondary_dist->record.signature.primary_type << std::endl; - // std::cout << "Secondary Types: "; - // for (const auto& type : secondary_dist->record.signature.secondary_types) { - // std::cout << type << " "; - // } - // std::cout << std::endl; - - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - // this try-catch clock is to debug the secondary process LE no available interaction error - try{ - siren::dataclasses::InteractionRecord secondary_record = SampleSecondaryProcess(*secondary_dist); - std::shared_ptr secondary_datum = tree.add_entry(secondary_record, parent); - add_secondaries(secondary_datum); - } catch (const std::exception& e) { - std::cerr << "Error occurred: " << e.what() << std::endl; - - // Print the primary type and secondary types for debugging - std::cerr << "Primary Type: " << secondary_dist->record.signature.primary_type << std::endl; - std::cerr << "Secondary Types: "; - for (const auto& type : secondary_dist->record.signature.secondary_types) { - std::cerr << type << " "; - } - std::cerr << std::endl; - - // Print the primary momentum - std::cerr << "Primary Momentum: "; - for (double component : secondary_dist->record.primary_momentum) { - std::cerr << component << " "; - } - std::cerr << std::endl; - - // Print the secondary IDs - std::cerr << "Secondary IDs: "; - for (const auto& id : secondary_dist->record.secondary_ids) { - std::cerr << id << " "; - } - std::cerr << std::endl; - throw; - } catch (...) { - std::cerr << "Unknown exception caught!" << std::endl; - throw; - } - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - + siren::dataclasses::InteractionRecord secondary_record = SampleSecondaryProcess(*secondary_dist); + std::shared_ptr secondary_datum = tree.add_entry(secondary_record, parent); + add_secondaries(secondary_datum); + } - // // std::cout << "while loop 1-2" << std::endl; - } injected_events += 1; return tree; diff --git a/projects/injection/private/Weighter.cxx b/projects/injection/private/Weighter.cxx index afd11cfb..bd166909 100644 --- a/projects/injection/private/Weighter.cxx +++ b/projects/injection/private/Weighter.cxx @@ -114,7 +114,7 @@ double Weighter::EventWeight(siren::dataclasses::InteractionTree const & tree) c double generation_probability = injectors[idx]->EventsToInject();//GenerationProbability(tree); for(auto const & datum : tree.tree) { // skip weighting if record contains hadronization - if (isQuark(datum->record.signature.primary_type)) { + if (datum->record.signature.target_type == siren::dataclasses::Particle::ParticleType::Hadronization) { continue; } std::tuple bounds; diff --git a/projects/interactions/private/CharmDISFromSpline.cxx b/projects/interactions/private/CharmDISFromSpline.cxx index d2474056..42189f5a 100644 --- a/projects/interactions/private/CharmDISFromSpline.cxx +++ b/projects/interactions/private/CharmDISFromSpline.cxx @@ -50,16 +50,10 @@ bool kinematicallyAllowed(double x, double y, double E, double M, double m) { //the numerator of b (or b*d) double bd = sqrt(term * term - ((m * m) / (E * E))); - // also try the D-Meson Mass? double s = 2 * M * E; double Q2 = s * x * y; double Mc = siren::utilities::Constants::D0Mass; - // if ((Q2 / (1 - x) + pow(M, 2) < pow(M + Mc, 2))) { - // std::cout << "SIREN D Meson constraint is trigged!" << std::endl; - // } return ((ad - bd) <= d * y and d * y <= (ad + bd)) && (Q2 / (1 - x) + pow(M, 2) >= pow(M + Mc, 2)); //Eq. 7 - // return ((ad - bd) <= d * y and d * y <= (ad + bd)); //Eq. 7 - } } @@ -197,8 +191,6 @@ void CharmDISFromSpline::ReadParamsFromSplineTable() { // returns true if successfully read minimum Q2 bool q2_good = differential_cross_section_.read_key("Q2MIN", minimum_Q2_); - // std::cout << "reading results: " << mass_good << " " << int_good << " " << q2_good << std::endl; - if(!int_good) { // assume DIS to preserve compatability with previous versions interaction_type_ = 1; @@ -221,9 +213,7 @@ void CharmDISFromSpline::ReadParamsFromSplineTable() { } } else { - // // std::cout << "mass and int both not good" << std::endl; if(differential_cross_section_.get_ndim() == 3) { - // std::cout << "dim is 3" << std::endl; target_mass_ = siren::utilities::Constants::isoscalarMass; } else if(differential_cross_section_.get_ndim() == 2) { @@ -233,8 +223,6 @@ void CharmDISFromSpline::ReadParamsFromSplineTable() { } } } - - // std::cout << "target mass is " << target_mass_ << std::endl; } void CharmDISFromSpline::InitializeSignatures() { @@ -293,7 +281,6 @@ double CharmDISFromSpline::TotalCrossSection(dataclasses::InteractionRecord cons rk::P4 p1(geom3::Vector3(interaction.primary_momentum[1], interaction.primary_momentum[2], interaction.primary_momentum[3]), interaction.primary_mass); double primary_energy; primary_energy = interaction.primary_momentum[0]; - // std::cout << "primary energy is " << primary_energy << std::endl; // if we are below threshold, return 0 if(primary_energy < InteractionThreshold(interaction)) return 0; @@ -304,7 +291,6 @@ double CharmDISFromSpline::TotalCrossSection(siren::dataclasses::ParticleType pr if(not primary_types_.count(primary_type)) { throw std::runtime_error("Supplied primary not supported by cross section!"); } - // std::cout << "now in real sample total xsec func" << std::endl; double log_energy = log10(primary_energy); if(log_energy < total_cross_section_.lower_extent(0) @@ -316,9 +302,7 @@ double CharmDISFromSpline::TotalCrossSection(siren::dataclasses::ParticleType pr } int center; - // std::cout << "maybe problem is here?" << std::endl; total_cross_section_.searchcenters(&log_energy, ¢er); - // std::cout << "maybe problem is here??" << std::endl; double log_xs = total_cross_section_.ndsplineeval(&log_energy, ¢er, 0); @@ -355,7 +339,6 @@ double CharmDISFromSpline::DifferentialCrossSection(dataclasses::InteractionReco } double CharmDISFromSpline::DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass, double Q2) const { - bool check_criteria = false; // this is used to gauge kinematic constraint in xsec and SIREN impementations, will eventually be deleted double log_energy = log10(energy); // check preconditions if(log_energy < differential_cross_section_.lower_extent(0) @@ -375,12 +358,9 @@ double CharmDISFromSpline::DifferentialCrossSection(double energy, double x, dou if(Q2 < minimum_Q2_) // cross section not calculated, assumed to be zero return 0; - if (!check_criteria) { - if(!kinematicallyAllowed(x, y, energy, target_mass_, secondary_lepton_mass)) { - return 0; - } + if(!kinematicallyAllowed(x, y, energy, target_mass_, secondary_lepton_mass)) { + return 0; } - std::array coordinates{{log_energy, log10(x), log10(y)}}; std::array centers; if(!differential_cross_section_.searchcenters(coordinates.data(), centers.data())) @@ -388,30 +368,6 @@ double CharmDISFromSpline::DifferentialCrossSection(double energy, double x, dou double result = pow(10., differential_cross_section_.ndsplineeval(coordinates.data(), centers.data(), 0)); assert(result >= 0); - if (check_criteria) { - // this is a check of kinematic constraint implementation - if (result == 0) { - if(kinematicallyAllowed(x, y, energy, target_mass_, secondary_lepton_mass)) { - std::cout << "xsec gives 0 but kinematically allowed!" << std::endl; - } - } - - if(!kinematicallyAllowed(x, y, energy, target_mass_, secondary_lepton_mass)) { - // check if this is due to charm production constraint - double M = target_mass_; - double E = energy; - double s = 2 * M * E; - double Q2 = s * x * y; - double Mc = siren::utilities::Constants::D0Mass; - if ((Q2 / (1 - x) + pow(M, 2) < pow(M + Mc, 2))) { // if so check result - if (result != 0) { - std::cout << "SIREN constraint not passed but xsec does not give 0!" << std::endl; - } - } - return 0; - } - } - return unit * result; @@ -428,7 +384,6 @@ void CharmDISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionR if (differential_cross_section_.get_ndim() != 3) { throw std::runtime_error("I expected 3 dimensions in the cross section spline, but got " + std::to_string(differential_cross_section_.get_ndim()) +". Maybe your fits file doesn't have the right 'INTERACTION' key?"); } - // std::cout << "in sample final state of charm DIS from spline" << std::endl; rk::P4 p1(geom3::Vector3(record.primary_momentum[1], record.primary_momentum[2], record.primary_momentum[3]), record.primary_mass); rk::P4 p2(geom3::Vector3(0, 0, 0), record.target_mass); @@ -450,11 +405,9 @@ void CharmDISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionR double m = GetLeptonMass(record.signature.secondary_types[lepton_index]); double m1 = record.primary_mass; - // std::cout << "getting mass of primary: " << m1 << std::endl; double m3 = m; double E1_lab = p1_lab.e(); double E2_lab = p2_lab.e(); - // std::cout << "getting energy of primary: " << E1_lab << std::endl; // The out-going particle always gets at least enough energy for its rest mass double yMax = 1 - m / primary_energy; @@ -494,16 +447,12 @@ void CharmDISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionR // rejection sample a point which is kinematically allowed by calculation limits double trialQ; double trials = 0; - // std::cout << "do loop 1" << std::endl; do { - // std::cout << "do loop 2" << std::endl; if (trials >= 100) throw std::runtime_error("too much trials"); trials += 1; kin_vars[1] = random->Uniform(logXMin,0); kin_vars[2] = random->Uniform(logYMin,logYMax); trialQ = (2 * E1_lab * E2_lab) * pow(10., kin_vars[1] + kin_vars[2]); - // std::cout << kin_vars[1] << " " << kin_vars[2] << " " << trialQ << " " << minimum_Q2_ << std::endl; - // std::cout << primary_energy << " " << target_mass_ << " " << m << std::endl; } while(trialQ CharmHadronization::GetPossibleSi std::vector signatures; dataclasses::InteractionSignature signature; signature.primary_type = primary; - signature.target_type = siren::dataclasses::Particle::ParticleType::Decay; + signature.target_type = siren::dataclasses::Particle::ParticleType::Hadronization; signature.secondary_types.resize(2); signature.secondary_types[0] = siren::dataclasses::Particle::ParticleType::Hadrons; diff --git a/projects/interactions/private/CharmMesonDecay.cxx b/projects/interactions/private/CharmMesonDecay.cxx index 7f6520ab..5b92f0f1 100644 --- a/projects/interactions/private/CharmMesonDecay.cxx +++ b/projects/interactions/private/CharmMesonDecay.cxx @@ -148,11 +148,9 @@ double CharmMesonDecay::TotalDecayWidthForFinalState(dataclasses::InteractionRec siren::dataclasses::Particle::ParticleType::EPlus, siren::dataclasses::Particle::ParticleType::NuE}; if (primary == siren::dataclasses::Particle::ParticleType::DPlus && secondaries == k0_eplus_nue) { - // branching_ratio = 0.089; branching_ratio = 1; tau = 1040 * (1e-15); } else if (primary == siren::dataclasses::Particle::ParticleType::D0 && secondaries == kminus_eplus_nue) { - // branching_ratio = 0.03538; branching_ratio = 1; tau = 410.1 * (1e-15); } @@ -296,34 +294,6 @@ void CharmMesonDecay::computeDiffGammaCDF(std::vector constants, double cdf_vector.push_back(1); pdf_vector.push_back(0); - // for debugging and plotting, print the pdf and cdf tables - // for (size_t i = 0; auto& element : cdf_Q2_nodes) { - // std::cout << element; - // // Print comma if it's not the last element - // if (++i != cdf_Q2_nodes.size()) { - // std::cout << ", "; - // } - // } - // std::cout << std::endl; - - // for (size_t i = 0; auto& element : cdf_vector) { - // std::cout << element; - // // Print comma if it's not the last element - // if (++i != cdf_vector.size()) { - // std::cout << ", "; - // } - // } - // std::cout << std::endl; - - // for (size_t i = 0; auto& element : pdf_vector) { - // std::cout << element; - // // Print comma if it's not the last element - // if (++i != pdf_vector.size()) { - // std::cout << ", "; - // } - // } - // std::cout << std::endl; - // set the spline table siren::utilities::TableData1D inverse_cdf_data; inverse_cdf_data.x = cdf_vector; @@ -339,17 +309,11 @@ void CharmMesonDecay::computeDiffGammaCDF(std::vector constants, double void CharmMesonDecay::SampleFinalState(dataclasses::CrossSectionDistributionRecord & record, std::shared_ptr random) const { // first obtain the constants needed for further computation from the signature - // std::cout<<"b1"< constants = FormFactorFromRecord(record); double mD = particleMass(record.signature.primary_type); double mK = particleMass(record.signature.secondary_types[0]); - // std::cout << "input masses: " << mD << " " << mK << std::endl; - // first sample a q^2 - //////////////////////////////////////////// - // computeDiffGammaCDF(constants, mD, mK);// - //////////////////////////////////////////// double rand_value_for_Q2 = random->Uniform(0, 1); double Q2 = inverseCdf(rand_value_for_Q2); @@ -358,8 +322,6 @@ void CharmMesonDecay::SampleFinalState(dataclasses::CrossSectionDistributionReco double sinTheta = std::sin(std::acos(cosTheta)); // set the x axis to be the D direction geom3::UnitVector3 x_dir = geom3::UnitVector3::xAxis(); - // std::cout<<"b2"<W+K/Pi decay, now treat the W->l+nu decay double ml = particleMass(record.signature.secondary_types[1]); double mnu = 0; @@ -400,11 +359,6 @@ void CharmMesonDecay::SampleFinalState(dataclasses::CrossSectionDistributionReco rk::P4 p4l_Wrest(P * geom3::Vector3(W_cosTheta, W_sinTheta, 0), ml); rk::P4 p4nu_Wrest(P * geom3::Vector3(-W_cosTheta, -W_sinTheta, 0), 0); - // std::cout << "momentums: " << p4l_Wrest << " " << p4nu_Wrest << std::endl; - // std::cout << "check mass of l and nu: " << p4l_Wrest.m() << " " << p4nu_Wrest.m() << std::endl; - //now rotate so they are defined wrt the lab frame W direction - // std::cout<<"b5"< & secondaries = record.GetSecondaryParticleRecords(); siren::dataclasses::SecondaryParticleRecord & kpi = secondaries[0]; siren::dataclasses::SecondaryParticleRecord & lepton = secondaries[1]; @@ -434,56 +387,6 @@ void CharmMesonDecay::SampleFinalState(dataclasses::CrossSectionDistributionReco neutrino.SetMass(p4nu_lab.m()); neutrino.SetHelicity(record.primary_helicity); - // finally, we can populate the record, for implementation in prometheus, maybe add treatment of hadrons, but could be implemnted on p side - // record.secondary_momenta.resize(3); - // record.secondary_masses.resize(3); - // record.secondary_helicity.resize(3); // 0 is the hadron, 1 is the lepton, 2 is the neutrino - // // the K/pi - // record.secondary_momenta[0][0] = p4K_lab.e(); - // record.secondary_momenta[0][1] = p4K_lab.px(); - // record.secondary_momenta[0][2] = p4K_lab.py(); - // record.secondary_momenta[0][3] = p4K_lab.pz(); - // record.secondary_masses[0] = p4K_lab.m(); - // record.secondary_helicity[0] = 0; - // // the lepton - // record.secondary_momenta[1][0] = p4l_lab.e(); - // record.secondary_momenta[1][1] = p4l_lab.px(); - // record.secondary_momenta[1][2] = p4l_lab.py(); - // record.secondary_momenta[1][3] = p4l_lab.pz(); - // record.secondary_masses[1] = p4l_lab.m(); - // record.secondary_helicity[1] = 1; - // // the neutrino - // record.secondary_momenta[2][0] = p4nu_lab.e(); - // record.secondary_momenta[2][1] = p4nu_lab.px(); - // record.secondary_momenta[2][2] = p4nu_lab.py(); - // record.secondary_momenta[2][3] = p4nu_lab.pz(); - // record.secondary_masses[2] = p4nu_lab.m(); - // record.secondary_helicity[2] = 1; - - //for debug purposes - // double p4w_rest_Q2 = pow(p4W_Drest.e(), 2) - pow(p4W_Drest.px(), 2) - - // pow(p4W_Drest.py(), 2) - pow(p4W_Drest.pz(), 2); - // double p4w_lab_Q2 = pow(p4W_lab.e(), 2) - pow(p4W_lab.px(), 2) - - // pow(p4W_lab.py(), 2) - pow(p4W_lab.pz(), 2); - // std::cout << p4W_Drest.e() << " " << p4W_lab.e() << " " << PW << " " << p4W_Drest.p() << " " << p4W_Drest << std::endl; - // std::cout << p4K_Drest.e() << " " << p4K_lab.e() << " " << PK << " " << p4K_Drest.p() << " " << p4K_Drest << std::endl; - // std::cout << Q2 << " " << sqrt(Q2)<< std::endl; - // std::cout << "invariant mass of the W in two frames are " << p4w_lab_Q2 << " " << p4w_rest_Q2 << std::endl; - // std::cout << "check mass of W: " << p4W_lab.m() << " " << p4W_Drest.m() << std::endl; - // std::cout << "check mass of K: " << p4K_lab.m() << " " << p4K_Drest.m() << std::endl; - - - - // rk::P4 inv_mass_Wrest = p4l_Wrest + p4nu_Wrest; - // rk::P4 inv_mass_lab = p4l_lab + p4nu_lab; - // std::cout << "inv masses in two frames: " << pow(inv_mass_Wrest.m(), 2) << " " << pow(inv_mass_lab.m(), 2) << std::endl; - // std::cout << "using inv mass calculator: " << pow(invMass(p4l_Wrest, p4nu_Wrest), 2) << " " << pow(invMass(p4l_lab, p4nu_lab), 2) << std::endl; - // std::cout << "energy of l and nu and inv mass: " << El << " " << Enu << " " << pow(El+Enu, 2) << std::endl; - // std::cout << "momentums: " << p4l_Wrest << " " << p4nu_Wrest << std::endl; - // std::cout << "check mass of l and nu: " << p4l_Wrest.m() << " " << p4nu_Wrest.m() << std::endl; - // std::cout << "W energy in rest frame: " << pow(p4W_lab.m(), 2) << std::endl; - - } double CharmMesonDecay::FinalStateProbability(dataclasses::InteractionRecord const & record) const { diff --git a/projects/interactions/private/DMesonELoss.cxx b/projects/interactions/private/DMesonELoss.cxx index 867f7edf..5779d77f 100644 --- a/projects/interactions/private/DMesonELoss.cxx +++ b/projects/interactions/private/DMesonELoss.cxx @@ -14,9 +14,6 @@ #include // for P4, Boost #include // for Vector3 -#include // for splinetable -//#include - #include "SIREN/interactions/CrossSection.h" // for CrossSection #include "SIREN/dataclasses/InteractionRecord.h" // for Interactio... #include "SIREN/dataclasses/Particle.h" // for Particle @@ -131,18 +128,9 @@ double DMesonELoss::DifferentialCrossSection(dataclasses::InteractionRecord cons // rk::P4 p2(geom3::Vector3(interaction.target_momentum[1], interaction.target_momentum[2], interaction.target_momentum[3]), interaction.target_mass); double primary_energy; rk::P4 p1_lab; - // rk::P4 p2_lab; - // if(interaction.target_momentum[1] == 0 and interaction.target_momentum[2] == 0 and interaction.target_momentum[3] == 0) { primary_energy = interaction.primary_momentum[0]; p1_lab = p1; - // p2_lab = p2; - // } else { - // rk::Boost boost_start_to_lab = p2.restBoost(); - // p1_lab = boost_start_to_lab * p1; - // p2_lab = boost_start_to_lab * p2; - // primary_energy = p1_lab.e(); - // std::cout << "D Meson Diff Xsec: not in lab frame???" << std::endl; - // } + double final_energy = interaction.secondary_momenta[0][0]; double z = 1 - final_energy / primary_energy; @@ -170,8 +158,6 @@ double DMesonELoss::InteractionThreshold(dataclasses::InteractionRecord const & void DMesonELoss::SampleFinalState(dataclasses::CrossSectionDistributionRecord& interaction, std::shared_ptr random) const { - // std::cout << "In D Meson E Loss Sample Final State" << std::endl; - rk::P4 p1(geom3::Vector3(interaction.primary_momentum[1], interaction.primary_momentum[2], interaction.primary_momentum[3]), interaction.primary_mass); // rk::P4 p2(geom3::Vector3(interaction.target_momentum[1], interaction.target_momentum[2], interaction.target_momentum[3]), interaction.target_mass); @@ -184,24 +170,8 @@ void DMesonELoss::SampleFinalState(dataclasses::CrossSectionDistributionRecord& double primary_energy; double Dmass = interaction.primary_mass; rk::P4 p1_lab; - // rk::P4 p2_lab; - // if(interaction.target_momentum[1] == 0 and interaction.target_momentum[2] == 0 and interaction.target_momentum[3] == 0) { p1_lab = p1; - // p2_lab = p2; primary_energy = p1_lab.e(); - // } else { - // // this is currently not implemented - // // Rest frame of p2 will be our "lab" frame - // rk::Boost boost_start_to_lab = p2.restBoost(); - // p1_lab = boost_start_to_lab * p1; - // p2_lab = boost_start_to_lab * p2; - // primary_energy = p1_lab.e(); - // // std::cout << "D Meson Energy Loss: not in lab frame???" << std::endl; - // } - // following line is wrong but i dont want to change it now fuck it. - // std::cout << " " << interaction.primary_momentum[0] << " " << interaction.primary_momentum[1] << " " << interaction.primary_momentum[2] << " " << interaction.primary_momentum[3]; - // std::cout << primary_energy << " " << pow(primary_energy, 2) - pow(Dmass, 2) << " " << - // sqrt(pow(interaction.primary_momentum[1], 2) +pow(interaction.primary_momentum[2], 2) +pow(interaction.primary_momentum[3], 2)) << std::endl; // sample an inelasticity from gaussian using Box-Muller Transform double sigma = 0.2; double z0 = 0.56; // for mesons only, for baryons it's 0.59, but not implemented yet @@ -218,8 +188,6 @@ void DMesonELoss::SampleFinalState(dataclasses::CrossSectionDistributionRecord& while (u1 == 0); u2 = random->Uniform(0, 1); double z = sigma * sqrt(-2.0 * log(u1)) * cos(2.0 * M_PI * u2) + z0; - // std::cout << z<< std::endl; - // now modify the energy of the charm hadron and the corresponding momentum final_energy = primary_energy * (1-z); if (pow(final_energy, 2) - pow(Dmass, 2) >= 0) { @@ -228,14 +196,9 @@ void DMesonELoss::SampleFinalState(dataclasses::CrossSectionDistributionRecord& accept = false; } } while (!accept); - // this might be an infinite loop??????? - // need to check if the cross section length is good enough, how to make some relevant plots? - // std::cout << final_energy << std::endl; double p3f = sqrt(pow(final_energy, 2) - pow(Dmass, 2)); double p3i = std::sqrt(std::pow(p1.px(), 2) + std::pow(p1.py(), 2) + std::pow(p1.pz(), 2)); double p_ratio = p3f / p3i; - - // std::cout << " " << p3f << " " << p3i << " " << p_ratio << std::endl; rk::P4 pf(p_ratio * geom3::Vector3(p1.px(), p1.py(), p1.pz()), Dmass); std::vector & secondaries = interaction.GetSecondaryParticleRecords(); @@ -246,17 +209,6 @@ void DMesonELoss::SampleFinalState(dataclasses::CrossSectionDistributionRecord& dmeson.SetMass(pf.m()); dmeson.SetHelicity(interaction.primary_helicity); - // interaction.secondary_momenta.resize(1); - // interaction.secondary_masses.resize(1); - // interaction.secondary_helicity.resize(1); - - // interaction.secondary_momenta[0][0] = pf.e(); // p3_energy - // interaction.secondary_momenta[0][1] = pf.px(); // p3_x - // interaction.secondary_momenta[0][2] = pf.py(); // p3_y - // interaction.secondary_momenta[0][3] = pf.pz(); // p3_z - // interaction.secondary_masses[0] = pf.m(); - - // interaction.secondary_helicity[0] = interaction.primary_helicity; } double DMesonELoss::FinalStateProbability(dataclasses::InteractionRecord const & interaction) const { diff --git a/resources/Examples/DMesonExample/DIS_D.py b/resources/Examples/DMesonExample/DIS_D.py index 8fb483a9..c89d25b0 100644 --- a/resources/Examples/DMesonExample/DIS_D.py +++ b/resources/Examples/DMesonExample/DIS_D.py @@ -1,30 +1,38 @@ import os - +import numpy as np import siren from siren.SIREN_Controller import SIREN_Controller +import nuflux # Number of events to inject -events_to_inject = 5 +events_to_inject = 10000 # Expeirment to run experiment = "IceCube" +# physical flux model to use +physical_flux = "atmos" + # Define the controller controller = SIREN_Controller(events_to_inject, experiment, seed = 1) # Particle to inject primary_type = siren.dataclasses.Particle.ParticleType.NuMu -cross_section_model = "CSMSDISSplines" +xs_option = "" # current choices are the empty string and cutoff-"" +xsfiledir = "/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/miaochenjin/CharmXS/xsec_splines/M_Muon-{}105MeV".format(xs_option) -xsfiledir = siren.utilities.get_cross_section_model_path(cross_section_model) +if xs_option == "": + spline_option = "M_Muon-105MeV" # this is to account for the fact that I put the output names of this particular spline incorrectly +else: + spline_option = "" # Cross Section Model target_type = siren.dataclasses.Particle.ParticleType.Nucleon DIS_xs = siren.interactions.CharmDISFromSpline( - os.path.join(xsfiledir, "dsdxdy_nu_CC_iso.fits"), - os.path.join(xsfiledir, "sigma_nu_CC_iso.fits"), + os.path.join(xsfiledir, "{}dsdxdynu-N-cc-HERAPDF15LO_EIG_central.fits".format(spline_option)), + os.path.join(xsfiledir, "{}sigmanu-N-cc-HERAPDF15LO_EIG_central.fits".format(spline_option)), [primary_type], [target_type], "m" ) @@ -41,9 +49,35 @@ primary_physical_distributions["mass"] = mass_dist # energy distribution -edist = siren.distributions.PowerLaw(2, 1e5, 1e10) +edist = siren.distributions.PowerLaw(2, 1e4, 1e10) primary_injection_distributions["energy"] = edist -primary_physical_distributions["energy"] = edist + +# make an atmospheric flux +flux = nuflux.makeFlux('honda2006') +nu_type=nuflux.NuMu +erange = np.logspace(2,6,100) +erange_atmo = np.logspace(2,6,100) +cosrange = np.linspace(0,1,100) +atmo_flux_tables = {} +particle = nuflux.NuMu +siren_key = siren.dataclasses.Particle.ParticleType(int(particle)) +atmo_flux_tables[siren_key] = np.zeros(len(erange)) +for i,e in enumerate(erange): + f = flux.getFlux(particle,e,cosrange) + atmo_flux_tables[siren_key][i] += 0.01*np.sum(f) * 1e4 * 2 * np.pi + +# this is for weighting the events to the astrophysical flux +if physical_flux == "astro": + edist_astro = siren.distributions.PowerLaw(2, 1e4, 1e10) + norm = 1e-18 * 1e4 * 4 * np.pi # GeV^-1 m^-2 s^-1 + edist_astro.SetNormalizationAtEnergy(norm,1e5) + primary_physical_distributions["energy"] = edist_astro +elif physical_flux == "atmos": + edist_atmo = siren.distributions.TabulatedFluxDistribution(erange_atmo,atmo_flux_tables[primary_type],True) + primary_physical_distributions["energy"] = edist_atmo +else: + primary_injection_distributions["energy"] = edist + # direction distribution direction_distribution = siren.distributions.IsotropicDirection() @@ -77,7 +111,7 @@ def add_secondary_to_controller(controller, secondary_type, secondary_xsecs, sec return secondary_collection -# secondary interactions +# # secondary interactions charms = siren.dataclasses.Particle.ParticleType.Charm DPlus = siren.dataclasses.Particle.ParticleType.DPlus D0 = siren.dataclasses.Particle.ParticleType.D0 @@ -90,14 +124,15 @@ def add_secondary_to_controller(controller, secondary_type, secondary_xsecs, sec secondary_DPlus_collection = add_secondary_to_controller(controller, DPlus, D_energy_loss, DPlus_decay) secondary_D0_collection = add_secondary_to_controller(controller, D0, D_energy_loss, D0_decay) +# secondary_DPlus_collection = add_secondary_to_controller(controller, DPlus, DPlus_decay) +# secondary_D0_collection = add_secondary_to_controller(controller, D0, D0_decay) + controller.SetInteractions(primary_xs, [secondary_charm_collection, secondary_D0_collection, secondary_DPlus_collection]) +# controller.SetInteractions(primary_xs, [secondary_charm_collection]) +# controller.SetInteractions(primary_xs, []) controller.Initialize() -# def stop(datum, i): -# secondary_type = datum.record.signature.secondary_types[i] -# return ((secondary_type != siren.dataclasses.Particle.ParticleType.Charm) and (secondary_type != siren.dataclasses.Particle.ParticleType.DPlus)) - def stop(datum, i): return False @@ -105,6 +140,12 @@ def stop(datum, i): events = controller.GenerateEvents() -os.makedirs("output", exist_ok=True) +print("finished generating events") + +outdir = "/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/miaochenjin/DBSearch/SIREN_outputs" +expname = "0819_LE_debug_CharmHadron_atmos" +savedir = os.path.join(outdir, expname) + +os.makedirs(savedir, exist_ok=True) -controller.SaveEvents("output/FullSim") +controller.SaveEvents("{}/{}_".format(savedir, expname)) diff --git a/resources/Examples/DMesonExample/make_plots.py b/resources/Examples/DMesonExample/make_plots.py deleted file mode 100644 index 52d4dad9..00000000 --- a/resources/Examples/DMesonExample/make_plots.py +++ /dev/null @@ -1,250 +0,0 @@ -import h5py -import numpy as np -from matplotlib import pyplot as plt -from matplotlib.colors import LogNorm -from mpl_toolkits.axes_grid1 import make_axes_locatable -from parse_output import analysis -import os - -pathname = "/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/miaochenjin/DBSearch/SIREN_outputs/" -# parquetname = "0708_test/0708_test_.parquet" -expname = "0709_astro_flux" -parquetname = "{}/{}_.parquet".format(expname, expname) - - -filename = os.path.join(pathname, parquetname) -savedir = os.path.join(pathname, "plots/") - -plt.style.use('paper.mplstyle') - -sim = analysis(filename) - - -c = 3 * 1e8 # m/s -m_D0 = 1.86962 # GeV -m_Dp = 1.86484 -t_Dp = 1040 * 1e-15 # s -t_D0 = 410 * 1e-15 -m_ice = 18.02 # g/mol -N = 6.02214 * 1e23 #mol^-1 -rho = 0.917 # g/cm^3 - - -def normalize(hist, xbins, ybins): - normed_hist = np.zeros_like(hist) - for i in range(len(xbins) - 1): - tot = 0 - for j in range(len(ybins) - 1): - tot += hist[i][j] - for j in range(len(ybins) - 1): - if tot != 0: - normed_hist[i][j] = hist[i][j] / tot - else: - normed_hist[i][j] = 0 - return normed_hist - -def analytic_decay_length(E, t, m): - return E * t / ((m/(c ** 2)) * c) - -def xsec(E): - return (np.exp(1.891 + 0.2095 * np.log10(E)) - 2.157 + 1.263 * np.log10(E)) * 1e-27 # convert to cm^2 - -def analytic_free_path(E): - return (m_ice / (rho * N * xsec(E))) / 100 # convert to m - -def plot_separation_distribution(analysis_, dim = 2): - D0_energies, D0_separations, Dp_energies, Dp_separations, D0_weights, Dp_weights = analysis_.separation_analysis() - min_eng = 1e1 - max_eng = 1e9 - energy_bins = np.logspace(np.log10(min_eng), np.log10(max_eng), 20) - - energy_bins_centers = np.zeros((len(energy_bins) - 1,)) - for i in range(len(energy_bins_centers)): - energy_bins_centers[i] = np.sqrt(energy_bins[i] * energy_bins[i + 1]) - D0_analytic_lengths = analytic_decay_length(energy_bins_centers, t_D0, m_D0) - Dp_analytic_lengths = analytic_decay_length(energy_bins_centers, t_Dp, m_Dp) - - min_sep = 1e-3 - max_sep = 50000 - log_separation_bins = np.logspace(np.log10(min_sep), np.log10(max_sep), 20) - sep_bin_widths = np.sqrt(log_separation_bins[1:] * log_separation_bins[:-1]) - - if dim == 2: - - X2, Y2 = np.meshgrid(energy_bins, log_separation_bins) - log_hist_D0, _, _ = np.histogram2d(D0_energies, D0_separations, bins = (energy_bins, log_separation_bins), weights = D0_weights) - log_hist_Dp, _, _ = np.histogram2d(Dp_energies, Dp_separations, bins = (energy_bins, log_separation_bins), weights = Dp_weights) - log_hist_D0 = normalize(log_hist_D0, energy_bins, log_separation_bins) - log_hist_Dp = normalize(log_hist_Dp, energy_bins, log_separation_bins) - - fig, axes = plt.subplots(nrows = 1, ncols = 2, figsize = (11, 5)) - - log_im1 = axes[0].pcolor(X2, Y2, log_hist_D0.T, cmap="plasma", alpha = 0.7, vmin=0, vmax=1) - log_im2 = axes[1].pcolor(X2, Y2, log_hist_Dp.T, cmap="plasma", alpha = 0.7, vmin=0, vmax=1) - - # divider1 = make_axes_locatable(axes[0]) - # cax1 = divider1.append_axes('right', size='5%', pad=0.05) - divider2 = make_axes_locatable(axes[1]) - cax2 = divider2.append_axes('right', size='5%', pad=0.05) - fig.colorbar(log_im2, cax=cax2, orientation='vertical', alpha = 0.7) - # fig.colorbar(log_im1, cax=cax1, orientation='vertical', alpha = 0.7) - - axes[0].set_title(r"$D^0$ Separation") - axes[1].set_title(r"$D^+$ Separation") - - axes[0].set_xlabel(r"$E_{D^0}$ [GeV]") - axes[1].set_xlabel(r"$E_{D^+}$ [GeV]") - - axes[0].set_ylabel("Separation Length [m]") - - axes[0].set_xscale('log') - axes[1].set_xscale('log') - axes[0].set_yscale('log') - axes[1].set_yscale('log') - - axes[0].set_ylim(min_sep, max_sep) - axes[1].set_ylim(min_sep, max_sep) - axes[0].set_xlim(min_eng, max_eng) - axes[1].set_xlim(min_eng, max_eng) - - # also plot the analytic lines - axes[0].plot(energy_bins_centers, D0_analytic_lengths, color = '#FEF3E8', alpha = 0.7) - axes[1].plot(energy_bins_centers, Dp_analytic_lengths, label = r"$d = \frac{E \tau}{mc}$", color = '#FEF3E8', alpha = 0.7) - - legend = axes[1].legend(loc = 'upper left') - for text in legend.get_texts(): - text.set_color('#FEF3E8') - - savename = os.path.join(savedir, "Separation_Length_Distribution") - - elif dim == 1: - - D0_separations.extend(Dp_weights) - D0_weights.extend(Dp_weights) - - sep_hist, _ = np.histogram(D0_separations, log_separation_bins, weights = D0_weights) - fig, ax = plt.subplots(1, 1, figsize = (8, 6)) - ax.hist(log_separation_bins[:-1], bins = log_separation_bins, weights = sep_hist / sep_bin_widths, \ - label = r"$\phi \sim E^{-2}$", \ - alpha = 0.9, color = '#D06C9D', histtype = 'step') - - savename = os.path.join(savedir, "PowerLaw_2_Separation_Length_Distribution") - ax.legend(loc = 'upper right') - ax.set_xscale('log') - ax.set_yscale('log') - ax.set_xlabel(r'$d_{\textrm{Sep}}$ [m]') - ax.set_ylabel('Normalized Weights Event Count') - - fig.savefig(savename, bbox_inches = 'tight') - -plot_separation_distribution(sim, dim = 1) -plot_separation_distribution(sim, dim = 2) - - -def plot_2d_energy_loss(analysis_): - E_D0, E_Dp, n_D0, n_Dp = analysis_.energy_loss_analysis_2d() - energy_bins = np.logspace(np.log10(min(min(E_D0), min(E_Dp))), np.log10(max(max(E_D0), max(E_Dp))), 20) - num_bins = np.linspace(-0.01, 7.99, 9) - X, Y = np.meshgrid(energy_bins, num_bins) - hist_D0, _, _ = np.histogram2d(E_D0, n_D0, bins = (energy_bins, num_bins)) - hist_Dp, _, _ = np.histogram2d(E_Dp, n_Dp, bins = (energy_bins, num_bins)) - - hist_D0 = normalize(hist_D0, energy_bins, num_bins) - hist_Dp = normalize(hist_Dp, energy_bins, num_bins) - - fig, axes = plt.subplots(nrows = 1, ncols = 2, figsize = (11, 5)) - im1 = axes[0].pcolor(X, Y, hist_D0.T, cmap="plasma", alpha = 0.7) - im2 = axes[1].pcolor(X, Y, hist_Dp.T, cmap="plasma", alpha = 0.7) - divider2 = make_axes_locatable(axes[1]) - cax2 = divider2.append_axes('right', size='5%', pad=0.05) - fig.colorbar(im2, cax=cax2, orientation='vertical', alpha = 0.7) - axes[0].axvline(x = 53 * 1e3, color = '#FEF3E8', alpha = 0.7, label = r"$d_{D^0} = l_{D^0}$") - axes[1].axvline(x = 22 * 1e3, color = '#FEF3E8', alpha = 0.7, label = r"$d_{D^+} = l_{D^+}$") - - axes[0].set_title(r"$D^0-p$ Collision") - axes[1].set_title(r"$D^+-p$ Collision") - - axes[0].set_xlabel(r"$E_{D^0}$ [GeV]") - axes[1].set_xlabel(r"$E_{D^+}$ [GeV]") - - axes[0].set_ylim(0, 8) - axes[1].set_ylim(0, 8) - - axes[0].set_ylabel(r"$n_{\textrm{Elastic Collision}}$") - axes[0].set_xscale('log') - axes[1].set_xscale('log') - legend0 = axes[0].legend() - legend1 = axes[1].legend() - for text in legend0.get_texts(): - text.set_color('#FEF3E8') - for text in legend1.get_texts(): - text.set_color('#FEF3E8') - - - fig.savefig("./plots/Energy_loss_2d_Distribution", bbox_inches = 'tight') - -# plot_2d_energy_loss(sim) -# exit(0) - -def plot_free_path_distribution(): - D0_E_list, D0_free_path_list, Dp_E_list, Dp_free_path_list = analysis_.free_path_analysis() - energy_bins = np.logspace(1.5, 9, 20) - distance_bins = np.logspace(-3, np.log10(5000), 20) - X, Y = np.meshgrid(energy_bins, distance_bins) - hist_D0, _, _ = np.histogram2d(D0_E_list, D0_free_path_list, bins = (energy_bins, distance_bins)) - hist_Dp, _, _ = np.histogram2d(Dp_E_list, Dp_free_path_list, bins = (energy_bins, distance_bins)) - - hist_D0 = normalize(hist_D0, energy_bins, distance_bins) - hist_Dp = normalize(hist_Dp, energy_bins, distance_bins) - - energy_bins_centers = np.zeros((len(energy_bins) - 1,)) - for i in range(len(energy_bins_centers)): - energy_bins_centers[i] = np.sqrt(energy_bins[i] * energy_bins[i + 1]) - D0_analytic_lengths = analytic_free_path(energy_bins_centers) - Dp_analytic_lengths = analytic_free_path(energy_bins_centers) - - - fig, axes = plt.subplots(nrows = 1, ncols = 2, figsize = (11, 5)) - im1 = axes[0].pcolor(X, Y, hist_D0.T, cmap="plasma", alpha = 0.7, vmin=0, vmax=1) - im2 = axes[1].pcolor(X, Y, hist_Dp.T, cmap="plasma", alpha = 0.7, vmin=0, vmax=1) - divider2 = make_axes_locatable(axes[1]) - cax2 = divider2.append_axes('right', size='5%', pad=0.05) - fig.colorbar(im2, cax=cax2, orientation='vertical', alpha = 0.7) - - axes[0].set_title(r"$D^0-p$ Free Path") - axes[1].set_title(r"$D^+-p$ Free Path") - - axes[0].set_xlabel(r"$E_{D^0}$ [GeV]") - axes[1].set_xlabel(r"$E_{D^+}$ [GeV]") - - axes[0].plot(energy_bins_centers, D0_analytic_lengths, color = '#FEF3E8', alpha = 0.7) - axes[1].plot(energy_bins_centers, Dp_analytic_lengths, label = r"$l = \frac{m_{\textrm{ice}}}{\rho N_A \sigma(E)}$", color = '#FEF3E8', alpha = 0.7) - - # also plot the decay lengths to explain low energy increase - D0_decay_analytic_lengths = analytic_decay_length(energy_bins_centers, t_D0, m_D0) - Dp_decay_analytic_lengths = analytic_decay_length(energy_bins_centers, t_Dp, m_Dp) - - axes[0].plot(energy_bins_centers, D0_decay_analytic_lengths, label = r"$d = \frac{E \tau}{mc}$", color = '#A597B6', alpha = 0.7) - axes[1].plot(energy_bins_centers, Dp_decay_analytic_lengths, color = '#A597B6', alpha = 0.7) - - axes[0].set_ylabel(r"$l_{\textrm{Free}}$") - axes[0].set_xscale('log') - axes[1].set_xscale('log') - axes[0].set_yscale('log') - axes[1].set_yscale('log') - - axes[0].set_ylim(1e-3, 5000) - axes[1].set_ylim(1e-3, 5000) - - legend0 = axes[0].legend(loc = 'upper left') - for text in legend0.get_texts(): - text.set_color('#A597B6') - - legend1 = axes[1].legend(loc = 'upper left') - for text in legend1.get_texts(): - text.set_color('#FEF3E8') - - fig.savefig("./plots/Free_Path_Distribution", bbox_inches = 'tight') - return - -# plot_free_path_distribution() \ No newline at end of file