ClassTraj.cpp

Example of use of the KVRungeKutta class

Implementation of a classical trajectory calculation for heavy ion collisions

Friction parameter

alpha=240 gives maximum fusion X-section for 40Ca+60Ni@6MeV/u (bfus<=5.25fm: 866mb)

alpha=750 gives 66mb of fusion for 84Kr(605 MeV)+238U

• measured fusion-fission Xsec from PRL 33, 502 (1974): 55+/-15 mb

//Created by KVClassFactory on Thu Jun 24 14:22:17 2010
//Author: John Frankland,,,
 
#ifndef __ClassTraj_H
#define __ClassTraj_H
 
#include "KVRungeKutta.h"
#include "KVWilckeReactionParameters.h"
 
#include "TArc.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TMarker.h"
#include "TTree.h"
#include "TFile.h"
#include "TLatex.h"
 
class ClassTraj : public KVRungeKutta {
   Double_t* coords;
   Double_t* RP;
   Double_t* RT;
   Double_t* VP;
   Double_t* VT;
   Double_t RCM[2];
   Double_t VCM[2];
 
   KVWilckeReactionParameters fWilcke;
   Double_t ELAB, E0;
   Double_t alpha;//friction parameter
 
   //CM trajectory
   TArc* ProjARC;
   TArc* TargARC;
   TGraph* ProjTRAJ;
   Int_t nTRAJPoints;
   TCanvas* CMPositionPad;
   TLatex* infox;
   //lab trajectory
   TArc* ProjARClab;
   TArc* TargARClab;
   TGraph* ProjTRAJlab;
   TCanvas* LabPositionPad;
   Int_t nTRAJPointslab;
   Double_t minDist;
   Double_t EtotDeb, EtotFin, TKEL, TKE, VPOT;
   Double_t ThetaProjLab, absThetaProjLab;
   Double_t EProjLab;
   Double_t ImpactParameter;
   Double_t AngMom, RotationAngle;
   Bool_t fused, inelastic, incoming, interacting, outgoing;
   Double_t tcon, tsep, RotationTime;
   Double_t randomX;
 
   TString sysName;
   TString colName;
 
   TFile* theFile;
   TTree* theTree;
 
   Bool_t fVisualization;//true = show trajectories
 
public:
   ClassTraj(const KVNucleus& proj, const KVNucleus& targ, Bool_t viz = kTRUE);
   virtual ~ClassTraj();
 
   void SetViz(Bool_t on = kTRUE)
   {
      fVisualization = on;
   }
 
   void InitTrajectory(Double_t b, Double_t e, Double_t a = 0);
 
   Double_t dr(Double_t* rp, Double_t* rt) const
   {
      return sqrt(pow(rp[0] - rt[0], 2) + pow(rp[1] - rt[1], 2));
   }
   Double_t dv(Double_t* vp, Double_t* vt) const
   {
      return sqrt(pow(vp[0] - vt[0], 2) + pow(vp[1] - vt[1], 2));
   }
   Double_t dr_x(Double_t* rp, Double_t* rt) const
   {
      // x-component of unit vector from target to projectile
      return (rp[0] - rt[0]) / dr(rp, rt);
   }
   Double_t dr_y(Double_t* rp, Double_t* rt) const
   {
      // y-component of unit vector from target to projectile
      return (rp[1] - rt[1]) / dr(rp, rt);
   }
   Double_t dv_x(Double_t* vp, Double_t* vt) const
   {
      // x-component of unit vector from target to projectile
      return (vp[0] - vt[0]) / dv(vp, vt);
   }
   Double_t dv_y(Double_t* vp, Double_t* vt) const
   {
      // y-component of unit vector from target to projectile
      return (vp[1] - vt[1]) / dv(vp, vt);
   }
   void UpdateARCPositions(Double_t time);
   Double_t ThetaLabProj()
   {
      return TMath::Sign(TMath::RadToDeg() * acos(
                            (VP[0] + VCM[0]) /
                            sqrt(pow(VP[0] + VCM[0], 2) + pow(VP[1] + VCM[1], 2))), VP[1] + VCM[1]);
   };
   Double_t ELabProj()
   {
      return 0.5 * fWilcke.GetAP() * KVNucleus::u() * (pow(VP[0] + VCM[0], 2) + pow(VP[1] + VCM[1], 2));
   };
   Double_t TotalAngularMomentum()
   {
      Double_t lz = dv_x(VP, VT) * dr_y(RP, RT) - dv_y(VP, VT) * dr_x(RP, RT);
      lz *= fWilcke.GetMu() * KVNucleus::u() * dr(RP, RT) * dv(VP, VT) / KVNucleus::hbar;
      return lz;
   }
   Double_t TotalEnergy() const;
 
