KVElasticCountRates.cpp

/*
$Id: KVElasticCountRates.cpp,v 1.9 2007/04/04 10:39:17 ebonnet Exp $
$Revision: 1.9 $
$Date: 2007/04/04 10:39:17 $
*/
 
//Created by KVClassFactory on Fri Nov 20 2015
//Author: John Frankland
 
#include "KVElasticCountRates.h"
#include "KVMultiDetArray.h"
#include "TH1F.h"
#include "KVDetector.h"
#include "KVGroup.h"
#include "KVTarget.h"
#include "KV2Body.h"
#include "KVNucleus.h"
#include "TObjArray.h"
#include "KVDatime.h"
#include <TH2F.h>
#include <vector>
#include <algorithm>
 
using namespace std;
 
ClassImp(KVElasticCountRate)
ClassImp(KVElasticCountRates)
 
 
 
 
void KVElasticCountRates::SetAngularRange(Double_t theta_min, Double_t theta_max, Double_t phi_min, Double_t phi_max)
{
   fAngularRange = {theta_min, theta_max, phi_min, phi_max};
   fVolume = (theta_max - theta_min) * (phi_max - phi_min) * TMath::DegToRad() * TMath::DegToRad();
}
 
 
 
 
bool KVElasticCountRates::initial_checks_and_reset()
{
   fArrayHistos.Delete();
   fRates.clear();
   fGlobalMap = new TH2F("hGlobalMap", "X-section in x(vert)-y(hori) plane",
                         500, -fAngularRange.GetThetaMax(), fAngularRange.GetThetaMax(),
                         500, -fAngularRange.GetThetaMax(), fAngularRange.GetThetaMax());
   return true;
}
 
 
 
 
void KVElasticCountRates::reset_before_new_scattering()
{
   det_sim.ClearHitGroups();
}
 
 
 
 
std::pair<double, double> KVElasticCountRates::get_random_angles_for_scattering(const KV2Body&)
{
   //set random direction of outgoing projectile
   double theta,phi;
   fAngularRange.GetRandomAngles(theta, phi, "random");
   return {theta,phi};
}
 
 
 
 
void KVElasticCountRates::detect_particle_fill_histograms(KVNucleus* ejectile, double theta, double phi, double xsec)
{
   //now detect particle in array
   auto detectors = det_sim.PropagateParticle(ejectile);
   FillHistograms(&detectors, theta, phi, xsec);
}
 
 
 
 
void KVElasticCountRates::end_of_run()
{
   PrintResults(fBeamIntensity);
}
 
 
 
 
void KVElasticCountRates::FillHistograms(const KVNameValueList* dets, double theta, double phi, double xsec)
{
   // parse the list dets
   // fill histograms with energy loss for all detectors
   // clear the detector energy losses
   // delete the list
 
   if (!dets) return;
 
   Int_t ndets = dets->GetNpar();
   for (int i = 0; i < ndets; i++) {
      TString detname = dets->GetNameAt(i);
      KVDetector* det = gMultiDetArray->GetDetector(detname);
      if (!det) continue;
      TH1F* histo = (TH1F*)fArrayHistos.FindObject(detname);
      if (!histo) {
         histo = new TH1F(detname, Form("Eloss in %s", detname.Data()), GetEbinning(), 0, 0);
         fArrayHistos.Add(histo);
      }
      double de = dets->GetDoubleValue(i);
      histo->Fill(de, xsec * sin(theta * TMath::DegToRad()));
      histo = (TH1F*)fArrayHistos.FindObject(detname + "_dW");
      if (!histo) {
         histo = new TH1F(detname + "_dW", Form("Solid angle of %s", detname.Data()), GetEbinning(), 0, 0);
         fArrayHistos.Add(histo);
      }
      histo->Fill(de, sin(theta * TMath::DegToRad()));
      TH2F* histo2 = (TH2F*)fArrayHistos.FindObject(detname + "_map");
      if (!histo2) {
         histo2 = new TH2F(detname + "_map", Form("Map of %s", detname.Data()), 100, 0, 0, 100, 0, 0);
         fArrayHistos.Add(histo2);
      }
      histo2->Fill(theta, phi, xsec);
      fGlobalMap->Fill(theta * sin(phi * TMath::DegToRad()), theta * cos(phi * TMath::DegToRad()), xsec);
      det->Clear();
   }
}
 
 
 
