Line data    Source code

       1             : #include "crpropa/Units.h"
       2             : #include "crpropa/Common.h"
       3             : #include "crpropa/ParticleID.h"
       4             : #include "crpropa/module/CandidateSplitting.h"
       5             : 
       6             : #include "gtest/gtest.h"
       7             : #include <stdexcept>
       8             : #include <cmath>
       9             : 
      10             : namespace crpropa {
      11             : 
      12           1 : TEST(testCandidateSplitting, SimpleTest) {
      13           1 :         int nSplit = 2;
      14             :         int nBins = 4;
      15             :         double minWeight = pow(1. / nSplit, 2);
      16             :         double Emin = 1; // dimensionless for testing
      17             :         double Emax = 10;       
      18             : 
      19           1 :         CandidateSplitting split_lin(nSplit, Emin, Emax, nBins, minWeight);
      20             :         double dE = (Emax - Emin) / nBins;
      21           1 :         EXPECT_DOUBLE_EQ(split_lin.getEnergyBins()[0], Emin);
      22           1 :         EXPECT_DOUBLE_EQ(split_lin.getEnergyBins()[1], Emin + dE);
      23             : 
      24           1 :         EXPECT_EQ(split_lin.getNsplit(), nSplit);
      25           1 :         EXPECT_DOUBLE_EQ(split_lin.getMinimalWeight(), minWeight);
      26             : 
      27           2 :         CandidateSplitting split_log(nSplit, Emin, Emax, nBins, minWeight, true);
      28             :         double dE_log = pow(Emax / Emin, 1. / (nBins - 1.0));
      29           1 :         EXPECT_DOUBLE_EQ(split_log.getEnergyBins()[0], Emin);
      30           1 :         EXPECT_DOUBLE_EQ(split_log.getEnergyBins()[1], Emin * dE_log);
      31             : 
      32             :         double spectralIndex = -2.;
      33           2 :         CandidateSplitting split_dsa(spectralIndex, Emin, nBins);
      34             :         double dE_dsa = pow(1. / 2, 1. / (spectralIndex + 1));
      35           1 :         EXPECT_DOUBLE_EQ(split_dsa.getEnergyBins()[0], Emin * dE_dsa);
      36           1 :         EXPECT_DOUBLE_EQ(split_dsa.getEnergyBins()[nBins - 1], Emin * pow(dE_dsa, nBins));
      37           1 : }
      38             : 
      39             : 
      40           1 : TEST(testCandidateSplitting, CheckSplits) {
      41             :         int nSplit = 2;
      42             :         int nBins = 3;
      43             :         double Emin = 1; // dimensionless for testing
      44             :         double Emax = 10;
      45             :         double minWeight = pow(1. / nSplit, 4);
      46             : 
      47           1 :         CandidateSplitting splitting(nSplit, Emin, Emax, nBins, minWeight);
      48           1 :         Candidate c(nucleusId(1,1),0.5);
      49             :         double weight = 1.0;
      50           1 :         double serial = c.getSerialNumber();
      51             :         
      52           1 :         splitting.process(&c); // no split
      53           1 :         EXPECT_DOUBLE_EQ(c.getWeight(), weight);
      54           1 :         EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial);
      55             : 
      56           1 :         c.current.setEnergy(2); 
      57           1 :         splitting.process(&c); // 1. split
      58             :         weight = weight/nSplit;
      59           1 :         EXPECT_DOUBLE_EQ(c.getWeight(), weight);
      60           1 :         EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 1);
      61           1 :         c.previous.setEnergy(2);
      62             : 
      63           1 :         c.current.setEnergy(6); 
      64           1 :         splitting.process(&c); // 2. split
      65             :         weight = weight/nSplit;
      66           1 :         EXPECT_DOUBLE_EQ(c.getWeight(), weight);
      67           1 :         EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 2);
      68           1 :         c.previous.setEnergy(6);
      69             : 
      70           1 :         c.current.setEnergy(0.5); 
      71           1 :         splitting.process(&c); // no split, cooling
      72           1 :         EXPECT_DOUBLE_EQ(c.getWeight(), weight);
      73           1 :         EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 2);
      74           1 :         c.previous.setEnergy(0.5);
      75             : 
      76           1 :         c.current.setEnergy(6); 
      77           1 :         splitting.process(&c); // 3. & 4. split, crosses two boundaries
      78             :         weight = weight/nSplit/nSplit;
      79           1 :         EXPECT_DOUBLE_EQ(c.getWeight(), weight);
      80           1 :         EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 4);
      81           1 :         c.previous.setEnergy(6);
      82             : 
      83           1 :         c.current.setEnergy(8); 
      84           1 :         splitting.process(&c); // no split, minimal weight reached
      85           1 :         EXPECT_DOUBLE_EQ(c.getWeight(), weight);
      86           1 :         EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 4);
      87           1 :         c.previous.setEnergy(8);
      88           1 : }
      89             : 
      90             : } //namespace crpropa