   virtual void CalcDerivs(Double_t /*X*/, Double_t* /*Y*/, Double_t* /*DY/DX*/);
   void Run(Double_t tmin, Double_t tmax, Double_t dt);
   void Reset()
   {
      // Reset visualization ready for new trajectories calculation
 
      if (CMPositionPad) delete CMPositionPad;
      if (LabPositionPad) delete LabPositionPad;
 
      if (!ProjARC) {
         ProjARC = new TArc(0, 0, fWilcke.GetCP());
         ProjARC->SetFillColor(kBlue);
      }
      if (!TargARC) {
         TargARC = new TArc(0, 0, fWilcke.GetCT());
         TargARC->SetFillColor(kRed);
      }
      if (!ProjARClab) {
         ProjARClab = new TArc(0, 0, 2 * fWilcke.GetCP());
         ProjARClab->SetFillColor(kBlue);
      }
      if (!TargARClab) {
         TargARClab = new TArc(0, 0, 2 * fWilcke.GetCT());
         TargARClab->SetFillColor(kRed);
      }
      CMPositionPad = new TCanvas("CMPOS", "CM Positions", 800, 500);
      CMPositionPad->DrawFrame(-50, -40, 50, 40);
      TMarker* cmark = new TMarker(0, 0, 2);
      cmark->SetMarkerSize(2);
      cmark->Draw();
      TArc* rint = new TArc(0, 0, fWilcke.GetRint());
      rint->SetLineStyle(2);
      rint->SetFillStyle(0);
      rint->Draw();
      if (infox) delete infox;
      infox = new TLatex;
      UpdateInfos(0, 0);
      infox->Draw();
      LabPositionPad = new TCanvas("LABPOS", "Lab Positions", 800, 500);
      LabPositionPad->DrawFrame(-50, -300, 550, 300); //(-100,-80, 100,80);
      cmark = new TMarker(0, 0, 2);
      cmark->SetMarkerSize(2);
      cmark->Draw();
   }
   void UpdateInfos(Double_t b, Double_t time)
   {
      if (infox) {
         if (outgoing)infox->SetText(-40, -30, Form("b=%.1ffm t=%.0ffm/c OUTGOING", b, time));
         else if (incoming)infox->SetText(-40, -30, Form("b=%.1ffm t=%.0ffm/c INCOMING", b, time));
         else if (interacting)infox->SetText(-40, -30, Form("b=%.1ffm t=%.0ffm/c INTERACTING", b, time));
         else infox->SetText(-40, -30, Form("b=%.1ffm t=%.0ffm/c", b, time));
      }
   }
   void CalcTrajectories(Double_t e, Double_t bmin, Double_t bmax, Double_t db, Double_t tmin, Double_t tmax, Double_t dt, Double_t friction = 0);
   void CalcTrajectories(Double_t e, Int_t ntraj, Double_t bmin, Double_t bmax, Double_t tmin, Double_t tmax, Double_t dt, Double_t friction = 0);
 
   void InitTree(Double_t e);
   void DoOneTraj(Double_t tmin, Double_t friction, Double_t b, Double_t e, Double_t tmax, Double_t dt);
   ClassDef(ClassTraj, 0) //Classical trajectory for heavy ion collisions
};
 
#endif

//Created by KVClassFactory on Thu Jun 24 14:22:17 2010
//Author: John Frankland,,,
 
#include "ClassTraj.h"
#include "Riostream.h"
#include "TRandom.h"
using namespace std;
#include "KVTreeAnalyzer.h"
 
ClassImp(ClassTraj)
 
ClassTraj::ClassTraj(const KVNucleus& proj, const KVNucleus& targ, Bool_t viz)
   : KVRungeKutta(8), fWilcke(proj, targ), fVisualization(viz)
{
   // 2-D problem
   // x-axis = beam
   // y-axis = impact parameter
   // projectile comes from -x direction
   alpha = 0;
   coords = new Double_t [8];
   RP = coords;
   RT = coords + 2;
   VP = coords + 4;
   VT = coords + 6;
 
   sysName.Form("%s + %s", proj.GetSymbol(), targ.GetSymbol());
   colName.Form("%s%s", proj.GetSymbol(), targ.GetSymbol());
 
   fWilcke.Print();
 
   CMPositionPad = LabPositionPad = 0;
   ProjARC = TargARC = ProjARClab = TargARClab = 0;
   infox = 0;
   if (fVisualization) Reset();
}
 