struct count_rate {
   TString detector;
   double counts;
   double energy;
   double theta;
   double phi;
   double fluence;
   double dissipation;
   double intXsec;
   count_rate(TString n, double c, double e, double t, double p, double f, double d, double i)
      : detector(n), counts(c), energy(e), theta(t), phi(p), fluence(f), dissipation(d), intXsec(i) {}
   void print()
   {
      printf("%s \t:  N=%8.2f/sec. \t <E>=%7.1f MeV \t Tot.Xsec=%9.3E barn \t fluence=%9.3E/sec./cm**2 \t dissip.=%9.3E MeV/sec./cm**2\n",
             detector.Data(), counts, energy, intXsec, fluence, dissipation);
   }
   void PutInList(KVNameValueList& nl)
   {
      if (nl.HasParameter("DetList")) {
         KVString tmp = nl.GetStringValue("DetList");
         tmp += ",";
         tmp += detector;
         nl.SetValue("DetList", tmp.Data());
      }
      else {
         nl.SetValue("DetList", detector.Data());
      }
      nl.SetValue(Form("%s.Counts(/sec)", detector.Data()), Form("%8.2f", counts));
 
      nl.SetValue(Form("%s.Counts(/sec)", detector.Data()), Form("%8.2f", counts));
      nl.SetValue(Form("%s.Emean(MeV)", detector.Data()), Form("%7.1f", energy));
      nl.SetValue(Form("%s.TotXsec(barn)", detector.Data()), Form("%9.3E", intXsec));
      nl.SetValue(Form("%s.Fluence(/sec./cm**2)", detector.Data()), Form("%9.3E", fluence));
      nl.SetValue(Form("%s.dissip(MeV/sec./cm**2)", detector.Data()), Form("%9.3E", dissipation));
   }
};
 
 
 
bool compare_count_rates(count_rate a, count_rate b)
{
   if (TMath::Abs(a.theta - b.theta) < 0.5) return a.phi < b.phi;
   return a.theta < b.theta;
}
 
 
 
 
void KVElasticCountRates::PrintResults(Double_t beam_intensity)
{
   // Print mean energy deposit & counting rate for given beam intensity in particles per second
 
   TIter it(&fArrayHistos);
   TH1F* h;
   fRates.clear();
 
   std::vector<count_rate> count_rates;
 
   while ((h = (TH1F*)it())) {
      TString name = h->GetName();
      if (!name.EndsWith("_dW") && !name.EndsWith("_map")) {
         TH2F* map = (TH2F*)fArrayHistos.FindObject(name + "_map");
         // integrated cross-section
         double intXsec = h->Integral() * fVolume / fNtirages;
         // counting rate
         double rate = fAtomicDensity * beam_intensity * intXsec;
         // mean energy
         double emean = h->GetMean();
         KVDetector* det = gMultiDetArray->GetDetector(name);
         double fluence = rate / det->GetEntranceWindowSurfaceArea();
         double dissipation = emean * rate / det->GetEntranceWindowSurfaceArea();
         count_rates.push_back(
            count_rate(name, rate, emean, map->GetMean(), map->GetMean(2), fluence, dissipation, intXsec)
         );
         fRates[name.Data()] = KVElasticCountRate(rate, emean, intXsec, fluence, dissipation);
      }
   }
   std::sort(count_rates.begin(), count_rates.end(), compare_count_rates);
 
   for (std::vector<count_rate>::iterator it = count_rates.begin(); it != count_rates.end(); ++it) {
      it->print();
   }
}
 
 
 
 
KVNameValueList KVElasticCountRates::PutResultsInList(Double_t beam_intensity)
{
   // Print mean energy deposit & counting rate for given beam intensity in particles per second
   KVNameValueList nl;
   nl.SetName("Generated by KVElasticCountRates::PutResultsInList method");
   KVDatime dt;
   nl.SetTitle(dt.AsSQLString());
   TIter it(&fArrayHistos);
   TH1F* h;
   fRates.clear();
   nl.SetValue("Intensity(pps)", beam_intensity);
 
   std::vector<count_rate> count_rates;
 
   while ((h = (TH1F*)it())) {
      TString name = h->GetName();
      if (!name.EndsWith("_dW") && !name.EndsWith("_map")) {
         TH2F* map = (TH2F*)fArrayHistos.FindObject(name + "_map");
         // integrated cross-section
         double intXsec = h->Integral() * fVolume / fNtirages;
         // counting rate
         double rate = fAtomicDensity * beam_intensity * intXsec;
         // mean energy
         double emean = h->GetMean();
         KVDetector* det = gMultiDetArray->GetDetector(name);
         double fluence = rate / det->GetEntranceWindowSurfaceArea();
         double dissipation = emean * rate / det->GetEntranceWindowSurfaceArea();
         count_rates.push_back(
            count_rate(name, rate, emean, map->GetMean(), map->GetMean(2), fluence, dissipation, intXsec)
         );
         fRates[name.Data()] = KVElasticCountRate(rate, emean, intXsec, fluence, dissipation);
      }
   }
   std::sort(count_rates.begin(), count_rates.end(), compare_count_rates);
 
   for (std::vector<count_rate>::iterator it = count_rates.begin(); it != count_rates.end(); ++it) {
      it->PutInList(nl);
   }
   return nl;
}
 
 
//__________________________________________________________________//