Skip to content
Snippets Groups Projects
TExogamPhysics.cxx 48.2 KiB
Newer Older
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(0.964079);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(14.6);
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(1.11207);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(13.64);
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(0.778904);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(12.94);
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(1.08587);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(10.21);
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(0.244698);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(7.58);
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(0.867378);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(4.25);
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(0.443965);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(2.82);
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(0.411116);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(2.23);
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(1.08974);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(1.73);
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(1.29914);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(1.62);
    Source_isotope.push_back("$^{152}$Eu");
    Source_E.push_back(1.21295);
    Source_Sig.push_back(0.0001);
    Source_branching_ratio.push_back(1.42);
  }
  else{
    std::cout << "Please enter a valid source for gamma ray calibration\nCurrently supported sources are 60Co and 152Eu\n";
    exit(1);
  }
}


/////////////////////////////////////////////////////////////////////
// FIXME Probably could be done better, currently a but inelegant
/////////////////////////////////////////////////////////////////////
void TExogamPhysics::CreateCalibrationEFiles(ofstream* calib_file,
  std::string Path = NPOptionManager::getInstance()->GetCalibrationOutputPath();
  std::string OutputName = NPOptionManager::getInstance()->GetOutputFile();
  if (OutputName.size() > 5) {
    if (OutputName.substr(OutputName.size() - 5, OutputName.size()) == ".root") {
      OutputName = OutputName.substr(0, OutputName.size() - 5);
    }
  }
  TString Filename = "Cal_Exogam_E";
  (*calib_file).open(((string)(Path + OutputName + "/Exogam_E/" + Filename + ".cal")).c_str());
  (*dispersion_file).open(((string)(Path + OutputName + "/Exogam_E/" +Filename + ".dispersion")).c_str());
}

void TExogamPhysics::CreateCalibrationEHGFiles(ofstream* calib_file,
  std::string Path = NPOptionManager::getInstance()->GetCalibrationOutputPath();
  std::string OutputName = NPOptionManager::getInstance()->GetOutputFile();
  if (OutputName.size() > 5) {
    if (OutputName.substr(OutputName.size() - 5, OutputName.size()) == ".root") {
      OutputName = OutputName.substr(0, OutputName.size() - 5);
    }
  }
  TString Filename = "Cal_Exogam_EHG";
  (*calib_file).open(((string)(Path + OutputName + "/Exogam_EHG/" + Filename + ".cal")).c_str());
  (*dispersion_file).open(((string)(Path + OutputName + "/Exogam_EHG/" +Filename + ".dispersion")).c_str());
}

void TExogamPhysics::CreateCalibrationOuterFiles(ofstream* calib_file,
  std::string Path = NPOptionManager::getInstance()->GetCalibrationOutputPath();
  std::string OutputName = NPOptionManager::getInstance()->GetOutputFile();
  if (OutputName.size() > 5) {
    if (OutputName.substr(OutputName.size() - 5, OutputName.size()) == ".root") {
      OutputName = OutputName.substr(0, OutputName.size() - 5);
    }
  }
  TString Filename = "Cal_Exogam_Outer";
  (*calib_file).open(((string)(Path + OutputName + "/Exogam_Outer/" + Filename + ".cal")).c_str());
  (*dispersion_file).open(((string)(Path + OutputName + "/Exogam_Outer/" +Filename + ".dispersion")).c_str());
}