void ClassTraj::UpdateARCPositions(Double_t t)
{
   ProjARC->SetX1(RP[0]);
   ProjARC->SetY1(RP[1]);
   TargARC->SetX1(RT[0]);
   TargARC->SetY1(RT[1]);
   ProjTRAJ->SetPoint(nTRAJPoints++, RP[0], RP[1]);
   ProjARClab->SetX1(RP[0] + RCM[0] + VCM[0]*t);
   ProjARClab->SetY1(RP[1] + RCM[1] + VCM[1]*t);
   TargARClab->SetX1(RT[0] + RCM[0] + VCM[0]*t);
   TargARClab->SetY1(RT[1] + RCM[1] + VCM[1]*t);
   ProjTRAJlab->SetPoint(nTRAJPointslab++, RP[0] + RCM[0] + VCM[0]*t, RP[1] + RCM[1] + VCM[1]*t);
}
 
 
ClassTraj::~ClassTraj()
{
   delete [] coords;
}
 
void ClassTraj::InitTrajectory(Double_t b, Double_t e, Double_t fric)
{
   // set up trajectory with impact parameter b
   // energy e mev/nucleon
 
   if (fVisualization)UpdateInfos(b, 0);
   alpha = fric;
   if (fVisualization) {
      ProjTRAJ = new TGraph;
      ProjTRAJ->SetLineWidth(2);
      ProjTRAJlab = new TGraph;
      ProjTRAJlab->SetLineWidth(2);
   }
   ELAB = e;
   Double_t R0 = minDist;
   // lab coordinates
   RP[0] = -R0;
   RP[1] = b;
   RT[0] = RT[1] = 0;
   for (int i = 0; i < 2; i++) {
      RCM[i] = (fWilcke.GetAP() * RP[i] + fWilcke.GetAT() * RT[i]) / (1.*fWilcke.GetAC());
      RP[i] -= RCM[i];
      RT[i] -= RCM[i];
   }
//    cout << "CENTRE OF MASS : (" << RCM[0] << "," << RCM[1] << ")" << endl;
//    cout << "INITIAL DISTANCE :" << dr(RP,RT) << " fm" << endl;
//    cout << "INITIAL POSITIONS :" << endl;
//    cout << "   PROJECTILE : (" << RP[0] << "," << RP[1] << ")" << endl;
//    cout << "   TARGET  : (" << RT[0] << "," << RT[1] << ")" << endl;
//    cout << "POTENTIAL ENERGY : " << fWilcke.GetTotalPotential()->Eval(dr(RP,RT)) << " MeV" << endl;
//    cout << "CM TOTAL ENERGY : " << fWilcke.ECM(ELAB) << " MeV" << endl;
   E0 = fWilcke.ECM(ELAB) - fWilcke.GetTotalPotential()->Eval(dr(RP, RT));
//    cout << "INITIAL CM KINETIC ENERGY : " << E0 << " MeV" << endl;
   Double_t V0 = sqrt(2.*E0 / fWilcke.GetMu() / KVNucleus::kAMU);
//    cout << "INITIAL PROJECTILE VELOCITY : " << V0*KVNucleus::C() << " cm/ns" << endl;
   // lab velocities
   VP[0] = V0;
   VP[1] = VT[0] = VT[1] = 0.;
   for (int i = 0; i < 2; i++) {
      VCM[i] = (fWilcke.GetAP() * VP[i] + fWilcke.GetAT() * VT[i]) / (1.*fWilcke.GetAC());
      VP[i] -= VCM[i];
      VT[i] -= VCM[i];
   }
//    cout << "CENTRE OF MASS VELOCITY : " << VCM[0]*KVNucleus::C() << " cm/ns" << endl;
//    cout << "INITIAL CM VELOCITIES :" << endl;
//    cout << "   PROJECTILE : " << VP[0]*KVNucleus::C() << " cm/ns"  << endl;
//    cout << "   TARGET  : " << VT[0]*KVNucleus::C() << " cm/ns" << endl;
 
   if (fVisualization) {
      nTRAJPoints = 0;
      nTRAJPointslab = 0;
      UpdateARCPositions(0);
      CMPositionPad->cd();
      ProjARC->Draw();
      TargARC->Draw();
      ProjTRAJ->Draw("l");
      CMPositionPad->Modified();
      CMPositionPad->Update();
      LabPositionPad->cd();
      ProjARClab->Draw();
      TargARClab->Draw();
      ProjTRAJlab->Draw("l");
      LabPositionPad->Modified();
      LabPositionPad->Update();
   }
}
 
