GRUTinizer
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TReaction.h
Go to the documentation of this file.
1 #ifndef TREACTION_H
2 #define TREACTION_H
3 
4 #include <iostream>
5 #include <string>
6 #include <cmath>
7 
8 #include "TNamed.h"
9 #include "TMath.h"
10 #include "TGraph.h"
11 
12 #include "TNucleus.h"
13 
14 #ifndef PI
15 #define PI (3.14159265358979312e+00)
16 #endif
17 
18 #ifndef R2D
19 #define R2D (5.72957795130823229e+01)
20 #endif
21 
22 #ifndef D2R
23 #define D2R (1.74532925199432955e-02)
24 #endif
25 
26 // *********** *********** *********** *********** *********** *********** *********** //
27 //
28 // __TReaction__
29 //
30 // A simple, intuitive two body reaction class.
31 // Lab frame variables are generally calculated by boosting from the CM frame
32 // Reactions have the form A(B,C)D where;
33 // [0] - B is "beam"
34 // [1] - A is "targ"
35 // [2] - C is "ejec"
36 // [3] - D is "reco"
37 //
38 // Uses planck units (c=1 and is dimensionless)
39 // M = mc^2
40 // E = full mass energy (inc. kinetic energy)
41 // T = kinetic energy
42 // V = magnitude of velocity (=beta)
43 // P = magnitude of momentum
44 // G = gamma factor (=[1-V^2]-0.5)
45 //
46 // * An excitation energy can be included in the final state heavy recoil nucleus using ex3
47 // in the reaction initialisation or by using SetExcEnergy(exc)
48 // * ALL ENERGIES ARE IN MEV - MASSES TOO !
49 // * currently only works for part==2, I think that part==3 just means theta_cm -> -theta_cm
50 //
51 // BUGS
52 // THERE ARE PERSISTANT NUMERICAL ERRORS (~0.1-1%) ...??! WHY?!!?!
53 // *********** *********** *********** *********** *********** *********** *********** //
54 
55 
56 
57 class TReaction : public TNamed {
58 public:
59  TReaction(const char *beam, const char *targ, const char *ejec, const char *reco, double ebeam=0.0, double ex3=0.0, bool inverse=false);
60 
61  void InitReaction();
62 
63  // returns reaction input parameters
64  TNucleus *GetNucleus(int part) { return fNuc[part];}
65  double GetM(int part) { return fM[part]; }
66  double GetExc() { return fExc; }
67  double GetQVal() { return fQVal; }
68  bool Inverse() { return fInverse; }
69  double GetTBeam(bool inverse);
70  double GetVBeam() { return fVLab[0]; }
71 
72  // CM frame properties
73  double GetInvariantMass() { return fInvariantMass; }
74  double GetCmE() { return fCmE; }
75  double GetCmTi() { return fCmTi; }
76  double GetCmTf() { return fCmTf; }
77  double GetCmV() { return fCmV; }
78  double GetCmP() { return fCmP; }
79  double GetCmG() { return fCmG; }
80 
81  // particle properties in CM frame
82  double GetECm(int part) { return fECm[part]; }
83  double GetTCm(int part) { return fTCm[part]; }
84  double GetVCm(int part) { return fVCm[part]; }
85  double GetPCm(int part) { return fPCm[part]; }
86  double GetGCm(int part) { return fGCm[part]; }
87 
88  // particle properties in LAB frame (default is beam)
89  double GetThetaMax(int part) {return fThetaMax[part]; }
90  // this stuff depends on the lab angle as the cm motion and the particle in the cm are coupled
91  double GetELab(double theta_lab=0.0, int part=0)
92  { return GetELabFromThetaCm(ConvertThetaLabToCm(theta_lab,part),part); }
93  double GetTLab(double theta_lab=0.0, int part=0)
94  { return GetTLabFromThetaCm(ConvertThetaLabToCm(theta_lab,part),part); }
95  double GetVLab(double theta_lab=0.0, int part=0)
96  { return GetVLabFromThetaCm(ConvertThetaLabToCm(theta_lab,part),part); }
97  double GetPLab(double theta_lab=0.0, int part=0)
98  { return GetPLabFromThetaCm(ConvertThetaLabToCm(theta_lab,part),part); }
99  double GetGLab(double theta_lab=0.0, int part=0)
100  { return GetGLabFromThetaCm(ConvertThetaLabToCm(theta_lab,part),part); }
101 
102  // this stuff depends on the CM angle as the cm motion and the particle in the cm are coupled
103  double GetELabFromThetaCm(double theta_cm=0.0, int part=0); // FULL MASS+KINETIC ENERGY
104  double GetTLabFromThetaCm(double theta_cm=0.0, int part=0); // KINETIC ENERGY
105  double GetVLabFromThetaCm(double theta_cm=0.0, int part=0);
106  double GetPLabFromThetaCm(double theta_cm=0.0, int part=0);
107  double GetGLabFromThetaCm(double theta_cm=0.0, int part=0);
108 
109  double GetExcEnergy(double ekin=0.00, double theta_lab=0.00, int part = 2);
110  void AnalysisAngDist(double ekin, double theta_lab, int part, double &exc, double &theta_cm, double &omega_lab2cm);
111  double AnalysisBeta(double ekin, int part);
112 
113  double GetRutherfordCm(double theta_cm, int part = 2, bool Units_mb = true);
114  double GetRutherfordLab(double theta_lab, int part = 2, bool Units_mb = true);
115 
116  // Conversion from LAB frame to CM frame
117  double ConvertThetaLabToCm(double theta_lab, int part = 2);
118  double ConvertOmegaLabToCm(double theta_lab, int part = 2);
119  void ConvertLabToCm(double theta_lab, double omega_lab, double &theta_cm, double &omega_cm, int part = 2);
120 
121  // Conversion from CM frame to LAB frame
122  double ConvertThetaCmToLab(double theta_cm, int part = 2);
123  double ConvertOmegaCmToLab(double theta_cm, int part = 2);
124  void ConvertCmToLab(double theta_cm, double omega_cm, double &theta_lab, double &omega_lab, int part = 2);
125 
126  // Graphs for conversions and kinematic/cross-section curves
127  // Frame_Lab -> TLab[ThetaLab] and Frame_Cm -> TLab[ThetaCm]
128  TGraph *KinVsTheta(double thmin = 0.0, double thmax = 180.0, int part = 2, bool Frame_Lab = true, bool Units_keV = true);
129  // Frame_Lab -> ThetaCm[ThetaLab] and Frame_Cm -> ThetaLab[ThetaCm]
130  TGraph *ThetaVsTheta(double thmin = 0.0, double thmax = 180.0, int part = 2, bool Frame_Lab = true);
131  // Frame_Lab -> dOmegaCm/dOmegaLab[ThetaLab] and Frame_Cm -> dOmegaLab/dOmegaCm[ThetaCm]
132  TGraph *OmegaVsTheta(double thmin = 0.0, double thmax = 180.0, int part = 2, bool Frame_Lab = true);
133  // Frame_Lab -> dSigma/dThetaLab[ThetaLab] and Frame_Cm -> dSigma/dThetaCm[ThetaCm]
134  TGraph *RutherfordVsTheta(double thmin = 1.0, double thmax = 179.0, int part = 2, bool Frame_Lab = true, bool Units_keV = true);
135 
136  void Print(Option_t *opt="");
137  void Clear(Option_t *opt="");
138 
139  void SetExcEnergy(double exc) { SetCmFrame(exc); }
140 
141 private:
142 
143  void SetCmFrame(double exc); // enables the reaction to be modified using excitation energy
144 
145  // USER INPUTS
147  double fTBeam;
148  bool fInverse;
149  double fExc;
150  double fM[4];
151 
152  // CM FRAME MOTION
153  double fQVal; // effective Q value (includes excitation)
154  double fS; // 'S' = M^2
155  double fInvariantMass;
156  double fCmTi;
157  double fCmTf;
158  double fCmE;
159  double fCmV;
160  double fCmP;
161  double fCmG;
162 
163  // PARTICLES IN CM FRAME
164  double fTCm[4];
165  double fECm[4];
166  double fPCm[4];
167  double fVCm[4];
168  double fGCm[4];
169 
170  // PARTICLE IN LAB FRAME
171  // Note that in the lab frame only the initial state (beam/targ) is fixed in the reaction
172  double fTLab[2];
173  double fELab[2];
174  double fPLab[2];
175  double fVLab[2];
176  double fGLab[2];
177  double fThetaMax[4]; // only nonzero for ejectile and recoil
178 
179  ClassDef(TReaction,1); // Calculates reaction parameters for scattering experiments
180 };
181 #endif