/* * Particle.h * * Author: * Oleg Kalashev * * Copyright (c) 2020 Institute for Nuclear Research, RAS * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal * in the Software without restriction, including without limitation the rights * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * copies of the Software, and to permit persons to whom the Software is * furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included in * all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN * THE SOFTWARE. */ #ifndef PARTICLE_H #define PARTICLE_H #include "Cosmology.h" #include <cstring> #include "Units.h" #include <math.h> namespace mcray { enum ParticleType {// order of masses in the array Particle::massesMeV must follow the order of particles below Electron = 0, Positron, Photon, NeutrinoE, NeutrinoM, NeutrinoT, NeutrinoAE, NeutrinoAM, NeutrinoAT, Neutron, Proton, EndLightParticle, StartNuclei = EndLightParticle, A02 = StartNuclei, A03, A04, A05, A06, A07, A08, A09, A10, A11, A12, A13, A14, A15, A16, A17, A18, A19, A20, A21, A22, A23, A24, A25, A26, A27, A28, A29, A30, A31, A32, A33, A34, A35, A36, A37, A38, A39, A40, A41, A42, A43, A44, A45, A46, A47, A48, A49, A50, A51, A52, A53, A54, A55, A56, EndNuclei, ParticleTypeEOF = EndNuclei //must be the last }; typedef cosmo_time coord_type; #define ParticleCustomDataSize 16 class Particle { private: inline void Reset() { memset(this,0,sizeof(Particle));Weight=1.; LastB = -1.; Pdir[2]=1.; } void GenerateId(); public: inline double Mass() const { return Mass(Type); } inline int ElectricCharge() const { return ElectricCharge(Type); } //Construct primary particle inline Particle(ParticleType aType, cosmo_time aZsource){ Reset(); Type=aType; fProductionTime.setZ(aZsource) ; Time.setZ(aZsource); SourceParticle = this; GenerateId(); } inline Particle& operator=(const Particle& aParticle) { memcpy(this, &aParticle, sizeof(Particle)); return *this; } inline coord_type beta() const { coord_type m=Mass(); coord_type gamma_1 = m/Energy; coord_type beta2 = 1. - gamma_1*gamma_1; return sqrt(beta2); } inline void SetParent(const Particle& aParticle){ GenerateId(); if(id != aParticle.id){ fPrevId = aParticle.id; fProductionTime = Time; } } static double Mass(ParticleType aType); static double ElectricCharge(ParticleType aType); void PropagateFreely(cosmo_time dt); ParticleType Type; double Weight; CosmoTime Time; double Energy; long Ninteractions;//number of interactions in the interaction chain coord_type X[3];//comoving coordinates, source location: (0,0,0) coord_type Pdir[3];//momentum direction, initial value: (0,0,1) unsigned long long id; unsigned long long fPrevId; // TODO: interaction specific attributes should be stored separately in custom structures which can // be accessed via pointers stored in interactionData field void* interactionData[ParticleCustomDataSize];//used to store interaction specific attributes as pointers static int ReserveInteractionDataSlot();//returns index to access and store interaction data or -1 if no space left const Particle* SourceParticle; double Deflection2;//uncorrelated deflection angle squared double CorrelatedDeflection;//current value of correlated deflection mutable double CorrelatedBpath;// travel distance within B-field coherence length CosmoTime fProductionTime; // time when the particle was produced either in source or on the way as secondary from interaction; CosmoTime fCascadeProductionTime; // time when primary EM cascade particle was produced by hadron double dt; // particle time delay; mutable double LastB; //last transverce magnetic field within B-field coherence length inline double JetOpenningAngle() const { return DeflectionAngle()-ObservationAngle(); } inline CosmoTime LastDeflectionTime() const { return ElectricCharge() ? Time : fProductionTime; } inline double ObservationAngle() const { double beta = DeflectionAngle(); if(beta==0.) return 0.; double sinBeta=sin(beta); double pathToLastDeflection = LastDeflectionTime().t()-SourceParticle->Time.t();//distance from the source for the parent electron/positron (valid for small z source and deflection) double distanceToSource = Time.t()-SourceParticle->Time.t();// distance between the source and the observer (valid for small z source) return asin(pathToLastDeflection/distanceToSource*sinBeta);// observation angle } inline double DeflectionAngle() const { return sqrt(Deflection2 + CorrelatedDeflection*CorrelatedDeflection); } inline double TimeDelay() const { if(ElectricCharge()) return 0.;//the approximation used below works only for neutral particles ( LastDeflectionTime() > Time.t() ) double sinBeta=sin(DeflectionAngle()); double pathToLastDeflection = LastDeflectionTime().t()-SourceParticle->Time.t();//xx - distance from the source for the parent electron/positron (valid for small z source and deflection) double distanceToSource = Time.t()-SourceParticle->Time.t();// distance between the source and the observer (valid for small z source) return 2*pathToLastDeflection*(1.0-pathToLastDeflection/distanceToSource)*sinBeta*sinBeta;// time delay } inline static const char* Name(ParticleType aType) { return ParticleNames[aType]; } std::string ToString(); private: static const double MassesMeV[EndLightParticle]; static const int ElectricCharges[EndLightParticle]; static const char* ParticleNames[ParticleTypeEOF]; static int fNextCustomDataIndex; static unsigned long long fMaxId; }; } #endif /* PARTICLE_H */