void TExogamPhysics::ReadDoCalibration(NPL::InputParser parser) {
  vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Exogam");

  vector<string> calibs = {"Threshold_E","Threshold_EHG","Threshold_Outers","FirstCr","LastCr","FirstOuter","LastOuter","FitOrder","Source"};

  for (unsigned int i = 0; i < blocks.size(); i++) {
    if (blocks[i]->HasTokenList(calibs)) {
      if (NPOptionManager::getInstance()->GetVerboseLevel())
        cout << endl << "////  Exogam Calibration" << endl;
      unsigned int FirstCr = blocks[i]->GetInt("FirstCr");
      unsigned int LastCr = blocks[i]->GetInt("LastCr");
      unsigned int FirstOuter = blocks[i]->GetInt("FirstOuter");
      unsigned int LastOuter = blocks[i]->GetInt("LastOuter");
      FitPolOrder = blocks[i]->GetInt("FitOrder");
      Source_name = blocks[i]->GetString("Source");
      Threshold_E_Cal = blocks[i]->GetInt("Threshold_E");
      Threshold_EHG_Cal = blocks[i]->GetInt("Threshold_EHG");
      Threshold_Outers_Cal = blocks[i]->GetInt("Threshold_Outers");
      for(unsigned int k = FirstCr; k <= LastCr; k++){
        DoCalibrationE[k] = true;
        DoCalibrationEHG[k] = true; 
        for(unsigned int p = FirstOuter; p <= LastOuter; p++){
          DoCalibrationOuter[k][p] = true; 
      }
    }
    else {
      cout << "ERROR: Missing token for Exogam DoCalibration blocks, check your "
/////////////////////////////////////////////////////////////////////
void TExogamPhysics::WriteHistogramsCalib() {
  std::cout << "Writing Exogam Histograms\n";
  WriteHistogramsE();
  RootHistogramsCalib::getInstance()->GetFile()->Close();
}

/////////////////////////////////////////////////////////////////////
void TExogamPhysics::WriteHistogramsE() {
  auto File = RootHistogramsCalib::getInstance()->GetFile();
  auto TH1Map = RootHistogramsCalib::getInstance()->GetTH1Map();
  auto TGraphMap = RootHistogramsCalib::getInstance()->GetTGraphMap();

  map<int, bool>::iterator it;
  map<int, map<int,bool>>::iterator it2;
  std::string hnameEXOE;

  if (!File->GetDirectory("Exogam"))
    File->mkdir("Exogam");
  File->cd("Exogam");

  for (it = DoCalibrationE.begin(); it != DoCalibrationE.end(); it++) {
    if (it->second) {
      hnameEXOE = Form("EXO_E%d", it->first);
      if (!gDirectory->GetDirectory(Form("EXO_Cr%d", it->first)))
        gDirectory->mkdir(Form("EXO_Cr%d", it->first));
      gDirectory->cd(Form("EXO_Cr%d", it->first));
      (*TH1Map)["Exogam"][hnameEXOE]->Write();
      if((*TGraphMap)["Exogam"][Form("Calib_Graph_EXO_E%d",it->first)]!= nullptr)
        (*TGraphMap)["Exogam"][Form("Calib_Graph_EXO_E%d",it->first)]->Write();
      if((*TGraphMap)["Exogam"][Form("Residue_Graph_EXO_E%d",it->first)]!= nullptr)
        (*TGraphMap)["Exogam"][Form("Residue_Graph_EXO_E%d",it->first)]->Write();
    }
    File->cd("Exogam");
  }
  for (it = DoCalibrationEHG.begin(); it != DoCalibrationEHG.end(); it++) {
    if (it->second) {
      hnameEXOE = Form("EXO_EHG%d", it->first);
      if (!gDirectory->GetDirectory(Form("EXO_Cr%d", it->first)))
        gDirectory->mkdir(Form("EXO_Cr%d", it->first));
      gDirectory->cd(Form("EXO_Cr%d", it->first));
      (*TH1Map)["Exogam"][hnameEXOE]->Write();
      if((*TGraphMap)["Exogam"][Form("Calib_Graph_EXO_EHG%d",it->first)]!= nullptr)
        (*TGraphMap)["Exogam"][Form("Calib_Graph_EXO_EHG%d",it->first)]->Write();
      if((*TGraphMap)["Exogam"][Form("Residue_Graph_EXO_EHG%d",it->first)]!= nullptr)
        (*TGraphMap)["Exogam"][Form("Residue_Graph_EXO_EHG%d",it->first)]->Write();
    }
    File->cd("Exogam");
  }
  for (it2 = DoCalibrationOuter.begin(); it2 != DoCalibrationOuter.end(); it2++) {
    for (it = (it2->second).begin(); it != (it2->second).end(); it++) {
      if (it->second) {
        hnameEXOE = Form("EXO_Outer%d_%d",it2->first,it->first);

        if (!gDirectory->GetDirectory(Form("EXO_Cr%d", it2->first)))
          gDirectory->mkdir(Form("EXO_Cr%d", it2->first));
        gDirectory->cd(Form("EXO_Cr%d", it2->first));

        if (!gDirectory->GetDirectory("Outers"))
          gDirectory->mkdir("Outer");
        gDirectory->cd("Outer");

        (*TH1Map)["Exogam"][hnameEXOE]->Write();
        if((*TGraphMap)["Exogam"][Form("Calib_Graph_%s",hnameEXOE.c_str())]!= nullptr)
          (*TGraphMap)["Exogam"][Form("Calib_Graph_%s",hnameEXOE.c_str())]->Write();
        if((*TGraphMap)["Exogam"][Form("Residue_Graph_%s",hnameEXOE.c_str())]!= nullptr)
          (*TGraphMap)["Exogam"][Form("Residue_Graph_%s",hnameEXOE.c_str())]->Write();
      }
      File->cd("Exogam");
///////////////////////////////////////////////////////////////////////////
  //	tranform an integer to a string
  double fEXO_E(const TExogamData* m_EventData, const unsigned int& i) {
    static string name;
    name = "EXO/E";
    name += NPL::itoa(m_EventData->GetExoCrystal(i));
    return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoE(i),1);
  double fEXO_EHG(const TExogamData* m_EventData, const unsigned int& i) {
    static string name;
    name = "EXO/EHG";
    name += NPL::itoa(m_EventData->GetExoCrystal(i));
    return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoEHG(i),1);
  double fEXO_T(const TExogamData* m_EventData, const unsigned int& i) {
    static string name;
    name = "EXOGAM/Cr_";
    name += NPL::itoa(m_EventData->GetExoCrystal(i));
    return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoTDC(i),1);
  double fEXO_Outer(const TExogamData* m_EventData, const unsigned int& i, const unsigned int OuterNumber) {
    static string name;
    name = "EXO/Outer";
    name += NPL::itoa(m_EventData->GetExoCrystal(i));
    name += "_";
      return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter1(i),1);
      return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter2(i),1);
      return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter3(i),1);
      return CalibrationManager::getInstance()->ApplyCalibration(name, m_EventData->GetExoOuter4(i),1);
      std::cout << "WARNING: Outer number != 0-3, something is wrong\n";
} // namespace EXOGAM_LOCAL

/////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//            Construct Method to be pass to the DetectorFactory              //
////////////////////////////////////////////////////////////////////////////////
NPL::VDetector* TExogamPhysics::Construct() { return (NPL::VDetector*)new TExogamPhysics(); }
NPL::VTreeReader* TExogamPhysics::ConstructReader() { return (NPL::VTreeReader*)new TExogamPhysicsReader(); }

////////////////////////////////////////////////////////////////////////////////
//            Registering the construct method to the factory                 //
////////////////////////////////////////////////////////////////////////////////
  class proxy_exogam {
    public:
      proxy_exogam() {
        NPL::DetectorFactory::getInstance()->AddToken("Exogam", "Exogam");
        NPL::DetectorFactory::getInstance()->AddDetector("Exogam", TExogamPhysics::Construct);
        NPL::DetectorFactory::getInstance()->AddDetectorReader("Exogam", TExogamPhysics::ConstructReader);
      }
  };