void ClassTraj::CalcDerivs(Double_t t, Double_t* Y, Double_t* DYDX)
{
   Double_t* rp = Y;
   Double_t* rt = Y + 2;
   Double_t* vp = Y + 4;
   Double_t* vt = Y + 6;
   Double_t* drp = DYDX;
   Double_t* drt = DYDX + 2;
   Double_t* dvp = DYDX + 4;
   Double_t* dvt = DYDX + 6;
 
   Double_t distance = dr(rp, rt);
   Double_t uv_x = dr_x(rp, rt);
   Double_t uv_y = dr_y(rp, rt);
 
   Double_t Force = -fWilcke.GetTotalPotential()->Derivative(distance);
   if (TMath::IsNaN(Force)) Force = 0;
 
   for (int i = 0; i < 2; i++) {
      drp[i] = vp[i];
      drt[i] = vt[i];
   }
   Double_t fp = Force / (fWilcke.GetAP() * KVNucleus::kAMU);
   dvp[0] = fp * uv_x;
   dvp[1] = fp * uv_y;
   Double_t ft = -Force / (fWilcke.GetAT() * KVNucleus::kAMU);
   dvt[0] = ft * uv_x;
   dvt[1] = ft * uv_y;
 
   if (fInitialDeriv) {
      if (distance <= minDist)minDist = distance;
      if (incoming && distance < fWilcke.GetRint() + 10.) {
         incoming = kFALSE;
         interacting = kTRUE;
      }
      else if (interacting && distance > fWilcke.GetRint() + 10.) {
         interacting = kFALSE;
         outgoing = kTRUE;
      }
      if (tcon < 0 && distance < fWilcke.GetRint()) tcon = t;
      if (tcon > 0 && tsep < 0) {
         if (t > tcon && distance > fWilcke.GetRint()) tsep = t;
      }
   }
   // friction force
   if (alpha > 0) {
      if (fInitialDeriv) {
         // draw random number with mean=0 and sigma=1
         randomX = gRandom->Gaus();
      }
      if (distance < fWilcke.GetRint()) {
         Double_t vrel = dv(vp, vt);
         Double_t vv_x = dv_x(vp, vt);
         Double_t vv_y = dv_y(vp, vt);
         Double_t vpar = vrel * (vv_x * uv_x + vv_y * uv_y);
         Double_t FRICTION = -alpha * vpar;
         Double_t fforce = FRICTION * randomX + FRICTION;
         fp = fforce / (fWilcke.GetAP() * KVNucleus::kAMU);
         dvp[0] += fp * (uv_x);
         dvp[1] += fp * (uv_y);
         ft = -fforce / (fWilcke.GetAT() * KVNucleus::kAMU);
         dvt[0] += ft * (uv_x);
         dvt[1] += ft * (uv_y);
      }
   }
}
 
void ClassTraj::Run(Double_t tmin, Double_t tmax, Double_t dt)
{
   // Perform the simulation from time 'tmin' to 'tmax'.
 
   Double_t Time = tmin;
   tcon = tsep = -1.;
   incoming = kTRUE;
   interacting = kFALSE;
   outgoing = kFALSE;
   while (Time < tmax) {
      //cout << "ETOT = " << TotalEnergy() << endl;
      // take Runge-Kutta step from Time to Time+dt
 
      if (outgoing) {
         dt *= 1.5;
         if (Time + dt > tmax) {
            dt = tmax - Time;
         }
      }
      Integrate(coords, Time, Time + dt, dt);
 
      Time += dt;
      if (fVisualization) {
         UpdateInfos(ImpactParameter, Time);
         UpdateARCPositions(Time);
         CMPositionPad->Modified();
         CMPositionPad->Update();
         LabPositionPad->Modified();
         LabPositionPad->Update();
      }
   }
 
   if (fVisualization) {
      ProjTRAJ->SetLineWidth(1);
      ProjTRAJlab->SetLineWidth(1);
      ProjTRAJ->SetLineStyle(2);
      ProjTRAJlab->SetLineStyle(2);
   }
}
 
