//____________________________________________________________________________ /*! \class genie::GGenerateEvent.cxx \brief A base generation driver for an under-water neutrino telescope \author Carla Distefano LNS-INFN, Catania Included modifications by Alfonso Garcia NIKHEF, Amsterdam (January 2019) Included JPP propagation by Victor Carretero IFIC, Valencia (October 2020) \created July 19, 2017 \cpright Copyright (c) 2017-2022, The KM3NeT Collaboration For the full text of the license see $GSEAGEN/LICENSE */ //____________________________________________________________________________ #include #include #ifdef _GENIEVERSIONGTEQ3__ #include "Framework/Conventions/Constants.h" #else #include "Conventions/Constants.h" #endif #include "SeaNuDrivers/GGenerateEvent.h" #include "utils/matrix_operations.h" #include #ifdef _MUSIC_ENABLED__ extern "C" { void initialize_music_sr_(int *muin);} extern "C" { void initialize_music_sw_(int *muin);} extern "C" { void muon_transport_sr_(double *VX,double *VY,double *VZ,double *CX,double *CY,double *CZ,double *E,double *Depth,double *T,double *Rho,int *flag1,int *flag2,double *VXF,double *VYF,double *VZF,double *CXF,double *CYF,double *CZF,double *EF,double *DepthF,double *TF);} extern "C" { void muon_transport_sw_(double *VX,double *VY,double *VZ,double *CX,double *CY,double *CZ,double *E,double *Depth,double *T,double *Rho,int *flag1,int *flag2,double *VXF,double *VYF,double *VZF,double *CXF,double *CYF,double *CZF,double *EF,double *DepthF,double *TF);} #endif #ifdef _HETAU_ENABLED__ #include "Tauola/Tauola.h" #include "Tauola/TauolaHEPEVTParticle.h" using namespace Tauolapp; int TAUMODEL=1; int ITFLAG=1; double lifetime = 2.91e-13; //s double cSpeed = kLightSpeed/(units::meter/units::second); extern "C" { void initialize_tausic_sr_(int *tauin);} extern "C" { void initialize_tausic_sw_(int *tauin);} extern "C" { void tau_transport_sr_(double *VX,double *VY,double *VZ,double *CX,double *CY,double *CZ,double *E,double *Depth,double *T,double *Rho,int *flag1,int *flag2,double *VXF,double *VYF,double *VZF,double *CXF,double *CYF,double *CZF,double *EF,double *DepthF,double *TF, int *IDEC, int *ITFLAG, double *TAUTI, double *TAUTF);} extern "C" { void tau_transport_sw_(double *VX,double *VY,double *VZ,double *CX,double *CY,double *CZ,double *E,double *Depth,double *T,double *Rho,int *flag1,int *flag2,double *VXF,double *VYF,double *VZF,double *CXF,double *CYF,double *CZF,double *EF,double *DepthF,double *TF, int *IDEC, int *ITFLAG, double *TAUTI, double *TAUTF);} #endif //___________________________________________________________________________ void copy( const GSeaTrack& from , PropaTrack& to ) { to.Id = from.Id; to.Status = from.Status; to.Pdg = from.PdgCode; to.MotherId = from.MotherId; to.D1 = from.D1; to.D2 = from.D2; to.D3 = from.D3; to.V1 = from.V1; to.V2 = from.V2; to.V3 = from.V3; to.E = from.E; to.Length = from.Length; to.T = from.T; } //___________________________________________________________________________ void copy( const PropaTrack& from, GSeaTrack& to) { to.Id = from.Id; to.Status = from.Status; to.PdgCode = from.Pdg; to.MotherId = from.MotherId; to.D1 = from.D1; to.D2 = from.D2; to.D3 = from.D3; to.V1 = from.V1; to.V2 = from.V2; to.V3 = from.V3; to.E = from.E; to.Length = from.Length; to.T = from.T; } //___________________________________________________________________________ GGenerateEvent::GGenerateEvent(GenParam * GenPar) { fGenPar = GenPar; fSeaEvent = new GSeaEvent(); fNEvt=0; fNTot=0; EoF = false; // get flux driver flux_driver_diffuse=0; flux_driver_point=0; fFluxDriver=0; fFileDriver=0; #ifdef _ANTARES_ENABLED__ fAntaresEvt=0; #endif this->CalcMediaParams(); if(fGenPar->GenCode.compare("GENIE")==0){ if(fGenPar->GenMode.compare("DIFFUSE")==0){ flux_driver_diffuse = this->GetFluxDiffuse(); fFluxDriver=dynamic_cast(flux_driver_diffuse); } else if(fGenPar->GenMode.compare("POINT")==0) { flux_driver_point = this->GetFluxPoint(); fFluxDriver=dynamic_cast(flux_driver_point); } else { LOG("GGenerateEvent", pFATAL) << "Unknown generation mode "<GenMode<<", exiting..."; exit(1); } } else if(fGenPar->GenCode.compare("CORSIKA")==0){ if (fGenPar->GenMode.compare("BIN")==0) { LOG("GGenerateEvent", pDEBUG) <<"Generating an event from CORSIKA input"; GSeaCORSIKAFileFlux * corsika_file_driver = new GSeaCORSIKAFileFlux(fGenPar); fFileDriver = dynamic_cast(corsika_file_driver); } #ifdef _ANTARES_ENABLED__ else if(fGenPar->GenMode.compare("EVT")==0){ GSeaAntaresFileFlux * antares_file_driver = new GSeaAntaresFileFlux(fGenPar); fFileDriver = dynamic_cast(antares_file_driver); fAntaresEvt=antares_file_driver->GetEvent(); } #endif else { LOG("GGenerateEvent", pFATAL) << "Unknown generation mode "<GenMode<<", exiting..."; exit(1); } } if(fGenPar->TrackEvts>0){ // we initialize the muon propagator only if we use them if(fGenPar->PropCode.compare("PropaMuon")==0){ PropaMuonSW = new PropaMuon(fGenPar->SeaWaterComp,fGenPar->RhoSW); PropaMuonSR = new PropaMuon(fGenPar->RockComp,fGenPar->RhoSR); } #ifdef _MUSIC_ENABLED__ else if(fGenPar->PropCode.compare("MUSIC")==0){ string music_path; if(gSystem->Getenv("MUSICPATH"))music_path = gSystem->Getenv("MUSICPATH"); else { music_path = gSystem->Getenv("GSEAGEN"); music_path.append("/dat"); music_path = BOOST_PP_STRINGIZE(_MUSICDATA__); while(1){ size_t pos = music_path.find("\""); if(pos==string::npos)break; music_path.erase(pos,1); } } string cp_command("cp "); cp_command.append(music_path); cp_command.append("/music*.dat ."); system(cp_command.c_str()); int muin =2; initialize_music_sr_(&muin); initialize_music_sw_(&muin); } #endif #ifdef _PROPOSAL_ENABLED__ else if(fGenPar->PropCode.compare("PROPOSAL")==0){ this->ConfigProposal(); } #endif } #ifdef _HETAU_ENABLED__ if (fGenPar->TrackEvts==2) { Tauola::setSameParticleDecayMode(2); Tauola::setOppositeParticleDecayMode(2); } Tauola::initialize(); if (fGenPar->TrackEvts==3) Tauola::setTauBr(2,0.); string tausic_path; if(gSystem->Getenv("TAUSICPATH"))tausic_path = gSystem->Getenv("TAUSICPATH"); else { tausic_path = gSystem->Getenv("GSEAGEN"); tausic_path.append("/dat"); tausic_path = BOOST_PP_STRINGIZE(_TAUSICDATA__); while(1){ size_t pos = tausic_path.find("\""); if(pos==string::npos)break; tausic_path.erase(pos,1); } } string cp_command("cp "); cp_command.append(tausic_path); cp_command.append("/tausic*.dat ."); system(cp_command.c_str()); int tauin =2; initialize_tausic_sr_(&tauin); initialize_tausic_sw_(&tauin); #endif } //___________________________________________________________________________ GGenerateEvent::~GGenerateEvent(){ if(fFluxDriver)delete fFluxDriver; #ifdef _ANTARES_ENABLED__ if(fFileDriver)delete fFileDriver; #endif if(fGenPar->TrackEvts){ if(fGenPar->PropCode.compare("PropaMuon")==0){ if(PropaMuonSW)delete PropaMuonSW; if(PropaMuonSR)delete PropaMuonSR; } #ifdef _PROPOSAL_ENABLED__ if(fGenPar->PropCode.compare("PROPOSAL")==0){ if(ProposalMu)delete ProposalMu; } #endif } if(fSeaEvent)delete fSeaEvent; } //___________________________________________________________________________ void GGenerateEvent::CalcMediaParams(void){ // reads the medium compositions and calculate the molar masses, the number of nucleons and electrons of the media // fill the ComponentSW and ComponentSR if Proposal if enabled at compiling and choosen at run-times map AtomicMass; map ElementName; // load isotope table string filename = string(gSystem->Getenv("GENIE")) + string("/data/evgen/catalogues/iso/natural-isotopes.data"); bool is_accessible = ! (gSystem->AccessPathName( filename.c_str() )); if (!is_accessible) { LOG("gSeaGen", pFATAL) << "Can not read file: " << filename <<", exiting... "; exit(1); } // load the natural isotopes .txt file string input_buf; ifstream input(filename.c_str()); if (input.is_open()){ //skip first 8 lines (comments) for(int i=0; i<8; i++){ string buffer; getline(input, buffer); } int Z = -1, Z_previous = -1, nelements = 0, pdgcode = 0; double atomicmass = 0, abundance = 0; string elementname, subelementname; while( !input.eof() ) { //read in naturally occuring element info input >> Z; input >> elementname; input >> nelements; // check not re-reading same element if(Z!=Z_previous){ LOG("NatIsotop", pDEBUG) << "Reading entry for Z = " << Z; for(int n=0 ; n < nelements; n++){ input >> subelementname; input >> pdgcode; input >> atomicmass; input >> abundance; AtomicMass.insert(map::value_type(pdgcode,atomicmass)); ElementName.insert(map::value_type(pdgcode,elementname)); } } Z_previous = Z; } //!eof } else { LOG("gSeaGen", pFATAL) << "Can not open file: " << filename <<", exiting... "; exit(1); } // Br80 is not in the table but used in the default compositions, we add by hand if(AtomicMass.find(1000350800) == AtomicMass.end()){ AtomicMass.insert(map::value_type(1000350800,79.9185293)); ElementName.insert(map::value_type(1000350800,"Br")); } // SeaWater int TgtPdgCode; double TgtMassFrac; double AtomInMolecule; double MoleMin,mass; int nAtoms; double TgtMassFracMax=0.; double MoliFracMass=1.; vector AtomMass,Mole,MassFrac; vector Z,A; vector AtomName; #ifdef _JPP_ENABLED__ JPHYSICS::JRadiationSource::source_type EErad(&JPHYSICS::JRadiation::TotalCrossSectionEErad, &JPHYSICS::JRadiation::EfromEErad); JPHYSICS::JRadiationSource::source_type Brems(&JPHYSICS::JRadiation::TotalCrossSectionBrems, &JPHYSICS::JRadiation::EfromBrems); JPHYSICS::JRadiationSource::source_type GNrad(&JPHYSICS::JRadiation::TotalCrossSectionGNrad, &JPHYSICS::JRadiation::EfromGNrad); JPHYSICS::JRadiationSource::source_type Ion (&JPHYSICS::JRadiation::CalculateACoeff , NULL); #endif map::iterator iter = fGenPar->SeaWaterComp.begin(); for ( ; iter != fGenPar->SeaWaterComp.end(); ++iter) { TgtPdgCode=iter->first; TgtMassFrac=iter->second; if(AtomicMass.find(TgtPdgCode) != AtomicMass.end()){ if(TgtMassFrac>TgtMassFracMax){ TgtMassFracMax=TgtMassFrac; MoliFracMass=TgtMassFrac/AtomicMass.find(TgtPdgCode)->second; } MassFrac.push_back(TgtMassFrac); mass=AtomicMass.find(TgtPdgCode)->second; AtomMass.push_back(AtomicMass.find(TgtPdgCode)->second); AtomName.push_back(ElementName.find(TgtPdgCode)->second); A.push_back(pdg::IonPdgCodeToA(TgtPdgCode)); Z.push_back(pdg::IonPdgCodeToZ(TgtPdgCode)); Mole.push_back(TgtMassFrac/mass); } else { LOG("gSeaGen", pFATAL) << "Element "<< TgtPdgCode<<" not found in file " << filename <<", exiting... "; exit(1); } } MoleMin=*min_element(Mole.begin(),Mole.end()); fGenPar->MolarMassSW=0.; fGenPar->NNucleonsSW=0; fGenPar->NElectronsSW=0; for(unsigned iComp=0; iCompMolarMassSW+=(double)nAtoms*AtomMass[iComp]; fGenPar->NNucleonsSW+=nAtoms*A[iComp]; fGenPar->NElectronsSW+=nAtoms*Z[iComp]; AtomInMolecule=MassFrac[iComp]/AtomMass[iComp]/MoliFracMass; #ifdef _PROPOSAL_ENABLED__ if(fGenPar->PropCode.compare("PROPOSAL")==0)ComponentSW.push_back(PROPOSAL::Components::Component(AtomName[iComp].c_str(),Z[iComp],AtomMass[iComp],AtomInMolecule)); #endif #ifdef _JPP_ENABLED__ if(fGenPar->PropCode.compare("JPP")==0){ JPHYSICS::JRadiation jppelement(Z[iComp],A[iComp],40,0.01,0.1,0.1); // (Z,A,integration steps, minimal energy for brems, minimal energy for epair, minimal energy for nuclear interaction) JLANG::JSharedPointer jppEl(new JPHYSICS::JRadiationFunction(jppelement, 300, 0.2, 1.0e11)); //(element, number of energy bins, min energy, max energy) radiationSW.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSW * MassFrac[iComp], EErad)); radiationSW.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSW * MassFrac[iComp], Brems)); radiationSW.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSW * MassFrac[iComp], GNrad)); ionizationSW.push_back(new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSW * MassFrac[iComp], Ion)); } #endif } // Rock TgtMassFracMax=0.; MoliFracMass=1.; vector().swap(AtomMass); vector().swap(Mole); vector().swap(MassFrac); vector().swap(Z); vector().swap(A); vector().swap(AtomName); iter = fGenPar->RockComp.begin(); for ( ; iter != fGenPar->RockComp.end(); ++iter) { TgtPdgCode=iter->first; TgtMassFrac=iter->second; if(AtomicMass.find(TgtPdgCode) != AtomicMass.end()){ if(TgtMassFrac>TgtMassFracMax){ TgtMassFracMax=TgtMassFrac; MoliFracMass=TgtMassFrac/AtomicMass.find(TgtPdgCode)->second; } MassFrac.push_back(TgtMassFrac); mass=AtomicMass.find(TgtPdgCode)->second; AtomMass.push_back(AtomicMass.find(TgtPdgCode)->second); AtomName.push_back(ElementName.find(TgtPdgCode)->second); A.push_back(pdg::IonPdgCodeToA(TgtPdgCode)); Z.push_back(pdg::IonPdgCodeToZ(TgtPdgCode)); Mole.push_back(TgtMassFrac/mass); } else { LOG("gSeaGen", pFATAL) << "Element "<< TgtPdgCode<<" not found in file " << filename <<", exiting... "; exit(1); } } MoleMin=*min_element(Mole.begin(),Mole.end()); fGenPar->MolarMassSR=0.; fGenPar->NNucleonsSR=0; fGenPar->NElectronsSR=0; for(unsigned iComp=0; iCompMolarMassSR+=(double)nAtoms*AtomMass[iComp]; fGenPar->NNucleonsSR+=nAtoms*A[iComp]; fGenPar->NElectronsSR+=nAtoms*Z[iComp]; AtomInMolecule=MassFrac[iComp]/AtomMass[iComp]/MoliFracMass; #ifdef _PROPOSAL_ENABLED__ if(fGenPar->PropCode.compare("PROPOSAL")==0)ComponentSR.push_back(PROPOSAL::Components::Component(AtomName[iComp].c_str(),Z[iComp],AtomMass[iComp],AtomInMolecule)); #endif #ifdef _JPP_ENABLED__ if(fGenPar->PropCode.compare("JPP")==0){ JPHYSICS::JRadiation jppelement(Z[iComp],A[iComp],40,0.01,0.1,0.1); // (Z,A,integration steps, minimal energyfor brems, minimal energy for epair, minimal energy for nuclear interaction) JLANG::JSharedPointer jppEl(new JPHYSICS::JRadiationFunction(jppelement, 300, 0.2, 1.0e11)); //(element, number of energy bins, min energy, max energy) radiationSR.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSR * MassFrac[iComp], EErad)); radiationSR.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSR * MassFrac[iComp], Brems)); radiationSR.push_back( new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSR * MassFrac[iComp], GNrad)); ionizationSR.push_back(new JPHYSICS::JRadiationSource(jppEl, fGenPar->RhoSR * MassFrac[iComp], Ion)); } #endif } // Mantle vector().swap(AtomMass); vector().swap(Mole); vector().swap(MassFrac); vector().swap(Z); vector().swap(A); vector().swap(AtomName); iter = fGenPar->MantleComp.begin(); for ( ; iter != fGenPar->MantleComp.end(); ++iter) { TgtPdgCode=iter->first; TgtMassFrac=iter->second; if(AtomicMass.find(TgtPdgCode) != AtomicMass.end()){ if(TgtMassFrac>TgtMassFracMax){ TgtMassFracMax=TgtMassFrac; MoliFracMass=TgtMassFrac/AtomicMass.find(TgtPdgCode)->second; } MassFrac.push_back(TgtMassFrac); mass=AtomicMass.find(TgtPdgCode)->second; AtomMass.push_back(AtomicMass.find(TgtPdgCode)->second); AtomName.push_back(ElementName.find(TgtPdgCode)->second); A.push_back(pdg::IonPdgCodeToA(TgtPdgCode)); Z.push_back(pdg::IonPdgCodeToZ(TgtPdgCode)); Mole.push_back(TgtMassFrac/mass); } else { LOG("gSeaGen", pFATAL) << "Element "<< TgtPdgCode<<" not found in file " << filename <<", exiting... "; exit(1); } } MoleMin=*min_element(Mole.begin(),Mole.end()); fGenPar->MolarMassMantle=0.; fGenPar->NNucleonsMantle=0; fGenPar->NElectronsMantle=0; for(unsigned iComp=0; iCompMolarMassMantle+=(double)nAtoms*AtomMass[iComp]; fGenPar->NNucleonsMantle+=nAtoms*A[iComp]; fGenPar->NElectronsMantle+=nAtoms*Z[iComp]; } // Core vector().swap(AtomMass); vector().swap(Mole); vector().swap(MassFrac); vector().swap(Z); vector().swap(A); vector().swap(AtomName); iter = fGenPar->CoreComp.begin(); for ( ; iter != fGenPar->CoreComp.end(); ++iter) { TgtPdgCode=iter->first; TgtMassFrac=iter->second; if(AtomicMass.find(TgtPdgCode) != AtomicMass.end()){ if(TgtMassFrac>TgtMassFracMax){ TgtMassFracMax=TgtMassFrac; MoliFracMass=TgtMassFrac/AtomicMass.find(TgtPdgCode)->second; } MassFrac.push_back(TgtMassFrac); mass=AtomicMass.find(TgtPdgCode)->second; AtomMass.push_back(AtomicMass.find(TgtPdgCode)->second); AtomName.push_back(ElementName.find(TgtPdgCode)->second); A.push_back(pdg::IonPdgCodeToA(TgtPdgCode)); Z.push_back(pdg::IonPdgCodeToZ(TgtPdgCode)); Mole.push_back(TgtMassFrac/mass); } else { LOG("gSeaGen", pFATAL) << "Element "<< TgtPdgCode<<" not found in file " << filename <<", exiting... "; exit(1); } } MoleMin=*min_element(Mole.begin(),Mole.end()); fGenPar->MolarMassCore=0.; fGenPar->NNucleonsCore=0; fGenPar->NElectronsCore=0; for(unsigned iComp=0; iCompMolarMassCore+=(double)nAtoms*AtomMass[iComp]; fGenPar->NNucleonsCore+=nAtoms*A[iComp]; fGenPar->NElectronsCore+=nAtoms*Z[iComp]; } vector().swap(AtomMass); vector().swap(Mole); vector().swap(MassFrac); vector().swap(Z); vector().swap(A); vector().swap(AtomName); return; } //___________________________________________________________________________ void GGenerateEvent::ConfigProposal(void){ // configure Proposal if enabled at compiling and choosen at run-times #ifdef _PROPOSAL_ENABLED__ double global_ecut_inside = 500; double global_ecut_infront = -1; double global_ecut_behind = -1; double global_vcut_inside = -1; double global_vcut_infront = 0.05; //0.001; double global_vcut_behind = -1; bool global_cont_inside = false; bool global_cont_infront = true; bool global_cont_behind = false; bool do_interpolation = true; string path_to_tables_readonly=PROPOSALRES; string path_to_tables=PROPOSALRES; if(!gSystem->Getenv("PROPOSALCONF")){ struct stat info; if( stat( path_to_tables.c_str(), &info ) != 0 ){ LOG("GGenerateEvent", pFATAL) <<"Error: invalid PROPOSAL table path "<RanSeed); PROPOSAL::InterpolationDef interpolation_def; interpolation_def.path_to_tables=path_to_tables; interpolation_def.path_to_tables_readonly=path_to_tables_readonly; unique_ptr sec_def_global(new PROPOSAL::Sector::Definition()); if(gSystem->Getenv("PROPOSALCONF")){ string config_file = gSystem->Getenv("PROPOSALCONF"); if(!(access(config_file.c_str(), F_OK) != -1)){ LOG("GGenerateEvent", pFATAL) << config_file <<" does not exist --> check env. parameter PROPOSALCONF; exiting..."; exit(1); } else LOG("GGenerateEvent", pINFO) <<"Reading PROPOSAL inputs from file "<< config_file; // check if table paths exist before proposal parses the file to avoid abort struct stat info; string pathname; boost::property_tree::ptree json_config_check; boost::property_tree::json_parser::read_json(config_file, json_config_check); vector option; option.push_back("global.interpolation.path_to_tables"); option.push_back("global.interpolation.path_to_tables_readonly"); boost::optional optional_param; for(unsigned i=0; i(option[i]); if(optional_param){ pathname = optional_param.get(); if( stat( pathname.c_str(), &info ) != 0 ){ LOG("GGenerateEvent", pFATAL) <<"Error in PROPOSAL input file "<< config_file<<": invalid path "<> json_config; nlohmann::json json_global; if(json_config.contains("global")){ json_global = json_config["global"]; if (json_global.contains("interpolation")){ nlohmann::json json_interpol = json_global["interpolation"]; interpolation_def = PROPOSAL::InterpolationDef(json_interpol); do_interpolation = json_interpol.value("do_interpolation", do_interpolation); } sec_def_global.reset(new PROPOSAL::Sector::Definition(json_global)); } } shared_ptr med_air = PROPOSAL::CreateMedium("Air", 1.); shared_ptr med_sw = make_shared("gSeaWater", 1, 75.0, // I -3.5017, // C 0.09116, // a 3.4773, // m 0.2400, // X0 2.8004, // X1 0, // d0 fGenPar->RhoSW, // massDensity ComponentSW); shared_ptr med_sr = make_shared("gRock", 1, 149.0, // I -5.053, // C 0.078, // a 3.645, // m 0.288, // X0 3.196, // X1 0, // d0 fGenPar->RhoSR, // massDensity ComponentSR); shared_ptr geoCan; geoCan = make_shared(PROPOSAL::Vector3D(0, 0, (fGenPar->Can1+fGenPar->Can2)*0.5), fGenPar->Can3, 0., (fGenPar->Can2-fGenPar->Can1)); shared_ptr geometrySW; geometrySW = make_shared(PROPOSAL::Vector3D(0, 0, (-fGenPar->SeaBottomRadius)), (fGenPar->SeaBottomRadius+fGenPar->SiteDepth), fGenPar->SeaBottomRadius); geometrySW->SetHierarchy(1); shared_ptr geoSeaWater=geometrySW; shared_ptr geometrySR; geometrySR = make_shared(PROPOSAL::Vector3D(0, 0, (-fGenPar->SeaBottomRadius)), fGenPar->SeaBottomRadius, 0.); geometrySR->SetHierarchy(1); shared_ptr geoRock=geometrySR; shared_ptr geometryAir; geometryAir = make_shared(PROPOSAL::Vector3D(0, 0, (-fGenPar->SeaBottomRadius)), 1000000000.,(fGenPar->SeaBottomRadius+fGenPar->SiteDepth)); geometryAir->SetHierarchy(1); shared_ptr geoAir=geometryAir; vector sectors; PROPOSAL::Sector::Definition sec_def_sector=*sec_def_global; sec_def_sector.location = PROPOSAL::Sector::ParticleLocation::InfrontDetector; sec_def_sector.cut_settings.SetEcut(global_ecut_infront); sec_def_sector.cut_settings.SetVcut(global_vcut_infront); sec_def_sector.do_continuous_randomization = global_cont_infront; sec_def_sector.SetMedium(med_air); sec_def_sector.SetGeometry(geoAir); sectors.push_back(sec_def_sector); sec_def_sector.SetMedium(med_sw); sec_def_sector.SetGeometry(geoSeaWater); sectors.push_back(sec_def_sector); sec_def_sector.SetMedium(med_sr); sec_def_sector.SetGeometry(geoRock); sectors.push_back(sec_def_sector); sec_def_sector.location = PROPOSAL::Sector::ParticleLocation::InsideDetector; sec_def_sector.cut_settings.SetEcut(global_ecut_inside); sec_def_sector.cut_settings.SetVcut(global_vcut_inside); sec_def_sector.do_continuous_randomization = global_cont_inside; sec_def_sector.SetMedium(med_air); sec_def_sector.SetGeometry(geoAir); sectors.push_back(sec_def_sector); sec_def_sector.SetMedium(med_sw); sec_def_sector.SetGeometry(geoSeaWater); sectors.push_back(sec_def_sector); sec_def_sector.SetMedium(med_sr); sec_def_sector.SetGeometry(geoRock); sectors.push_back(sec_def_sector); sec_def_sector.location = PROPOSAL::Sector::ParticleLocation::BehindDetector; sec_def_sector.cut_settings.SetEcut(global_ecut_behind); sec_def_sector.cut_settings.SetVcut(global_vcut_behind); sec_def_sector.do_continuous_randomization = global_cont_behind; sec_def_sector.SetMedium(med_air); sec_def_sector.SetGeometry(geoAir); sectors.push_back(sec_def_sector); sec_def_sector.SetMedium(med_sw); sec_def_sector.SetGeometry(geoSeaWater); sectors.push_back(sec_def_sector); sec_def_sector.SetMedium(med_sr); sec_def_sector.SetGeometry(geoRock); sectors.push_back(sec_def_sector); if(do_interpolation) ProposalMu = new PROPOSAL::Propagator(PROPOSAL::MuMinusDef::Get(), sectors, geoCan,interpolation_def); else ProposalMu = new PROPOSAL::Propagator(PROPOSAL::MuMinusDef::Get(), sectors, geoCan); #endif return; } //________________________________________________________________________________________ GSeaRealAtmoFlux * GGenerateEvent::GetFluxDiffuse(void) { // Instantiate the diffuse flux driver GSeaRealAtmoFlux * atmo_flux_driver = new GSeaRealAtmoFlux(fGenPar); // set rotation for coordinate tranformation from the topocentric horizontal // system to a user-defined coordinate system: if(!fGenPar->Rotation.IsIdentity()) { atmo_flux_driver->SetUserCoordSystem(fGenPar->Rotation); } fGenPar->GlobalGenWeight=atmo_flux_driver->GetGlobalGenWeight(); return atmo_flux_driver; } //________________________________________________________________________________________ GSeaRealPointFlux* GGenerateEvent::GetFluxPoint(void) { // Instantiate the point flux driver GSeaRealPointFlux * point_flux_driver = new GSeaRealPointFlux(fGenPar); // set rotation for coordinate tranformation from the topocentric horizontal // system to a user-defined coordinate system: if(!fGenPar->Rotation.IsIdentity()) { point_flux_driver->SetUserCoordSystem(fGenPar->Rotation); } fGenPar->GlobalGenWeight=point_flux_driver->GetGlobalGenWeight(); return point_flux_driver; } //________________________________________________________________________________________ int GGenerateEvent::PosInCan(double VertX, double VertY, double VertZ){ // check if the interaction vertex is inside the detector can int VerInCan = 0; if (sqrt(VertX*VertX+VertY*VertY)<=fGenPar->Can3) { if(VertZ>= fGenPar->Can1 && VertZ<=fGenPar->Can2)VerInCan = 1; else if(VertZ>= -fGenPar->HBedRock && VertZ<0)VerInCan = 2; } return VerInCan; } //________________________________________________________________________________________ void GGenerateEvent::FillEvent() { if ( fGenPar->GenCode.compare("CORSIKA")==0 ) { // CORSIKA fSeaEvent->Tracks.reserve( fFileDriver->OtherTracks.size() + fFileDriver->SecondaryTracks.size() ); // preallocate memory fSeaEvent->Tracks.insert( fSeaEvent->Tracks.end(), fFileDriver->OtherTracks.begin(), fFileDriver->OtherTracks.end() ); // first pass OtherTracks vector().swap(fFileDriver->OtherTracks); // we no longer need those fSeaEvent->Tracks.insert( fSeaEvent->Tracks.end(), fFileDriver->SecondaryTracks.begin(), fFileDriver->SecondaryTracks.end() ); // then pass SecondaryTracks vector().swap(fFileDriver->SecondaryTracks); // we no longer need those } else { fSeaEvent->Tracks = gWrittenTracks; vector().swap(gWrittenTracks); // clear memory } fSeaEvent->NTracks=fSeaEvent->Tracks.size(); } //________________________________________________________________________________________ bool GGenerateEvent::PropagationMu(GSeaTrack *t) { // wrapper function to call propagation with regular GSeaTrack PropaTrack pt; copy( *t, pt ); bool r = PropagationMu ( &pt ); copy( pt, *t ); return r; } //________________________________________________________________________________________ int GGenerateEvent::DistToCan(GSeaTrack *t) { // wrapper function to call DistToCan with regular GSeaTrack PropaTrack pt; copy( *t, pt ); int r = DistToCan ( &pt ); copy( pt, *t ); return r; } //________________________________________________________________________________________ bool GGenerateEvent::PropagationMu(PropaTrack *trackf){ if (trackf->E<0.16) return false; // muon at rest (no propagation) PropaTrack track=*trackf; int iFlagCan=-1; bool LepCan=false; double Rho2Upp,Rho2Low; double REarthHcan = fGenPar->SeaBottomRadius+fGenPar->Can2; double REarthHcan2 = REarthHcan*REarthHcan; double DepthW,DepthR; // depths left to travel through water and rock double range=(pow(10.,fGenPar->RangeSW.Eval(log10(track.E)))/(fGenPar->RhoSW)) + fGenPar->MuonRangeTolerance; #ifdef _MUSIC_ENABLED__ int MUMODEL2=1; double DepthsF; double RhoSW=fGenPar->RhoSW; double RhoSR=fGenPar->RhoSR; #endif iFlagCan=DistToCan(&track); // check if the muon track intercepts the can double Rho = sqrt(track.V1*track.V1 + track.V2*track.V2); double HMargin = fGenPar->PropagationTolerance*sqrt((1.-track.D3)*(1.+track.D3)); // (1.-track.D3)*(1.+track.D3) is more efficient than (1.-track.D3*track.D3) Rho2Upp = (Rho + HMargin)*(Rho + HMargin); Rho2Low = (Rho - HMargin)*(Rho - HMargin); if (fGenPar->GenCode.compare("CORSIKA")==0){ // check whether the can top level was reached (only needed for CORSIKA) if ( track.V3 > sqrt(REarthHcan2-Rho2Upp) - fGenPar->SeaBottomRadius + fGenPar->PropagationTolerance*fabs(track.D3) ) AboveCan=true; // track reached the can top height (within Tolerance) else AboveCan=false; } if (iFlagCan==1) LepCan = true; // the track is inside the can (no propagation) else if (iFlagCan==2 || iFlagCan==-1) LepCan = false; // the track does not intercept the can (no propagation) else if (range < (fDistaCan-fGenPar->PropagationTolerance) ) LepCan = false; // muon does not have chance to reach the can (no propagation) else { #ifdef _JPP_ENABLED__ double JPP_Emin=1; //GeV Minimal energy, if energy is below 1 GeV, only ionization is computed until stop. #endif while( !LepCan ){ // propagate until we make it to the can (or fail to) #ifdef _JPP_ENABLED__ JSIRENE::JVertex vertex(0, track.T, track.E); //The "Z" to propagate in jpp is the first component. #endif bool InRock=false; DepthR=fabs(-track.V3/track.D3); if(track.V3<=0){ if (DepthR>fGenPar->PropagationTolerance) InRock=true; else InRock=false; } if(InRock){ // Transport in rock if(DepthRDistAfterFirstStep)DepthR=min(DepthR,fGenPar->PropagationStep); else DepthR=max(DepthR-fGenPar->DistAfterFirstStep,fGenPar->PropagationTolerance); // Transport the muon in rock if(fGenPar->PropCode.compare("PropaMuon")==0){ if(!PropaMuonSR->Propagate(track,DepthR)) return false; *trackf = PropaMuonSR->GetPropaTrack(); } #ifdef _MUSIC_ENABLED__ else if(fGenPar->PropCode.compare("MUSIC")==0){ DepthR*=100.; track.V1*=100.; track.V2*=100.; track.V3*=100.; muon_transport_sr_(&track.V1,&track.V2,&track.V3,&track.D1,&track.D2,&track.D3,&track.E,&DepthR,&track.T,&RhoSR,&MUMODEL2,&MUMODEL2,&trackf->V1,&trackf->V2,&trackf->V3,&trackf->D1,&trackf->D2,&trackf->D3,&trackf->E,&DepthsF,&trackf->T); trackf->V1/=100.; trackf->V2/=100.; trackf->V3/=100.; } #endif #ifdef _JPP_ENABLED__ else if(fGenPar->PropCode.compare("JPP")==0){ while (vertex.getE() >= JPP_Emin && vertex.getZ() < DepthR) { const int Nr = radiationSR.size(); double li[Nr]; // inverse interaction lengths double ls = 1.0e-5; // minimal total inverse interaction length (100 km)^-1 for (int i = 0; i != Nr; ++i) { ls += li[i] = radiationSR[i]->getInverseInteractionLength(vertex.getE()); } const int Ni = ionizationSR.size(); double Ai[Ni]; double As=0; for (int i = 0; i != Ni; ++i) { As +=Ai[i]= ionizationSR[i]->getA(vertex.getE()); } const double ds = min(DepthR-vertex.getZ(),min(gRandom->Exp(1.0) / ls, vertex.getRange(As))); vertex.step(ds,As); //Ionization continuous losses if (vertex.getE() >= JPP_Emin) { double Es = 0.0; // shower energy [GeV] double y = gRandom->Uniform(ls); for (int i = 0; i != Nr; ++i) { y -= li[i]; if (y < 0.0) { Es = radiationSR[i]->getEnergyOfShower(vertex.getE()); // shower energy [GeV] break; } } vertex.applyEloss(Es); } } // add muon end point const int Ni = ionizationSR.size(); double Ai[Ni]; double As=0; for (int i = 0; i != Ni; ++i) { As +=Ai[i]= ionizationSR[i]->getA(vertex.getE()); } if (vertex.getE() < JPP_Emin && vertex.getRange(As) > 0.0 && vertex.getZ() < DepthR) { vertex.step(min(vertex.getRange(),DepthR-vertex.getZ()),As); } trackf->E=vertex.getE(); //Projecting the distance traveled into the track direction. trackf->V1+=vertex.getZ()*track.D1; trackf->V2+=vertex.getZ()*track.D2; trackf->V3+=vertex.getZ()*track.D3; trackf->T=vertex.getT(); } #endif #ifdef _PROPOSAL_ENABLED__ else if(fGenPar->PropCode.compare("PROPOSAL")==0){ PROPOSAL::Vector3D direction(track.D1,track.D2,track.D3); direction.CalculateSphericalCoordinates(); PROPOSAL::DynamicData particle_in(PROPOSAL::MuMinusDef::Get().particle_type); particle_in.SetEnergy(track.E*1E3); // [MeV] particle_in.SetPropagatedDistance(0); particle_in.SetPosition(PROPOSAL::Vector3D(track.V1*1E2,track.V2*1E2,track.V3*1E2)); particle_in.SetDirection(direction); particle_in.SetTime(track.T*1e-09); // in sec PROPOSAL::Secondaries sec = ProposalMu->Propagate(particle_in,DepthR*100.); PROPOSAL::DynamicData particle_out = sec.GetModifyableSecondaries().back(); if (particle_out.GetPropagatedDistance()>0) trackf->E = particle_out.GetEnergy()/1.E3; else trackf->E=0.; // muon decayed trackf->V1=particle_out.GetPosition().GetX()/1.E2; trackf->V2=particle_out.GetPosition().GetY()/1.E2; trackf->V3=particle_out.GetPosition().GetZ()/1.E2; trackf->D1=particle_out.GetDirection().GetX(); trackf->D2=particle_out.GetDirection().GetY(); trackf->D3=particle_out.GetDirection().GetZ(); trackf->T =particle_out.GetTime()/1e-09; } #endif trackf->DirectionCosines(); } else{ // Transport in sea water if (AboveCan && fGenPar->GenCode.compare("CORSIKA")==0 ) { // propagation in 2 steps for CORSIKA showers // First we propagate from sea to the level of the top cap of the can. // Here we need to use eq. for intersection of a line and a sphere (R=R_Earth+h_can) // the equations can be found e.g. here: http://paulbourke.net/geometry/circlesphere/ // as the two points on the line we pick x1 = track.V1 and x2 = track.V1+1*track.D1 (and the same for V2 and V3) // This simplifies a to: a = track.D1*track.D1 + track.D2*track.D2 + track.D3*track.D3 = 1 double z = fGenPar->SeaBottomRadius+track.V3; double b = 2.*( track.D1*track.V1 + track.D2*track.V2 + track.D3*z ); double rho2 = track.V1*track.V1 + track.V2*track.V2; double c = rho2 + z*z - REarthHcan2; DepthW = 0.5*(-b-sqrt(b*b-4.*c)); // equations simplify and t= DepthW. We choose the solution // closer to the track position, since the other intersection point // would be on the other side of the Earth // We first only travel to the can top level. } else DepthW=fDistaCan; // for GENIE we directly propagate all the way to the can // if ( fGenPar->PropCode.compare("PROPOSAL") !=0 ) { // not necessary for PROPOSAL: if( DepthW < fGenPar->DistAfterFirstStep ) DepthW=min(DepthW,fGenPar->PropagationStep); else DepthW=max(DepthW-fGenPar->DistAfterFirstStep,fGenPar->PropagationTolerance); // } if( fGenPar->PropCode.compare("PropaMuon")==0 ){ if( ! PropaMuonSW->Propagate(track,DepthW) ){LepCan=false; break;} *trackf = PropaMuonSW->GetPropaTrack(); } #ifdef _MUSIC_ENABLED__ else if(fGenPar->PropCode.compare("MUSIC")==0){ DepthW*=100.; track.V1*=100.; track.V2*=100.; track.V3*=100.; muon_transport_sw_(&track.V1,&track.V2,&track.V3,&track.D1,&track.D2,&track.D3,&track.E,&DepthW,&track.T,&RhoSW,&MUMODEL2,&MUMODEL2,&trackf->V1,&trackf->V2,&trackf->V3,&trackf->D1,&trackf->D2,&trackf->D3,&trackf->E,&DepthsF,&trackf->T); trackf->V1/=100.; trackf->V2/=100.; trackf->V3/=100.; } #endif #ifdef _JPP_ENABLED__ else if(fGenPar->PropCode.compare("JPP")==0){ while (vertex.getE() >= JPP_Emin && vertex.getZ() < DepthW) { const int Nr = radiationSW.size(); double li[Nr]; // inverse interaction lengths double ls = 1.0e-5; // minimal total inverse interaction length (100 km)^-1 for (int i = 0; i != Nr; ++i) { ls += li[i] = radiationSW[i]->getInverseInteractionLength(vertex.getE()); } const int Ni = ionizationSW.size(); double Ai[Ni]; double As=0; for (int i = 0; i != Ni; ++i) { As +=Ai[i]= ionizationSW[i]->getA(vertex.getE()); } const double ds = min(DepthW-vertex.getZ(),min(gRandom->Exp(1.0) / ls, vertex.getRange(As))); vertex.step(ds,As); //Ionization continuous losses if (vertex.getE() >= JPP_Emin) { double Es = 0.0; // shower energy [GeV] double y = gRandom->Uniform(ls); for (int i = 0; i != Nr; ++i) { y -= li[i]; if (y < 0.0) { Es = radiationSW[i]->getEnergyOfShower(vertex.getE()); // shower energy [GeV] break; } } vertex.applyEloss(Es); } } // add muon end point const int Ni = ionizationSW.size(); double Ai[Ni]; double As=0; for (int i = 0; i != Ni; ++i) { As +=Ai[i]= ionizationSW[i]->getA(vertex.getE()); } if (vertex.getE() < JPP_Emin && vertex.getRange(As) > 0.0 && vertex.getZ()< DepthW) { vertex.step(min(DepthW-vertex.getZ(),vertex.getRange()),As); } trackf->E=vertex.getE(); //Projecting the distance traveled into the track direction. trackf->V1+=vertex.getZ()*track.D1; trackf->V2+=vertex.getZ()*track.D2; trackf->V3+=vertex.getZ()*track.D3; trackf->T=vertex.getT(); } #endif #ifdef _PROPOSAL_ENABLED__ else if(fGenPar->PropCode.compare("PROPOSAL")==0){ PROPOSAL::Vector3D direction(track.D1,track.D2,track.D3); direction.CalculateSphericalCoordinates(); PROPOSAL::DynamicData particle_in(PROPOSAL::MuMinusDef::Get().particle_type); particle_in.SetEnergy(track.E*1E3); // [MeV] particle_in.SetPropagatedDistance(0); particle_in.SetPosition(PROPOSAL::Vector3D(track.V1*1E2,track.V2*1E2,track.V3*1E2)); particle_in.SetDirection(direction); particle_in.SetTime(track.T*1e-09); // in sec PROPOSAL::Secondaries sec = ProposalMu->Propagate(particle_in,DepthW*100.); PROPOSAL::DynamicData particle_out = sec.GetModifyableSecondaries().back(); if (particle_out.GetPropagatedDistance()>0) trackf->E =particle_out.GetEnergy()/1.E3; else trackf->E=0.; // muon decayed trackf->V1=particle_out.GetPosition().GetX()/1.E2; trackf->V2=particle_out.GetPosition().GetY()/1.E2; trackf->V3=particle_out.GetPosition().GetZ()/1.E2; trackf->D1=particle_out.GetDirection().GetX(); trackf->D2=particle_out.GetDirection().GetY(); trackf->D3=particle_out.GetDirection().GetZ(); trackf->T =particle_out.GetTime()/1e-09; } #endif trackf->DirectionCosines(); } if (trackf->E<0.16) { track=*trackf; LepCan=false; break;} // muon at rest (no propagation) else if ((trackf->D3>0.0 || trackf->V3>track.V3) && (trackf->V3>fGenPar->Can2)) { // muon scattered and it now upgoing (and above can)! track=*trackf; LepCan=false; break; } else{ track=*trackf; // update the track iFlagCan = DistToCan(&track); if (AboveCan && fGenPar->GenCode.compare("CORSIKA")==0 ) { double Rho = sqrt(track.V1*track.V1 + track.V2*track.V2); Rho2Upp = (Rho + HMargin)*(Rho + HMargin); Rho2Low = (Rho - HMargin)*(Rho - HMargin); if ( ( track.V3 <= sqrt(REarthHcan2-Rho2Upp) - fGenPar->SeaBottomRadius + fGenPar->PropagationTolerance*fabs(track.D3) ) & ( track.V3 > sqrt(REarthHcan2-Rho2Low) - fGenPar->SeaBottomRadius - fGenPar->PropagationTolerance*fabs(track.D3) ) ) { // track reached the can top height (within Tolerance) ; we use an approximation here, locally neglecting the curvature AboveCan=false; if (iFlagCan==-1){ LepCan=false; } else if (iFlagCan==0) { LepCan=false; NeedToPropagate=true; } else if (iFlagCan==1) { LepCan=true; } // the track is inside the can (no propagation) else if (iFlagCan==2) { LepCan=false; NeedToShift=true; } // the track missed the can (needs shifting) break; } else continue; // did not reach the right depth yet } else { // track is propagated below the can top height if (iFlagCan==-1){ LepCan=false; break; } else if (iFlagCan==0) { LepCan=false; continue; } else if (iFlagCan==1) { LepCan=true; break; } // the track is inside the can (no propagation) else if (iFlagCan==2) { LepCan=false; break; } // the track missed the can (needs shifting) } } } // close while(1) } if (LepCan) trackf->SetStatus(1); return LepCan; } //________________________________________________________________________________________ #ifdef _HETAU_ENABLED__ double GGenerateEvent::PropagationTau(GSeaTrack *t) { // wrapper function to call propagation with regular GSeaTrack PropaTrack pt; copy( *t, pt ); double r = PropagationTau ( &pt ); copy( pt, *t ); return r; } double GGenerateEvent::PropagationTau(PropaTrack *trackf){ if(trackf->E<1.8) return -1; // tau at rest (no propagation) int iFlagCan = DistToCan(trackf); // We want to know if the tau track intercepts the can RandomGen * rnd = RandomGen::Instance(); if (iFlagCan==2 || iFlagCan==-1) return -1; // the track does not intercept the can (no propagation) else { //the track is outside the can and pointing double tauti = -log( rnd->RndGen().Rndm() ) * lifetime*1e9; //ns for tausic PropaTrack track=*trackf; double RhoSW=fGenPar->RhoSW; double RhoSR=fGenPar->RhoSR; double DepthW,DepthR; double DepthsF; int idec; double tautf = -1.; bool InRock=false; while(1){ DepthR=fabs(-track.V3/track.D3); if(track.V3<=0.) InRock = (DepthR>fGenPar->PropagationTolerance*fabs(track.D3)) ? true : false; if (tautf>0.) tauti = tautf; if(InRock){ // I'm in rock DepthR*=100.; track.V1*=100.; track.V2*=100.; track.V3*=100.; tau_transport_sr_(&track.V1,&track.V2,&track.V3,&track.D1,&track.D2,&track.D3,&track.E,&DepthR,&track.T,&RhoSR,&TAUMODEL,&TAUMODEL,&trackf->V1,&trackf->V2,&trackf->V3,&trackf->D1,&trackf->D2,&trackf->D3,&trackf->E,&DepthsF,&trackf->T,&idec,&ITFLAG,&tauti,&tautf); trackf->V1/=100.; trackf->V2/=100.; trackf->V3/=100.; } else{ DepthW=fDistaCan*100.; track.V1*=100.; track.V2*=100.; track.V3*=100.; tau_transport_sw_(&track.V1,&track.V2,&track.V3,&track.D1,&track.D2,&track.D3,&track.E,&DepthW,&track.T,&RhoSW,&TAUMODEL,&TAUMODEL,&trackf->V1,&trackf->V2,&trackf->V3,&trackf->D1,&trackf->D2,&trackf->D3,&trackf->E,&DepthsF,&trackf->T,&idec,&ITFLAG,&tauti,&tautf); trackf->V1/=100.; trackf->V2/=100.; trackf->V3/=100.; } track=*trackf; iFlagCan = DistToCan(&track); if (track.E<1.8) return 0.; //tau at rest before reaching the can if (idec==1) return 0.; //tau has decayed before reaching the can if (fDistaCan>0 && fDistaCan<=fGenPar->PropagationTolerance) { trackf->SetStatus(1); return tautf; } //tau has reached the can before decaying if (iFlagCan==1) { trackf->SetStatus(1); return tautf; } //tau has reached the can before decaying if (iFlagCan==2 || iFlagCan==-1) return -1; //tau not pointing at detector } // close while(1) } return -1; } //________________________________________________________________________________________ vector GGenerateEvent::DecayTau(GSeaTrack tau, double tauti) { vector decayed; GSeaTrack propagated_tau = tau; if (tauti>0) { double RhoSW=fGenPar->RhoSW; double Depth=1e10; double DepthsF; int idec; double tautf = -1.; tau.V1*=100.; tau.V2*=100.; tau.V3*=100.; tau_transport_sw_(&tau.V1,&tau.V2,&tau.V3,&tau.D1,&tau.D2,&tau.D3,&tau.E,&Depth,&tau.T,&RhoSW,&TAUMODEL,&TAUMODEL,&propagated_tau.V1,&propagated_tau.V2,&propagated_tau.V3,&propagated_tau.D1,&propagated_tau.D2,&propagated_tau.D3,&propagated_tau.E,&DepthsF,&propagated_tau.T,&idec,&ITFLAG,&tauti,&tautf); propagated_tau.V1/=100.; propagated_tau.V2/=100.; propagated_tau.V3/=100.; } //now we decay the propagated tau double momf = sqrt( propagated_tau.E*propagated_tau.E - kTauMass2 ); TauolaHEPEVTEvent * Tauola_evt = new TauolaHEPEVTEvent(); TauolaHEPEVTParticle *Tauola_tau = new TauolaHEPEVTParticle( propagated_tau.PdgCode, 1, propagated_tau.D1*momf, propagated_tau.D2*momf, propagated_tau.D3*momf, propagated_tau.E, kTauMass, -1, -1, -1, -1 ); Tauola_evt->addParticle(Tauola_tau); Tauola::decayOne(Tauola_tau,true,0.,0.,1); for ( int sec=1; secgetParticleCount(); sec++ ) { GSeaTrack DecTrack; DecTrack.Status = Tauola_evt->getParticle(sec)->getStatus(); DecTrack.PdgCode = Tauola_evt->getParticle(sec)->getPdgID(); DecTrack.E = Tauola_evt->getParticle(sec)->getE(); double px = Tauola_evt->getParticle(sec)->getPx(); double py = Tauola_evt->getParticle(sec)->getPy(); double pz = Tauola_evt->getParticle(sec)->getPz(); double mom = sqrt(px*px + py*py + pz*pz); DecTrack.D1 = px/mom; DecTrack.D2 = py/mom; DecTrack.D3 = pz/mom; DecTrack.V1 = propagated_tau.V1; DecTrack.V2 = propagated_tau.V2; DecTrack.V3 = propagated_tau.V3; DecTrack.T = propagated_tau.T; DecTrack.MotherId = propagated_tau.Id; DecTrack.Length = 0.; if (this->PosInCan(DecTrack.V1,DecTrack.V2,DecTrack.V3)==0) DecTrack.Status *= -1; decayed.push_back(DecTrack); } return decayed; } #endif //___________________________________________________________________________ int GGenerateEvent::DistToCan(PropaTrack *track) { // iFlag = -1 something has gone wrong! // iFlag = 0 if the track intercepts the can // iFlag = 1 if the point is inside the can // iFlag = 2 if the track never intercepts the can int iFlag=-1; fDistaCan = -1; // double Zmin= fGenPar->Can1; // can base // double Zmax= fGenPar->Can2; // can cap // double Radius = fGenPar->Can3; // can radius double utit=sqrt(track->D1*track->D1+track->D2*track->D2+track->D3*track->D3); track->D1/=utit; track->D2/=utit; track->D3/=utit; if(track->V3<=fGenPar->Can2 && track->V3>=fGenPar->Can1 && track->V1*track->V1+track->V2*track->V2<=fGenPar->Can3*fGenPar->Can3){ // inside the can iFlag = 1; fDistaCan = 0.; } else if(fabs(track->D3)==1){ // vertical muon if (sqrt(track->V1*track->V1+track->V2*track->V2) > fGenPar->Can3) iFlag = 2; // off in xy direction else if(track->D3==1 && track->V3>fGenPar->Can2) iFlag = 2; else if(track->D3==-1 && track->V3Can1) iFlag = 2; else{ iFlag = 0; if(track->D3==1) fDistaCan = fabs(fGenPar->Can1-track->V3); // up-going else fDistaCan = fabs(track->V3-fGenPar->Can2); // down-going } } else { // inclined muon double a = track->D1*track->D1+track->D2*track->D2; double b = 2.*(track->V1*track->D1+track->V2*track->D2); double c = track->V1*track->V1+track->V2*track->V2-fGenPar->Can3*fGenPar->Can3; double delta=b*b-4.*a*c; if (delta<0) iFlag = 2; else { delta = sqrt(delta); double t1=(-b-delta)/(2.*a); double t2=(-b+delta)/(2.*a); double z1=track->V3+track->D3*t1; double z2=track->V3+track->D3*t2; if((z1Can1 && z2Can1) || (z1>fGenPar->Can2 && z2>fGenPar->Can2)) { iFlag = 2; } else if(z1>=fGenPar->Can1 && z1<=fGenPar->Can2){ if(t1>=0) { iFlag = 0; fDistaCan = t1; } else iFlag = 2; // intersection before the vertex } else if(z1Can1 && track->D3>0){ // up-going iFlag = 0; t1=(fGenPar->Can1-track->V3)/track->D3; if(t1>=0) fDistaCan=t1; else iFlag = 2; // intersection before the vertex } else if(z1>fGenPar->Can2 && track->D3<0){ // down-going iFlag = 0; t1=(fGenPar->Can2-track->V3)/track->D3; if(t1>=0) fDistaCan=t1; else iFlag = 2; // intersection before the vertex } } } if (iFlag==0 && fDistaCan>=0 && fDistaCan<=fGenPar->PropagationTolerance) iFlag=1; // the muon is close enough return iFlag; } //___________________________________________________________________________ void GGenerateEvent::RotateAroundEarth(GSeaTrack *Track, const array,3> &RotMatrix) { // This function rotates position (V1,V2,V3) and direction (D1,D2,D3) // of a given Track around the center of the Earth // by a specified RotMatrix Track->V3 += fGenPar->SeaBottomRadius; // we switch to the coordinate system centered at center of the Earth // position: array newVec = ApplyRotationMatrix({Track->V1,Track->V2,Track->V3}, RotMatrix); Track->V1 = newVec[0]; Track->V2 = newVec[1]; Track->V3 = newVec[2]; Track->V3 -= fGenPar->SeaBottomRadius; // We go back to the coordinate system centered at the bottom of the can // direction: newVec = ApplyRotationMatrix({Track->D1,Track->D2,Track->D3}, RotMatrix); Track->D1 = newVec[0]; Track->D2 = newVec[1]; Track->D3 = newVec[2]; // normalise the direction cosines to D1^2+D2^2+D3^3=1 double DNorm = sqrt(Track->D1*Track->D1 + Track->D2*Track->D2 + Track->D3*Track->D3); Track->D1 /= DNorm; Track->D2 /= DNorm; Track->D3 /= DNorm; } //___________________________________________________________________________ array,3> GGenerateEvent::ComputeRotation(GSeaTrack *Track, array PointAtCan) { // This function computes the rotation matrix for a given Track around the center of the Earth // needed to intersect a specified PointAtCan (x,y,z) // for convenience, we switch to the coordinate system centered at center of the Earth: PointAtCan[2] += fGenPar->SeaBottomRadius; Track->V3 += fGenPar->SeaBottomRadius; double Track2 = Track->V1*Track->V1 + Track->V2*Track->V2 + Track->V3*Track->V3; double Point2 = PointAtCan[0]*PointAtCan[0] + PointAtCan[1]*PointAtCan[1] + PointAtCan[2]*PointAtCan[2]; // First we find the point, where the initial primary trajectory intersects // the PointAtCan level. We use eq. for intersection of a line and a sphere // ( R=sqrt((R_Earth+h_can)^2+x_can^2+y_can^2) ) // The equations can be found e.g. at http://paulbourke.net/geometry/circlesphere/ // as the two points on the line we pick x1 = Track->V1 and x2 = Track->V1+1*Track->D1 (and the same for V2 and V3) // This simplifies a to: a = Track->D1*Track->D1 + Track->D2*Track->D2 + Track->D3*Track->D3 = 1 double b = 2.*( Track->D1*Track->V1 + Track->D2*Track->V2 + Track->D3*Track->V3 ); double c = Track2 - Point2; double dist = 0.5*(-b-sqrt(b*b-4.*c)); // equations simplify and t=dist. We choose the solution // closer to the track position, since the other intersection point // would be on the other side of the Earth // intersection at same level as PointAtCan: double x = Track->V1 + Track->D1*dist; double y = Track->V2 + Track->D2*dist; double z = Track->V3 + Track->D3*dist; // Now we want to find the rotation matrix rotating this point to the specified PointAtCan. // For this we need the vector ort, orthogonal to vectors pointing at the two rotated points (u and v). // It will be the rotation axis. ort = u x v (cross-product). // based on: https://math.stackexchange.com/questions/2754627/rotation-matrices-between-two-points-on-a-sphere array ort = { y*PointAtCan[2]-z*PointAtCan[1], z*PointAtCan[0]-x*PointAtCan[2], x*PointAtCan[1]-y*PointAtCan[0], }; double ort_norma = sqrt(ort[0]*ort[0] + ort[1]*ort[1] + ort[2]*ort[2]); // for the rotation, we only need the sine and cosine of the rotation angle and the ort vector: double u_len = sqrt(x*x + y*y + z*z); double v_len = sqrt(Point2); double uv_len_product = u_len * v_len; double s = ort_norma/uv_len_product; c = (x*PointAtCan[0]+y*PointAtCan[1]+z*PointAtCan[2])/uv_len_product; double one_m_c = 1.-c; // To obtain the rotation matrix, we need the unit vector of ort: ort[0] /= ort_norma; ort[1] /= ort_norma; ort[2] /= ort_norma; // The rotation matrix is the obtained according to // https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle double ort_xy = ort[0]*ort[1]; double ort_xz = ort[0]*ort[2]; double ort_yz = ort[1]*ort[2]; array,3> RotMatrix = {{ // matrix of rotation {c+ort[0]*ort[0]*one_m_c, ort_xy*one_m_c-ort[2]*s, ort_xz*one_m_c+ort[1]*s}, {ort_xy*one_m_c+ort[2]*s, c+ort[1]*ort[1]*one_m_c, ort_yz*one_m_c-ort[0]*s}, {ort_xz*one_m_c-ort[1]*s, ort_yz*one_m_c+ort[0]*s, c+ort[2]*ort[2]*one_m_c}, }}; // we switch back to the coordinate system centered at the bottom of the can: PointAtCan[2] -= fGenPar->SeaBottomRadius; Track->V3 -= fGenPar->SeaBottomRadius; return RotMatrix; } //___________________________________________________________________________ array,3> GGenerateEvent::ShiftShower() { // This function computes the rotation necessary to have a primary // (V1,V2,V3) position and (D1,D2,D3) direction such that its trajectory would cross a // point randomly sampled on the surface of the can. RandomGen * rnd = RandomGen::Instance(); GSeaTrack primary = fFileDriver->OtherTracks[0]; // start with the original track read from CORSIKA array,3> RotMatrix; double RadDep = fGenPar->SeaBottomRadius+fGenPar->SiteDepth; double x,y,z; if (fGenPar->RTOpt.compare("proj")==0) { // proj option // we only compute for primary to rotate the shower as a whole: // the target point is the (0,0,h_can/2) RotMatrix = ComputeRotation(&primary, {0.,0.,0.5*fGenPar->Can2}); RotateAroundEarth(&primary, RotMatrix); // we perform the rotation on the primary track // the can cyllinder is increased by the xy and z projections of fDistaMax: double sinD3 = sqrt((1.-primary.D3)*(1.+primary.D3)); double height = (fGenPar->Can2-fGenPar->Can1) + fSeaEvent->DistaMax * sinD3; double radius = fGenPar->Can3 + fSeaEvent->DistaMax * fabs(primary.D3); double x_can,y_can,z_can; // position sampled at the can surface double y_can2; // the projected area of the cyllinder (according to e.g. eq. 10 in DOI:10.1364/AO.42.006710): double AreaTop = kPi * radius * radius * fabs(primary.D3) ; double AreaSide = 2. * height * radius * sinD3; fSeaEvent->Agen = AreaTop + AreaSide; // we either draw on top cap or on the side. The probability is weighted by the projected areas if ( rnd->RndFlux().Rndm() < AreaTop / fSeaEvent->Agen ) { // top cap of the cyllinder double r = radius * sqrt( rnd->RndFlux().Rndm() ); double phi = 2. * kPi * rnd->RndFlux().Rndm(); x = r*sin(phi); y = r*cos(phi); z = primary.D3 > 0 ? 0.0 : height; x_can = x; y_can = y; z_can = z; } else { // side of the cyllinder // this only randomly picks a point of the half of the side z_can = height * rnd->RndFlux().Rndm(); y_can = radius * 2. * (rnd->RndFlux().Rndm()-0.5); y_can2 = y_can*y_can; x_can = -sqrt( radius*radius - y_can2 ); // half because we only take one possible value for x // this rotates this half according to primary xy orientation (we can only hit the half that is in front of the primary) double phidir = atan2( primary.D2, primary.D1 ); double sphidir = sin(phidir); double cphidir = cos(phidir); x = cphidir*x_can - sphidir*y_can; y = sphidir*x_can + cphidir*y_can; x_can = x; y_can = y; } // Now that a random position on the can surface is picked, the actual // shifted primary 1st interaction position and direction can be computed: primary = fFileDriver->OtherTracks[0]; // we reset the primary to the nominal position and direction // (because later we want to rotate the muons in the same way!) RotMatrix = ComputeRotation(&primary, {x_can,y_can,z_can}); return RotMatrix; } else if (fGenPar->RTOpt.compare("can")==0) { // can option double zDif; double Can21 = fGenPar->Can2-fGenPar->Can1; double RT = 0.5*sqrt(Can21*Can21+4.*fGenPar->Can3*fGenPar->Can3); // diagonal of the can double Xc=0.,Yc=0.; double Zc=-fGenPar->SeaBottomRadius; // generate the primary vertex (fX0,fY0,fZ0: detector -origin coordinates) double xDif = fFileDriver->GetfX0()-Xc; double yDif = fFileDriver->GetfY0()-Yc; zDif = fFileDriver->GetfZ0()-Zc; double b = fFileDriver->OtherTracks[0].D1*xDif + fFileDriver->OtherTracks[0].D2*yDif + fFileDriver->OtherTracks[0].D3*zDif; double c = xDif*xDif+yDif*yDif+zDif*zDif-RadDep*RadDep; double RL = fabs(-b-sqrt(b*b-c)); RT+=fSeaEvent->DistaMax; fSeaEvent->Agen = kPi*RT*RT; x += -RL * fFileDriver->OtherTracks[0].D1; // we have to multiply RL times the arrival direction cosines y += -RL * fFileDriver->OtherTracks[0].D2; z += -RL * fFileDriver->OtherTracks[0].D3; TVector3 vec(x,y,z); // vector towards selected point TVector3 dvec1 = vec.Orthogonal(); // orthogonal vector TVector3 dvec2 = dvec1; // second orthogonal vector dvec2.Rotate(-0.5*kPi,vec); // rotate second vector by 90deg, now forming a new orthogonal cartesian coordinate system double psi = 2.*kPi* rnd->RndFlux().Rndm(); // rndm angle [0,2pi] double random = rnd->RndFlux().Rndm(); // rndm number [0,1] dvec1.SetMag(sqrt(random)*RT*cos(psi)); dvec2.SetMag(sqrt(random)*RT*sin(psi)); x += dvec1.X() + dvec2.X(); y += dvec1.Y() + dvec2.Y(); z += dvec1.Z() + dvec2.Z(); x += fFileDriver->GetfX0() ; y += fFileDriver->GetfY0(); z += fFileDriver->GetfZ0(); RotMatrix[0]={x,y,z}; return RotMatrix; } else { LOG("GGenerateEvent", pFATAL) << "Unknown generation surface option "<RTOpt<<", exiting..."; exit(1); } } //________________________________________________________________________________________________________ bool GGenerateEvent::PropagateTracks() { TrksInCan = false; NeedToShift = false; NeedToPropagate = true; AboveCan = true; fSeaEvent->NRetries=0; bool FirstStage = true; // first stage of propagation (from sea to the can top) RandomGen * rnd = RandomGen::Instance(); double RadDep = fGenPar->SeaBottomRadius+fGenPar->SiteDepth; double PrimX=0,PrimY=0,PrimZ=0; // needed to maintain the -rt 'can' option (sub-optimal, -rt 'proj' should be used by default) double V1diff,V2diff,V3diff; vector TempTracks; array,3> RotMatrix,OriginRotMatrix; array xyzShift={0.,0.,0.}; if ( fGenPar->GenCode.compare("CORSIKA")==0 ) { // CORSIKA tracks TempTracks = fFileDriver->SecondaryTracks; OriginRotMatrix = ShiftShower(); // we compute the rotation of the shower if (fGenPar->RTOpt.compare("proj")==0) { // we perform the rotation on the secondary tracks only (the rest will be rotated if the event is succesful): for ( auto& track : TempTracks ) RotateAroundEarth(&track, OriginRotMatrix); } else if (fGenPar->RTOpt.compare("can")==0) { xyzShift = ShiftShower()[0]; for ( auto& track : TempTracks ) { // do shifting double Xc=0.,Yc=0.; double Zc=-fGenPar->SeaBottomRadius; double xOff = track.V1-Xc; double yOff = track.V2-Yc; double zOff = track.V3-Zc; double bmu = track.D1*xOff+ track.D2*yOff +track.D3*zOff; double cmu = xOff*xOff+yOff*yOff+zOff*zOff-RadDep*RadDep; double RLmu = -bmu-sqrt(bmu*bmu-cmu); track.V1 += xyzShift[0] + fFileDriver->GetfX0() + RLmu * track.D1; track.V2 += xyzShift[1] + fFileDriver->GetfY0() + RLmu * track.D2; track.V3 += xyzShift[2] + fFileDriver->GetfZ0() + RLmu * track.D3; } } } else TempTracks = fSeaEvent->Tracks; // GENIE tracks int nGenTracks = TempTracks.size(); do { if (TrksInCan && !NeedToPropagate) { // muons inside can and all are propagated, exit and save event LOG("GGenerateEvent", pDEBUG) << "Succesfully transported tracks into the can"; if ( fGenPar->GenCode.compare("CORSIKA")==0 ) { fFileDriver->SecondaryTracks = TempTracks; if (fGenPar->RTOpt.compare("proj")==0) { // do shifting for OtherTracks for ( auto& track : fFileDriver->OtherTracks ) RotateAroundEarth(&track, OriginRotMatrix); } else if (fGenPar->RTOpt.compare("can")==0) { xyzShift[0] = PrimX - fFileDriver->OtherTracks[0].V1; xyzShift[1] = PrimY - fFileDriver->OtherTracks[0].V2; xyzShift[2] = PrimZ - fFileDriver->OtherTracks[0].V3; for ( auto& track : fFileDriver->OtherTracks ) { // do shifting for OtherTracks double Xc=0.,Yc=0.; double Zc=-fGenPar->SeaBottomRadius; double xOff = track.V1-Xc; double yOff = track.V2-Yc; double zOff = track.V3-Zc; double bmu = track.D1*xOff+ track.D2*yOff +track.D3*zOff; double cmu = xOff*xOff+yOff*yOff+zOff*zOff-RadDep*RadDep; double RLmu = -bmu-sqrt(bmu*bmu-cmu); track.V1 += xyzShift[0] + fFileDriver->GetfX0() + RLmu * track.D1; track.V2 += xyzShift[1] + fFileDriver->GetfY0() + RLmu * track.D2; track.V3 += xyzShift[2] + fFileDriver->GetfZ0() + RLmu * track.D3; } } } else TempTracks = fSeaEvent->Tracks = TempTracks; // GENIE Tracks vector().swap(TempTracks); // clear memory if (fSeaEvent->NRetries>0) LOG("GGenerateEvent", pDEBUG) << "Event was retried " << fSeaEvent->NRetries << "/" << fGenPar->ChancesPerShower << " times."; return TrksInCan; } else if (NeedToPropagate) { // propagate muons LOG("GGenerateEvent", pDEBUG) << "(re-)propagating ["<<(fGenPar->ChancesPerShower-fSeaEvent->NRetries)<<"/"<ChancesPerShower<<"] Chances left"; NeedToPropagate = false; NeedToShift = false; // first propagate the taus (we need to propagate all of them) #ifdef _HETAU_ENABLED__ for ( int i=0; i decayed; if( TempTracks[i].Status==1 ) { //propagate all final state taus inside can decayed = DecayTau(TempTracks[i],-log( rnd->RndGen().Rndm() )*lifetime*1e9); double V1diff = decayed[0].V1-TempTracks[i].V1; double V2diff = decayed[0].V2-TempTracks[i].V2; double V3diff = decayed[0].V3-TempTracks[i].V3; TempTracks[i].Length = sqrt(V1diff*V1diff + V2diff*V2diff + V3diff*V3diff); } else if( TempTracks[i].Status==-1) { //outside the can auto propagated_tau = TempTracks[i]; double tauti = PropagationTau(&propagated_tau); if ( tauti==0 ) { // decayed before got to the can TempTracks[i].SetStatus(2001); decayed = DecayTau(propagated_tau,tauti); double V1diff = decayed[0].V1-TempTracks[i].V1; double V2diff = decayed[0].V2-TempTracks[i].V2; double V3diff = decayed[0].V3-TempTracks[i].V3; TempTracks[i].Length = sqrt(V1diff*V1diff + V2diff*V2diff + V3diff*V3diff); } else if ( tauti>0 ) { // didnt decay before got to the can TempTracks[i].SetStatus(1001); double V1diff = propagated_tau.V1-TempTracks[i].V1; double V2diff = propagated_tau.V2-TempTracks[i].V2; double V3diff = propagated_tau.V3-TempTracks[i].V3; TempTracks[i].Length = sqrt(V1diff*V1diff + V2diff*V2diff + V3diff*V3diff); propagated_tau.MotherId = TempTracks[i].Id; if ( fGenPar->GenCode.compare("CORSIKA")==0 ) { // CORSIKA tracks propagated_tau.Id = TempTracks.size() + fFileDriver->OtherTracks.size(); // post-increment (without OtherTracks there will be duplicate Ids!) } else propagated_tau.Id = TempTracks.size(); // post-increment decayed = DecayTau(propagated_tau,tauti); V1diff = decayed[0].V1-propagated_tau.V1; V2diff = decayed[0].V2-propagated_tau.V2; V3diff = decayed[0].V3-propagated_tau.V3; propagated_tau.Length = sqrt(V1diff*V1diff + V2diff*V2diff + V3diff*V3diff); TempTracks.push_back( propagated_tau ); } } for ( unsigned int k=0; kGenCode.compare("CORSIKA")==0 ) { // CORSIKA tracks decayed[k].Id = TempTracks.size() + fFileDriver->OtherTracks.size(); } else decayed[k].Id = TempTracks.size(); // post-increment TempTracks.push_back( decayed[k] ); } } nGenTracks = TempTracks.size(); // we update the size because there could be new muons from taus #endif // now go to muons for ( int i=0; iGenCode.compare("CORSIKA")==0 ) { // CORSIKA tracks TempTracks[i] = fFileDriver->SecondaryTracks[i]; // redo shifting: if (fGenPar->RTOpt.compare("proj")==0) RotateAroundEarth(&TempTracks[i], OriginRotMatrix); // redo the rotation else if (fGenPar->RTOpt.compare("can")==0) { TempTracks[i].V1 += xyzShift[0]; TempTracks[i].V2 += xyzShift[1]; TempTracks[i].V3 += xyzShift[2]; } } else TempTracks[i] = unpropagated_mu; // GENIE tracks V1diff = propagated_mu.V1-TempTracks[i].V1; V2diff = propagated_mu.V2-TempTracks[i].V2; V3diff = propagated_mu.V3-TempTracks[i].V3; propagated_mu.Length = sqrt(V1diff*V1diff + V2diff*V2diff + V3diff*V3diff); propagated_mu.MotherId = TempTracks[i].Id; if ( fGenPar->GenCode.compare("CORSIKA")==0 ) { // CORSIKA tracks propagated_mu.Id = TempTracks.size() + fFileDriver->OtherTracks.size(); // post-increment (without OtherTracks there will be duplicate Ids!) } else propagated_mu.Id = TempTracks.size(); // post-increment TempTracks.push_back( propagated_mu ); TempTracks[i].Length = 0.; TempTracks[i].SetStatus(1001); // propagated track that got to the can TrksInCan = true; } else if (FirstStage && fGenPar->GenCode.compare("CORSIKA")==0) TempTracks[i] = propagated_mu; // we save the progress } else if ( TempTracks[i].Status==1 ) { TrksInCan = true; if ( abs(TempTracks[i].PdgCode)==311 ) { // K0 (K0Bar) is detected as K0S or KOL TempTracks[i].PdgCode = (rnd->RndFlux().Rndm()>0.5) ? 130 : 310; } } } FirstStage = false; // there should never be muons above can top after propagation } else if (NeedToShift) { // do shifting (at can top) LOG("GGenerateEvent", pDEBUG) << "Trying to shift the event into the can "; NeedToShift = false; TrksInCan = false; NeedToPropagate = false; AboveCan = false; if ( fGenPar->GenCode.compare("CORSIKA")==0 ) { // CORSIKA tracks while (!TrksInCan && !NeedToPropagate) { // we try until we hit the right point fSeaEvent->NRetries++; // we compute the shower shift to a different position (and hence, orientation) // we undo the previous rotation (because the shower is not in the starting position): RotMatrix = OriginRotMatrix; OriginRotMatrix = ShiftShower(); RotMatrix = MultiplyMatrices(OriginRotMatrix,InvertMatrix(RotMatrix)); // here we shift secondaries first before we update the primary position: if (fGenPar->RTOpt.compare("proj")==0) { nGenTracks = TempTracks.size(); for ( int i=0; iSecondaryTracks[i]; // reset the track back to the sea surface RotateAroundEarth(&TempTracks[i], OriginRotMatrix); // redo the rotation (we use OriginRotMatrix because the tracks are reset) V1diff = reshifted_mu.V1-TempTracks[i].V1; V2diff = reshifted_mu.V2-TempTracks[i].V2; V3diff = reshifted_mu.V3-TempTracks[i].V3; reshifted_mu.Length = sqrt(V1diff*V1diff + V2diff*V2diff + V3diff*V3diff); reshifted_mu.MotherId = TempTracks[i].Id; reshifted_mu.Id = TempTracks.size() + fFileDriver->OtherTracks.size(); // post-increment (without OtherTracks there will be duplicate Ids!) TempTracks[i].Length = 0.; TempTracks[i].SetStatus(1001); // track that got to the can TempTracks.push_back( reshifted_mu ); break; } // the track is inside the can case 0: { NeedToPropagate=true; break; } } } } else if (fGenPar->RTOpt.compare("can")==0) { xyzShift = ShiftShower()[0]; for ( int i=0; iSeaBottomRadius; double xOff = TempTracks[i].V1-Xc; double yOff = TempTracks[i].V2-Yc; double zOff = TempTracks[i].V3-Zc; double bmu = TempTracks[i].D1*xOff+ TempTracks[i].D2*yOff +TempTracks[i].D3*zOff; double cmu = xOff*xOff+yOff*yOff+zOff*zOff-RadDep*RadDep; double RLmu = -bmu-sqrt(bmu*bmu-cmu); TempTracks[i].V1 += xyzShift[0] + fFileDriver->GetfX0() + RLmu * TempTracks[i].D1; TempTracks[i].V2 += xyzShift[1] + fFileDriver->GetfY0() + RLmu * TempTracks[i].D2; TempTracks[i].V3 += xyzShift[2] + fFileDriver->GetfZ0() + RLmu * TempTracks[i].D3; if (!TrksInCan && !NeedToPropagate) { // do not check if we already know there are muons in can and muons still to be propagated switch (DistToCan(&TempTracks[i])) { // check if the muon track intercepts the can case 1: { TrksInCan = true; TempTracks[i].SetStatus(1); break; } // the track is inside the can case 0: { NeedToPropagate=true; break; } } } } PrimX = fFileDriver->OtherTracks[0].V1 + xyzShift[0]; PrimY = fFileDriver->OtherTracks[0].V2 + xyzShift[1]; // we only shift horizontally at can top } if (TrksInCan || NeedToPropagate) { // the shifting was succesful LOG("GGenerateEvent", pDEBUG) << "Succesfully shifted!"; break; } else if (fGenPar->ChancesPerShower==fSeaEvent->NRetries) { // no more retries left LOG("GGenerateEvent", pDEBUG) << "No success shifting ..."; break; } else if (fSeaEvent->NRetries%5==0) { // 5 shifting retries and if we fail, move to re-propagating LOG("GGenerateEvent", pDEBUG) << "No success shifting ..."; break; } } } else { // GENIE tracks LOG("GGenerateEvent", pDEBUG) << "Shifting not (yet) implemented for GENIE tracks"; // TODO: implement shifting for GENIE tracks } } else { // propagation and/or shifting failed, need to start over LOG("GGenerateEvent", pDEBUG) << "Propagation and/or shifting failed, starting over ..."; fSeaEvent->NRetries++; TrksInCan = false; NeedToShift = false; NeedToPropagate = true; AboveCan = true; FirstStage = true; // we go back to the first stage // reset the tracks: if ( fGenPar->GenCode.compare("CORSIKA")==0 ) { // CORSIKA tracks TempTracks = fFileDriver->SecondaryTracks; // we rotate the shower to a different position (and hence, orientation): OriginRotMatrix = ShiftShower(); if (fGenPar->RTOpt.compare("proj")==0) { for ( auto& track : TempTracks ) RotateAroundEarth(&track, OriginRotMatrix); // do rotation } else if (fGenPar->RTOpt.compare("can")==0) { xyzShift = ShiftShower()[0]; PrimX = fFileDriver->OtherTracks[0].V1 + xyzShift[0]; PrimY = fFileDriver->OtherTracks[0].V2 + xyzShift[1]; PrimZ = fFileDriver->OtherTracks[0].V3 + xyzShift[2]; for ( auto& track : TempTracks ) { // do shifting double Xc=0.,Yc=0.; double Zc=-fGenPar->SeaBottomRadius; double xOff = track.V1-Xc; double yOff = track.V2-Yc; double zOff = track.V3-Zc; double bmu = track.D1*xOff+ track.D2*yOff +track.D3*zOff; double cmu = xOff*xOff+yOff*yOff+zOff*zOff-RadDep*RadDep; double RLmu = -bmu-sqrt(bmu*bmu-cmu); track.V1 += xyzShift[0] + fFileDriver->GetfX0() + RLmu * track.D1; track.V2 += xyzShift[1] + fFileDriver->GetfY0() + RLmu * track.D2; track.V3 += xyzShift[2] + fFileDriver->GetfZ0() + RLmu * track.D3; } } } else TempTracks = fSeaEvent->Tracks; // GENIE tracks } } while ( fGenPar->ChancesPerShower >= fSeaEvent->NRetries ); if ( fGenPar->GenCode.compare("CORSIKA")==0 ) { // CORSIKA tracks if ( !(fFileDriver->FullRun) && (fSeaEvent->NRetries==fGenPar->ChancesPerShower) ) { // we tried, but failed LOG("GGenerateEvent", pDEBUG) << "Event was retried " << fGenPar->ChancesPerShower << "/" << fGenPar->ChancesPerShower << " times, but still did not reach the can."; } } return TrksInCan; } //________________________________________________________________________________________________________ bool GGenerateEvent::WriteTracks(bool TrksInCan) { //Add particles to the gWriteTracks depending on WriteParticlesMode. if ( fGenPar->GenCode.compare("CORSIKA")==0 ) { if ( fGenPar -> WriteParticlesMode == 2 ) { LOG("GGenerateEvent", pDEBUG) << "-write 2 option selected, saving everything"; } else if ( fGenPar -> WriteParticlesMode == 1 ) { if (!TrksInCan) { vector().swap(fFileDriver->OtherTracks); vector().swap(fFileDriver->SecondaryTracks); } } else if ( fGenPar -> WriteParticlesMode == 0 ) { if (TrksInCan) { for ( auto& track : fFileDriver->SecondaryTracks ) track.MotherId=0; // since there is no other parent, primary is set as mother for all secondaries fFileDriver->OtherTracks.erase(fFileDriver->OtherTracks.begin()+1,fFileDriver->OtherTracks.end()); // we remove everything except the primary from OtherTracks } else { vector().swap(fFileDriver->OtherTracks); vector().swap(fFileDriver->SecondaryTracks); } } return fFileDriver->OtherTracks.size() > 0; } else { vector().swap(gWrittenTracks); // clear memory if ( fGenPar -> WriteParticlesMode == 2 ) { gWrittenTracks = fSeaEvent->Tracks; } else if ( fGenPar -> WriteParticlesMode == 1 ) { if (TrksInCan) gWrittenTracks = fSeaEvent->Tracks; } else if ( fGenPar -> WriteParticlesMode == 0 && TrksInCan ) { for ( auto& track : fSeaEvent->Tracks ) { bool keep = track.Id==0 || track.Status == 5 || track.Status>1000 || track.Status==1 ; if (!keep) // find all track's children and set their MotherId to track's mother { for ( auto& t : fSeaEvent->Tracks ) { if ( t.MotherId == track.Id && track.MotherId!=-1) t.MotherId = track.MotherId; } } } for ( auto& track : fSeaEvent->Tracks ) { bool keep = track.Id==0 || track.Status == 5 || track.Status>1000 || track.Status==1 ; // repeated code :( if (keep) gWrittenTracks.push_back(track); } } return gWrittenTracks.size() > 0; } } //___________________________________________________________________________ double GGenerateEvent::GetXSecWater(int neutrino_pdg, double Energy){ double XSec=0.; map::iterator iter = fGenPar->XSecWater.begin(); for ( ; iter != fGenPar->XSecWater.end(); ++iter) { if(neutrino_pdg==iter->first){ TGraph GrSecTot = iter->second; XSec=GrSecTot.Eval(Energy); break; } } return XSec; } //___________________________________________________________________________ double GGenerateEvent::GetXSecSW(int neutrino_pdg, double Energy){ double XSec=0.; map::iterator iter = fGenPar->XSecTotSW.begin(); for ( ; iter != fGenPar->XSecTotSW.end(); ++iter) { if(neutrino_pdg==iter->first){ TGraph GrSecTot = iter->second; XSec=GrSecTot.Eval(Energy); break; } } return XSec; }