KVElasticScatter.cpp

/*
$Id: KVElasticScatter.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 Thu Oct 19 22:31:02 2006
//Author: John Frankland
 
#include "KVElasticScatter.h"
#include "KVMultiDetArray.h"
#include "TH1F.h"
#include "KVDetector.h"
#include "KVGroup.h"
#include "KVTarget.h"
#include "KV2Body.h"
#include "TObjArray.h"
#include "KVExpDB.h"
#include "TRandom.h"
using namespace std;
 
ClassImp(KVElasticScatter)
 
 
 
KVElasticScatter::KVElasticScatter(Int_t run)
{
   //Default constructor
   if (!gMultiDetArray)
      Error("KVElasticScatter", "First define the multidetector array for the simulation");
   SetRun(run);
   gMultiDetArray->SetSimMode(kTRUE);
}
 
 
 
 
 
KVElasticScatter::~KVElasticScatter()
{
   //Destructor
   if (fDepth)
      delete fDepth;
   if (fTheta)
      delete fTheta;
   gMultiDetArray->SetSimMode(kFALSE);
}
 
 
 
 
 
void KVElasticScatter::SetRun(Int_t run)
{
   //Set beam, detector parameters, target, etc. for run
   //
   // The default reaction kinematics is that of the KVDBSystem associated to the run.
   //
   // The default scattering position in the target is random.
   gMultiDetArray->SetParameters(run);
   fTarget = gMultiDetArray->GetTarget();
   fTarget->SetRandomized();
   fAtomicDensity = fTarget->GetAtomsPerCM2() * 1.e-24; //in barn^-1 units
   fIntLayer = 0;
   fMultiLayer = (fTarget->NumberOfLayers() > 1);
   // get reaction propoerties from system
   auto sys = gExpDB->GetDBRun(run)->GetSystem();
   fKinematics = sys->GetKinematics();
}
 
 
 
 
 
void KVElasticScatter::SetDetector(const Char_t* det)
{
   // Give name of detector towards which elastic scattering will occur.
   //
   // \note this does not necessarily mean that the particle will actually
   // arrive in the given detector (it may stop before it), but the aim is to
   // define the trajectory of the particle in the array.
 
   fDetector = gMultiDetArray->GetDetector(det);
   if (!fDetector) {
      Error("SetDetector", "Detector %s is unknown!", det);
      return;
   }
   // get trajectory for detector
   if (fDetector->GetNode()->GetNTraj() > 1) {
      Warning("SetDetector", "Ambiguous trajectory choice: detector %s has %d trajectories", det, fDetector->GetNode()->GetNTraj());
   }
   auto traj = (KVGeoDNTrajectory*)fDetector->GetNode()->GetTrajectories()->First();
   Info("SetDetector", "Using trajectory: %s", traj->GetName());
 
   //we store the association between detector type and index in list
   fDetInd.Clear();
   fAlignedDetectors.Clear();
 
   traj->IterateBackFrom();
 
   Int_t i = 0;
   KVGeoDetectorNode* gdn = nullptr;
   while ((gdn = traj->GetNextNode())) {
      KVDetector* d = gdn->GetDetector();
      fAlignedDetectors.Add(d);
      fDetInd.SetValue(d->GetLabel(), i);
      ++i;
   }
   //store number of aligned detectors
   fNDets = i;
   fHistos.Expand(fNDets);
   fHistos.SetOwner();      //will delete histos it stores
}
 
 
 
 
 
void KVElasticScatter::SetTargetScatteringLayer(const Char_t* name)
{
   //For multilayer targets, use this method to choose in which
   //layer the scattering will take place.
   //
   //If name="", reset any previous choice so that scattering
   //can take place in any layer
   if (!fTarget) {
      cout <<
           "<KVElasticScatter::SetTargetScatteringLayer> : No target set. Set run first."
           << endl;
      return;
   }
   fIntLayer = fTarget->GetLayerIndex(name);
   if (fIntLayer)
      fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
}
 
 
 
 
 
void KVElasticScatter::SetEbinning(Int_t nbins)
{
   //Set binning of the GetEnergy histogram
   // Default value is 500
   fBinE = nbins;
}
 
 
 
 
bool KVElasticScatter::initial_checks_and_reset()
{
   if (!fDetector) {
      cout <<
           "<KVElasticScatter::CalculateScattering> : Set detector first" <<
           endl;
      return false;
   }
   fHistos.Clear();         //delete any previous histograms
 
   KVDetector* d;
   TIter n(&fAlignedDetectors);
   while ((d = (KVDetector*) n()))
      fHistos.Add(new TH1F(Form("hEloss_%s", d->GetName()), "Eloss (MeV)", fBinE, 0., 0.));
 
   return true;
}
 
 
 
 
std::pair<double, double> KVElasticScatter::get_random_angles_for_scattering(const KV2Body& scattering_kinematics)
{
   auto max_scattering_angle = scattering_kinematics.GetMaxAngleLab(fNucleusOfInterest);
   if (fDetector->GetEntranceWindow().GetThetaMin() > max_scattering_angle) {
      Warning("CalculateScattering",
              "Detector %s has minimum theta %f deg. > max scattering angle of %s, %f deg.: abandon calculation",
              fDetector->GetName(), fDetector->GetEntranceWindow().GetThetaMin(), scattering_kinematics.GetNucleus(fNucleusOfInterest)->GetSymbol(),
              max_scattering_angle);
      return {-1., -1.};
   }
   //set random direction of outgoing projectile
   double th, ph;
   th = ph = 0.;
   do {
      fDetector->GetRandomAngles(th, ph);
   }
   while (th > max_scattering_angle);
 
   return {th, ph};
}
 
 
 
 
void KVElasticScatter::detect_particle_fill_histograms(KVNucleus* ejectile, double theta, double phi, double xsec)
{
   int j = 0;
   TIter n(&fAlignedDetectors);
   KVDetector* d;
   while ((d = (KVDetector*) n())) {
      auto cew = d->GetCentreOfEntranceWindow();
      d->DetectParticle(ejectile, &cew);
      //fill histograms
      if(fFWHM) // take into account detector resolution
      {
         auto E = d->GetEnergy();
         ((TH1F*)fHistos[j++])->Fill(gRandom->Gaus(E, get_det_resolution(E)/2.35), xsec * sin(theta * TMath::DegToRad()));
      }
      else
         ((TH1F*)fHistos[j++])->Fill(d->GetEnergy(), xsec * sin(theta * TMath::DegToRad()));
      //prepare for next round: set energy loss to zero
      d->Clear();
   }
}
 
 
 
 
void KVElasticScatter::CalculateScattering(Int_t N)
{
   // Perform scattering 'N' times
 
   if (!initial_checks_and_reset())
      return;
 
   fNtirages = N;
 
   //set up histograms
   // -- histograms for debugging: depth in target and angular distribution
   if (fDepth)
      delete fDepth;
   if (fTheta)
      delete fTheta;
   fDepth =
      new TH1F("hDepth", "Depth (mg/cm2)", 500, 0., fTarget->GetTotalEffectiveThickness());
   fTheta = new TH1F("hTheta", "Theta (deg.)", 500, 0., 0.);
 
   //set random interaction point for scattering
   TVector3 IP;
   if (fIntLayer) {
      fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
      IP = fTarget->GetInteractionPoint();
   }
   else {
      IP = fTarget->GetInteractionPoint(fBeamDirection);
   }
 
   /* -------------------------------------------------------------------------------------------------------------------------- */
 
   if (fGetTargFromLay) {
      // get target nucleus properties from scattering layer ?
      KVMaterial* targ_mat = fTarget->GetLayer(IP);
      auto zt = targ_mat->GetZ();
      auto at = targ_mat->GetMass();
      // check for fractional Z or A - compound material!
      if ((zt != int(zt)) || at != int(at)) {
         Error("CalculateScattering", "Target has non-elemental layer with effective Z=%f and A=%f. Cannot use.",
               zt, at);
         fGetTargFromLay = false;
      }
      else
         fKinematics->SetTarget({int(zt), int(at)});
   }
 
   for (int i = 0; i < N; i++) {
 
      reset_before_new_scattering();
 
      auto Projectile = *(fKinematics->GetNucleus(KV2Body::projectile));
      auto Target = *(fKinematics->GetNucleus(KV2Body::target));
 
      //calculate slowing of incoming projectile
      fTarget->SetIncoming();
      fTarget->DetectParticle(&Projectile);
 
      KV2Body scattering{Projectile, Target, fExx};
      if (fOutgoingPL)
         scattering.SetOutgoing(*fOutgoingPL);
 
      if (fExx > 0)
         scattering.SetEDiss(fExx);
 
      scattering.CalculateKinematics();
 
      auto theta_phi = get_random_angles_for_scattering(scattering);
 
      //set energy of scattered nucleus
      auto elab = scattering.GetELab(fNucleusOfInterest, theta_phi.first);
      auto ejectile = scattering.GetNucleus(fNucleusOfInterest);
 
      auto kine_branch = fKineSol;
      if (fKineSol == KV2Body::low_E_branch) {
         if (!elab->second()) {
            Warning("CalculateScattering", "Low energy branch requested, but only high energy kinematics for theta = %f deg.",
                    theta_phi.first);
            ejectile->SetEnergy(*elab->first());
            kine_branch = KV2Body::high_E_branch;
         }
         else {
            ejectile->SetEnergy(*elab->second());
         }
      }
      else if (fKineSol == KV2Body::both_branches) {
         if (!elab->second()) {
            Warning("CalculateScattering", "Both branches requested, but only high energy kinematics for theta = %f deg.",
                    theta_phi.first);
            ejectile->SetEnergy(*elab->first());
            kine_branch = KV2Body::high_E_branch;
         }
         else {
            ejectile->SetEnergy(1.0);
            ejectile->SetTheta(theta_phi.first);
            auto angle_with_beam = TMath::RadToDeg() * ejectile->Angle(fBeamDirection);
            auto xsec_hi = TMath::Abs(scattering.GetXSecRuthLab(angle_with_beam, fNucleusOfInterest, KV2Body::high_E_branch));
            auto xsec_lo = TMath::Abs(scattering.GetXSecRuthLab(angle_with_beam, fNucleusOfInterest, KV2Body::low_E_branch));
            if (gRandom->Uniform() < xsec_hi / (xsec_hi + xsec_lo)) {
               // high
               ejectile->SetEnergy(*elab->first());
               kine_branch = KV2Body::high_E_branch;
            }
            else {
               // low
               ejectile->SetEnergy(*elab->second());
               kine_branch = KV2Body::low_E_branch;
            }
         }
      }
      else
         ejectile->SetEnergy(*elab->first());
 
      ejectile->SetTheta(theta_phi.first);
      ejectile->SetPhi(theta_phi.second);
 
      // get Rutherford Xsec for projectile scattered to theta degrees in lab
      // taking into account the initial direction of the beam (if beam direction = (0,0,1) this is just the same as theta)
      double angle_with_beam = TMath::RadToDeg() * ejectile->Angle(fBeamDirection);
      auto xsec = TMath::Abs(scattering.GetXSecRuthLab(angle_with_beam, fNucleusOfInterest, kine_branch));
 
      // fill angular distribution
      fTheta->Fill(theta_phi.first, xsec);
 
      //slowing of outgoing ejectile in target
      fTarget->SetOutgoing();
 
      fTarget->DetectParticle(ejectile);
      fDepth->Fill(IP.z());
 
      //now detect particle in detector(s)
      detect_particle_fill_histograms(ejectile, theta_phi.first, theta_phi.second, xsec);
 
      //set random interaction point for scattering
      if (fIntLayer) {
         fTarget->SetInteractionLayer(fIntLayer, fBeamDirection);
         IP = fTarget->GetInteractionPoint();
      }
      else {
         IP = fTarget->GetInteractionPoint(fBeamDirection);
      }
      if (fGetTargFromLay) {
         // get target nucleus properties from scattering layer ?
         KVMaterial* targ_mat = fTarget->GetLayer(IP);
         auto zt = targ_mat->GetZ();
         auto at = targ_mat->GetMass();
         // check for fractional Z or A - compound material!
         if ((zt != int(zt)) || at != int(at)) {
            Error("CalculateScattering", "Target has non-elemental layer with effective Z=%f and A=%f. Cannot use.",
                  zt, at);
            fGetTargFromLay = false;
         }
         else
            fKinematics->SetTarget({int(zt), int(at)});
      }
   }
   end_of_run();
}
 
 
 
 
 
 
TH1F* KVElasticScatter::GetEnergy(const Char_t* label)
{
   //Energy loss in detector with given 'label' through which scattered particle passes.
 
   return (fDetInd.HasParameter(label) ? GetEnergy(fDetInd.GetIntValue(label)) : 0);
}
 
 
 
 
 
