nucleilist.cc

/**
 * @file nucleilist.cc
 * @author Luca Maccione
 * @email luca.maccione@desy.de
 * @brief Class TNucleiList is implemented. @sa nucleilist.h
 */

#include "nucleilist.h"
#include "constants.h"
#include "utilities.h"
#include "input.h"

#include <algorithm>
#include <stdlib.h>
#include <fstream>
#include <string>

using namespace std;

/**
 * @fn bool comparator(int i, int j)
 * @brief Function defining a weak ordering of nuclei, so that they are propagated from the heaviest to the lightest, with the appropriate inversions to take into account beta decays.
 * @param i uid of the first nucleus.
 * @param j uid of the second nucleus.
 */
bool comparator(int i, int j) {
  int Ai = -1000;
  int Zi = -1000;
  Utility::id_nuc(i, Ai, Zi);
   
  int Aj = -1000;
  int Zj = -1000;
  Utility::id_nuc(j, Aj, Zj);
   
  if ( (Ai == Aj) && ( (Zi == 4 && Zj == 5) || (Zi == 6 && Zj == 7) ||(Zi == 14 && Zj == 16)) ) return true;
  if (Zi == Zj && (Ai>Aj)) return true;
  if (Zi>Zj) return true;
  return false;
}

DECMODE strtoDecMode(const string& str) {
  if (str == "EC") return EC;
  if (str == "B-") return BM;
  if (str == "BB") return BB;
  if (str == "IT") return IT;
  if (str == "ECB-") return ECBM;
  if (str == "ECB+") return ECBP;
  return STABLE;
}

void TNucleiList::DeleteNucleusFromList(const int& uid_){
  for (vector<int>::iterator it = uid.begin(); it != uid.end(); ++it)
    if (*it==uid_) uid.erase(it);
  
  return;
}

TNucleiList::TNucleiList(Input* in) {
  ifstream infile(BNLdata.c_str(), ios::in); // Read databasis
  
  char s[3000];
  infile.getline(s,3000);
  
  int A = 0;
  int Z = 0;
  double tau = 0.;
  double tau_naked = 0.;
  string Mode;
  string name;
  
  int id;
  while (infile >> A >> Z >> tau >> tau_naked >> Mode >> name) {
    
    if (Z > in->Zmax || Z < in->Zmin) continue; /**< Select only nuclei in the wished charge range. \sa constants.h */
    id = 1000*Z+A; // Compute uid
    if (tau > minlifetime || Mode == "STABLE") {
      uid.push_back(id);
      lifetime[id]       = tau/Myr;
      lifetime_naked[id] = tau_naked/Myr;
      //cout << "...adding to the list this nucleus: id = " << id << "; lifetime in Myr = " << tau/Myr << endl;
      betadec[id] = strtoDecMode(Mode);
    }

    else id_short.push_back(id); // Short lived nucleus; THIS LIST IS NOT USED. Short-lived nuclei are not propagated but are treated properly in the spallation network
  }
   
  infile.close();
   
  sort(uid.begin(), uid.end(), comparator); // Sort nuclei according to ordering defined by function comparator.
  sort(id_short.begin(), id_short.end(), comparator);
   
  if (in->prop_ap || in->prop_DMap) {
    uid.push_back(1000*(-1)+1); // Secondary antiprotons
    lifetime[1000*(-1)+1] = -1;
    betadec[1000*(-1)+1] = STABLE;
  }
  
  if (in->prop_lep || in->prop_extracomp || in->prop_DMel) {
    uid.push_back(1000*(-1)); // Primary electrons
    lifetime[1000*(-1)] = -1;
    betadec[1000*(-1)] = STABLE;
      
    uid.push_back(1000*(1)); // Secondary positrons
    lifetime[1000*(1)] = -1;
    betadec[1000*(1)] = STABLE;
  }
  
  if (in->prop_deuteron) {
    uid.push_back(1000*(-1)+2); // secondary antideuterons
    lifetime[1000*(-1)+2] = -1;
    betadec[1000*(-1)+2] = STABLE;
  }
  
  return ;
}