Double_t ClassTraj::TotalEnergy() const
{
   Double_t ekin = 0.5 * KVNucleus::u() * (
                      fWilcke.GetAP() * (VP[0] * VP[0] + VP[1] * VP[1])
                      + fWilcke.GetAT() * (VT[0] * VT[0] + VT[1] * VT[1]));
   Double_t epot = fWilcke.GetTotalPotential()->Eval(dr(RP, RT));
   return ekin + epot;
}
 
void ClassTraj::CalcTrajectories(Double_t e, Double_t bmin, Double_t bmax, Double_t db, Double_t tmin, Double_t tmax, Double_t dt, Double_t friction)
{
   alpha = friction;
   InitTree(e);
 
   for (Double_t b = bmin; b <= bmax; b += db) {
 
      DoOneTraj(tmin, friction, b, e, tmax, dt);
   }
   theFile->Write();
   theFile->Close();
}
 
void ClassTraj::InitTree(Double_t e)
{
   TString simTitle = Form("%s %f MeV/u (alpha=%f)", sysName.Data(), e, alpha);
   TString filename = Form("ClassTraj_%s_%fAMeV_alpha%f.root", colName.Data(), e, alpha);
   theFile = new TFile(filename, "recreate");
   theTree = new TTree("ClassTrajSim", simTitle);
   theTree->Branch("ImpactParameter", &ImpactParameter);
   theTree->Branch("AngMom", &AngMom);
   theTree->Branch("TKEL", &TKEL);
   theTree->Branch("TKE", &TKE);
   theTree->Branch("VPOT", &VPOT);
   theTree->Branch("minDist", &minDist);
   theTree->Branch("ThetaProjLab", &ThetaProjLab);
   theTree->Branch("absThetaProjLab", &absThetaProjLab);
   theTree->Branch("EProjLab", &EProjLab);
   theTree->Branch("fused", &fused);
   theTree->Branch("inelastic", &inelastic);
   theTree->Branch("RotationTime", &RotationTime);
   theTree->Branch("RotationAngle", &RotationAngle);
   theTree->Branch("tcon", &tcon);
   theTree->Branch("tsep", &tsep);
}
 
void ClassTraj::DoOneTraj(Double_t tmin, Double_t friction, Double_t b, Double_t e, Double_t tmax, Double_t dt)
{
   minDist = 500;
   InitTrajectory(b, e, friction);
 
   ImpactParameter = b;
   AngMom = fWilcke.k(e) * b;
 
   Run(tmin, tmax, dt);
 
   ThetaProjLab = ThetaLabProj();
   EProjLab = ELabProj();
   absThetaProjLab = abs(ThetaProjLab);
   TKE = TotalEnergy();
   TKEL = fWilcke.ECM(ELAB) - TKE;
   RotationTime = (tsep > 0 ? tsep - tcon : tmax - tcon);
   VPOT = fWilcke.GetTotalPotential()->Eval(dr(RP, RT));
   inelastic = fused = kFALSE;
 
   if (dr(RP, RT) < fWilcke.GetRint()) fused = kTRUE;
   if (TKEL > 1.e-4 && !fused) inelastic = kTRUE;
 
   //cout << "Trot="<<RotationTime<<"fm/c L="<<AngMom<<"hbar Rmin="<<minDist<<"fm"<<endl;
   Double_t Period = 2 * TMath::Pi() * fWilcke.GetMu() * KVNucleus::kAMU * minDist * minDist / (AngMom * KVNucleus::hbar);
   RotationAngle = 2 * TMath::Pi() * RotationTime / Period * TMath::RadToDeg();
   //cout<<"Period = "<<Period<<"fm/c  --  Angle turned = "<<RotationAngle<<"deg."<<endl;
 
   theTree->Fill();
}
 
void ClassTraj::CalcTrajectories(Double_t e, Int_t ntraj, Double_t bmin, Double_t bmax, Double_t tmin, Double_t tmax, Double_t dt, Double_t friction)
{
   alpha = friction;
   InitTree(e);
 
   while (ntraj--) {
 
      if (!(ntraj % 5000)) cout << ntraj << " more to go..." << endl;
 
      Double_t b = sqrt(gRandom->Uniform(bmin * bmin, bmax * bmax));
      DoOneTraj(tmin, friction, b, e, tmax, dt);
   }
   theFile->Write();
//    theFile->Close();
   theTree->StartViewer();
}