TH1F* KVElasticScatter::GetEnergy(Int_t index)
{
   //Energy loss in any detector through which scattered particle passes.
   //
   //The index corresponds to the order in which detectors are crossed by the
   //particle, beginning with 0 for the first detector, and ending with
   //(GetNDets()-1)
   return (TH1F*)fHistos[index];
}
 
 
 
 
void KVElasticScatter::Print() const
{
   // Print details of calculation to be performed
 
   Info("Print", "CALCULATION DETAILS:");
   if (fGetTargFromLay)
      Info("Print", "Projectile is %s %.1f MeV/u; target nucleus depends on scattering layer.",
           fKinematics->GetNucleus(KV2Body::projectile)->GetSymbol(),
           fKinematics->GetNucleus(KV2Body::projectile)->GetAMeV());
   else
      Info("Print", "Entrance channel is: %s", fKinematics->GetTitle());
   if (fExx > 0) {
      Info("Print", "Inelastic scattering: outgoing Ex=%f", fExx);
   }
   Info("Print", "Target used is:");
   fTarget->Print();
   Info("Print", "Scattering will be calculated for outgoing %s",
        fNucleusOfInterest == KV2Body::projectile_like ? "PROJECTILE_LIKE" : "TARGET_LIKE");
   Info("Print", "We will use %s kinematic solutions",
        fKineSol == KV2Body::both_branches ? "BOTH" : (fKineSol == KV2Body::high_E_branch ? "the HIGH ENERGY" : "the LOW ENERGY"));
   if (fDetector) {
      Info("Print", "Particles scattered towards %s: theta:[%f,%f] phi:[%f,%f]", fDetector->GetName(),
           fDetector->GetThetaMin(), fDetector->GetThetaMax(), fDetector->GetPhiMin(), fDetector->GetPhiMax());
   }
}
 
 
//__________________________________________________________________//