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            1 :         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            1 :         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            2 : }
      89              : 
      90              : } //namespace crpropa