Hall C ROOT/C++ Analyzer (hcana)
Loading...
Searching...
No Matches
THcSecondaryKine.cxx
Go to the documentation of this file.
1
7#include "THcSecondaryKine.h"
8#include "THcPrimaryKine.h"
10#include "THcGlobals.h"
11#include "THcParmList.h"
12#include "VarDef.h"
13#include "TLorentzVector.h"
14#include "TVector3.h"
15#include "TMath.h"
16#include "TRotation.h"
17
18#include <cstring>
19#include <cstdio>
20#include <cstdlib>
21#include <iostream>
22using namespace std;
23
25
26//_____________________________________________________________________________
27THcSecondaryKine::THcSecondaryKine( const char* name, const char* description,
28 const char* secondary_spectro,
29 const char* primary_kine,
30 Double_t secondary_mass ) :
31 THaPhysicsModule(name,description), fMX(secondary_mass),
32 fSpectroName(secondary_spectro), fSpectro(NULL),
33 fPrimaryName(primary_kine), fPrimary(NULL)
34{
35 // Constructor
36}
37
38//_____________________________________________________________________________
44
45//_____________________________________________________________________________
59
60//_____________________________________________________________________________
62{
63 // Define/delete global variables.
64
65 if( mode == kDefine && fIsSetup ) return kOK;
66 fIsSetup = ( mode == kDefine );
67
68 RVarDef vars[] = {
69 { "th_xq", "Polar angle of detected particle with q (rad)",
70 "fTheta_xq" },
71 { "ph_xq", "Azimuth of detected particle with scattering plane (rad)",
72 "fPhi_xq" },
73 { "th_bq", "Polar angle of recoil system with q (rad)", "fTheta_bq" },
74 { "ph_bq", "Azimuth of recoil system with scattering plane (rad)",
75 "fPhi_bq" },
76 { "xangle", "Angle of detected particle with scattered electron (rad)",
77 "fXangle" },
78 { "pmiss", "Missing momentum magnitude (GeV), nuclear physics "
79 "definition (-pB)", "fPmiss" },
80 { "pmiss_x", "x-component of p_miss wrt q (GeV)", "fPmiss_x" },
81 { "pmiss_y", "y-component of p_miss wrt q (GeV)", "fPmiss_y" },
82 { "pmiss_z", "z-component of p_miss, along q (GeV)", "fPmiss_z" },
83 { "emiss_nuc", "Missing energy (GeV), nuclear physics definition "
84 "omega-Tx-Tb", "fEmiss_nuc" },
85 { "emiss", "Missing energy (GeV), ENGINE definition "
86 "omega-Mt-Ex", "fEmiss" },
87 { "Mrecoil", "Invariant mass of recoil system (GeV)", "fMrecoil" },
88 { "MMpi", "Missing mass under assumption hadron is a pion (GeV)", "fMMpi" },
89 { "MMK", "Missing mass under assumption hadron is a kaon (GeV)", "fMMK" },
90 { "MMp", "Missing mass under assumption hadron is a proton (GeV)", "fMMp" },
91 { "Erecoil", "Total energy of recoil system (GeV)", "fErecoil" },
92 { "Prec_x", "x-component of recoil system in lab (GeV/c)", "fB.X()" },
93 { "Prec_y", "y-component of recoil system in lab (GeV/c)", "fB.Y()" },
94 { "Prec_z", "z-component of recoil system in lab (GeV/c)", "fB.Z()" },
95 { "tx", "Kinetic energy of detected particle (GeV)", "fTX" },
96 { "tb", "Kinetic energy of recoil system (GeV)", "fTB" },
97 { "px_cm", "Magnitude of X momentum in CM system (GeV)", "fPX_cm" },
98 { "thx_cm", "Polar angle of X in CM system wrt q (rad)", "fTheta_x_cm" },
99 { "phx_cm", "Azimuth of X in CM system wrt q (rad)", "fPhi_x_cm" },
100 { "thb_cm", "Polar angle of recoil systm in CM wrt q (rad)",
101 "fTheta_b_cm" },
102 { "phb_cm", "Azimuth of recoil system in CM wrt q (rad)", "fPhi_b_cm" },
103 { "tx_cm", "Kinetic energy of X in CM (GeV)", "fTX_cm" },
104 { "tb_cm", "Kinetic energy of B in CM (GeV)", "fTB_cm" },
105 { "t_tot_cm", "Total CM kinetic energy", "fTtot_cm" },
106 { "MandelS", "Mandelstam s for secondary vertex (GeV^2)", "fMandelS" },
107 { "MandelT", "Mandelstam t for secondary vertex (GeV^2)", "fMandelT" },
108 { "MandelU", "Mandelstam u for secondary vertex (GeV^2)", "fMandelU" },
109 { 0 }
110 };
111 return DefineVarsFromList( vars, mode );
112}
113
114//_____________________________________________________________________________
116{
117 // Initialize the module.
118 // Locate the kinematics module named in fPrimaryName, which describes
119 // the primary interaction kinematics, and save pointer to it.
120
121 // Standard initialization. Calls this object's DefineVariables().
122
123 // cout << "*************&&&&&&&&&&&&&&************" << endl;
124 // cout << "In THcSecondaryKine::Init() " << endl;
125
126 fStatus = kOK;
127
128 fSpectro = dynamic_cast<THcHallCSpectrometer*>
129 ( FindModule( fSpectroName.Data(), "THcHallCSpectrometer"));
130 if( !fSpectro ) {
132 return fStatus;
133 }
134
135 fPrimary = dynamic_cast<THcPrimaryKine*>
136 ( FindModule( fPrimaryName.Data(), "THcPrimaryKine"));
137 if(!fPrimary) {
139 return fStatus;
140 }
141
142 if( (fStatus=THaPhysicsModule::Init( run_time )) != kOK ) {
143 return fStatus;
144 }
145
146 return fStatus;
147}
148
149//_____________________________________________________________________________
151{
152 // Calculate the kinematics.
153 if( !IsOK() ) return -1;
154
155 //Get secondary particle mass
156 //fMX = dynamic_cast <THcHallCSpectrometer*> (fSpectro)->GetParticleMass();
157
158 // Tracking information from the secondary spectrometer
160 if( !trkifo || !trkifo->IsOK() ) return 1;
161
162 // Require valid input data
163 if( !fPrimary->DataValid() ) return 2;
164
165 // Measured momentum of secondary particle X in lab
166 Double_t xptar = trkifo->GetTheta() + fOopCentralOffset;
167 TVector3 pvect;
168 fSpectro->TransportToLab(trkifo->GetP(), xptar, trkifo->GetPhi(), pvect);
169
170 // 4-momentum of X
171 fX.SetVectM( pvect, fMX );
172
173 // 4-momenta of the the primary interaction
174 const TLorentzVector* pA = fPrimary->GetA(); // Initial target
175 const TLorentzVector* pA1 = fPrimary->GetA1(); // Final target
176 const TLorentzVector* pQ = fPrimary->GetQ(); // Momentum xfer
177 const TLorentzVector* pP1 = fPrimary->GetP1(); // Final electron
178 Double_t omega = fPrimary->GetOmega(); // Energy xfer
179
180 // 4-momentum of undetected recoil system B
181 fB = *pA1 - fX;
182
183 // Angle of X with scattered primary particle
184 fXangle = fX.Angle( pP1->Vect());
185
186 // Angles of X and B wrt q-vector
187 // xq and bq are the 3-momentum vectors of X and B expressed in
188 // the coordinate system where q is the z-axis and the x-axis
189 // lies in the scattering plane (defined by q and e') and points
190 // in the direction of e', so the out-of-plane angle lies within
191 // -90<phi_xq<90deg if X is detected on the downstream/forward side of q.
192 TRotation rot_to_q;
193 rot_to_q.SetZAxis( pQ->Vect(), pP1->Vect()).Invert();
194 TVector3 xq = fX.Vect();
195 TVector3 bq = fB.Vect();
196 xq *= rot_to_q;
197 bq *= rot_to_q;
198 fTheta_xq = xq.Theta(); //"theta_xq"
199 fPhi_xq = xq.Phi(); //"out-of-plane angle", "phi"
200 fTheta_bq = bq.Theta();
201 fPhi_bq = bq.Phi();
202
203 // Missing momentum and components wrt q-vector
204 // The definition of p_miss as the negative of the undetected recoil
205 // momentum is the standard nuclear physics convention.
206 TVector3 p_miss = -bq;
207 fPmiss = p_miss.Mag(); //=fB.P()
208 //The missing momentum components in the q coordinate system.
209 //FIXME: just store p_miss in member variable
210 fPmiss_x = p_miss.X();
211 fPmiss_y = p_miss.Y();
212 fPmiss_z = p_miss.Z();
213
214 // Invariant mass of the recoil system, a.k.a. "missing mass".
215 // This invariant mass equals MB(ground state) plus any excitation energy.
216 fMrecoil = fB.M();
217
218 // Kinetic energies of X and B
219 fTX = fX.E() - fMX;
220 fTB = fB.E() - fMrecoil;
221
222 // Standard nuclear physics definition of "missing energy":
223 // binding energy of X in the target (= removal energy of X).
224 // NB: If X is knocked out of a lower shell, the recoil system carries
225 // a significant excitation energy. This excitation is included in Emiss
226 // here, as it should, since it results from the binding of X.
227 fEmiss_nuc = omega - fTX - fTB;
228 // Target Mass + Beam Energy - scatt ele energy - hadron total energy
229 fEmiss = omega + pA->M() - fX.E();
230
231 // In production reactions, the "missing energy" is defined
232 // as the total energy of the undetected recoil system.
233 // This is the "missing mass", Mrecoil, plus any kinetic energy.
234 fErecoil = fB.E();
235
236 // SJDK - 15/04/21 - Added new missing mass calculations
237 // Determination of the missing mass under various assumptions, utilises database read in to establish which
238 // corrections to use to determine values
239 // It might be easier/nicer to just manipulate the missing energy in the 4-vector fB to determine these quantities
240 // Could also tie the calculation into the actual read in value a bit more too
241 // Hadron in standard.kinematics is a charged pion
242 if ( fMX > 0.12957018 && fMX < 0.14957018 ){
244 fMMK = sqrt(abs((pow(fEmiss+(sqrt((fMass_pi*fMass_pi)+(pow((pvect.Mag()), 2))))-(sqrt((fMass_K*fMass_K)+(pow((pvect.Mag()), 2)))), 2)-(fPmiss*fPmiss))));
245 fMMp = sqrt(abs((pow(fEmiss+(sqrt((fMass_pi*fMass_pi)+(pow((pvect.Mag()), 2))))-(sqrt((fMass_p*fMass_p)+(pow((pvect.Mag()), 2)))), 2)-(fPmiss*fPmiss))));
246 }
247 // Hadron in standard.kinematics is a charged kaon
248 else if (fMX > 0.483677 && fMX < 0.503677 ){
249 fMMpi = sqrt(abs((pow(fEmiss+(sqrt((fMass_K*fMass_K)+(pow((pvect.Mag()), 2))))-(sqrt((fMass_pi*fMass_pi)+(pow((pvect.Mag()), 2)))), 2)-(fPmiss*fPmiss))));
251 fMMp = sqrt(abs((pow(fEmiss+(sqrt((fMass_K*fMass_K)+(pow((pvect.Mag()), 2))))-(sqrt((fMass_p*fMass_p)+(pow((pvect.Mag()), 2)))), 2)-(fPmiss*fPmiss))));
252 }
253 // Hadron in standard.kinematics is a proton
254 else if (fMX > 0.92828 && fMX < 0.94828 ){
255 fMMpi = sqrt(abs((pow(fEmiss+(sqrt((fMass_p*fMass_p)+(pow((pvect.Mag()), 2))))-(sqrt((fMass_pi*fMass_pi)+(pow((pvect.Mag()), 2)))), 2)-(fPmiss*fPmiss))));
256 fMMK = sqrt(abs((pow(fEmiss+(sqrt((fMass_p*fMass_p)+(pow((pvect.Mag()), 2))))-(sqrt((fMass_K*fMass_K)+(pow((pvect.Mag()), 2)))), 2)-(fPmiss*fPmiss))));
258 }
259 // Any other condition, try and determine something from whatever the hell the specified mass was
260 else if (fMX != 0 ){
261 fMMpi = sqrt(abs((pow(fEmiss+(sqrt((fMX*fMX)+(pow((pvect.Mag()), 2))))-(sqrt((fMass_pi*fMass_pi)+(pow((pvect.Mag()), 2)))), 2)-(fPmiss*fPmiss))));
262 fMMK = sqrt(abs((pow(fEmiss+(sqrt((fMX*fMX)+(pow((pvect.Mag()), 2))))-(sqrt((fMass_K*fMass_K)+(pow((pvect.Mag()), 2)))), 2)-(fPmiss*fPmiss))));
263 fMMp = sqrt(abs((pow(fEmiss+(sqrt((fMX*fMX)+(pow((pvect.Mag()), 2))))-(sqrt((fMass_p*fMass_p)+(pow((pvect.Mag()), 2)))), 2)-(fPmiss*fPmiss))));
264 }
265 // If the secondary mass was 0 or negative for some reason, this is just the "anything else" junk variable output basically
266 // Note, add an extra condition for the mass being 0 if it is relevant for you for some reason!
267 else{
268 fMMpi = -1000;
269 fMMK = -1000;
270 fMMp = -1000;
271 }
272
273 // Calculate some interesting quantities in the CM system of A'.
274 // NB: If the target is initially at rest, the A'-vector (spatial part)
275 // is the same as the q-vector, so we could just reuse the q rotation
276 // matrix. The following is completely general, i.e. allows for a moving
277 // target.
278
279 // Boost of the A' system.
280 // NB: ROOT's Boost() boosts particle to lab if the boost vector points
281 // along the particle's lab velocity, so we need the negative here to
282 // boost from the lab to the particle frame.
283 Double_t beta = pA1->P()/(pA1->E());
284 TVector3 boost(0.,0.,-beta);
285
286 // CM vectors of X and B.
287 // Express X and B in the frame where q is along the z-axis
288 // - the typical head-on collision picture.
289 TRotation rot_to_A1;
290 rot_to_A1.SetZAxis( pA1->Vect(), pP1->Vect()).Invert();
291 TVector3 x_cm_vect = fX.Vect();
292 x_cm_vect *= rot_to_A1;
293 TLorentzVector x_cm(x_cm_vect,fX.E());
294 x_cm.Boost(boost);
295 TLorentzVector b_cm(-x_cm.Vect(), pA1->E()-x_cm.E());
296 fPX_cm = x_cm.P();
297 // pB_cm, by construction, is the same as pX_cm.
298
299 // CM angles of X and B in the A' frame
300 fTheta_x_cm = x_cm.Theta();
301 fPhi_x_cm = x_cm.Phi();
302 fTheta_b_cm = b_cm.Theta();
303 fPhi_b_cm = b_cm.Phi();
304
305 // CM kinetic energies of X and B and total kinetic energy
306 fTX_cm = x_cm.E() - fMX;
307 fTB_cm = b_cm.E() - fMrecoil;
309
310 // Mandelstam variables for the secondary vertex.
311 // These variables are defined for the reaction gA->XB,
312 // where g is the virtual photon (with momentum q), and A, X, and B
313 // are as before.
314
315 fMandelS = (*pQ+*pA).M2();
316 fMandelT = (*pQ-fX).M2();
317 fMandelU = (*pQ-fB).M2();
318
319 fDataValid = true;
320 return 0;
321}
322
323//_____________________________________________________________________________
325{
326 // cout << "In THcSecondaryKine::ReadDatabase() " << endl;
327 fMass_pi = 0.13957018;
328 fMass_K = 0.493677;
329 fMass_p = 0.93828;
330 char prefix[2];
331
332 prefix[0] = tolower(GetName()[0]);
333 prefix[1] = '\0';
334
335 fOopCentralOffset = 0.0;
336 DBRequest list[]={
337 {"_oopcentral_offset",&fOopCentralOffset,kDouble, 0, 1},
338 {"partmass", &fMX, kDouble },
339 {0}
340 };
341 gHcParms->LoadParmValues((DBRequest*)&list,prefix);
342 cout << "THcSecondaryKine particleMASS: " << fMX << endl;
343
344 return kOK;
345}
346
347//_____________________________________________________________________________
349{
350 if( !IsInit())
351 fMX = m;
352 else
353 PrintInitError("SetMX()");
354}
355
356//_____________________________________________________________________________
357void THcSecondaryKine::SetSpectrometer( const char* name )
358{
359 if( !IsInit())
360 fSpectroName = name;
361 else
362 PrintInitError("SetSpectrometer()");
363}
364
365//_____________________________________________________________________________
366void THcSecondaryKine::SetPrimary( const char* name )
367{
368 if( !IsInit())
369 fPrimaryName = name;
370 else
371 PrintInitError("SetPrimary()");
372}
int Int_t
const Data_t kBig
double Double_t
const char Option_t
Option_t Option_t TPoint TPoint const char mode
ClassImp(VDC::AnalyticTTDConv) using namespace std
R__EXTERN class THcParmList * gHcParms
Definition THcGlobals.h:11
static Int_t DefineVarsFromList(const void *list, EType type, EMode mode, const char *def_prefix, const TObject *obj, const char *prefix, const char *here, const char *comment_subst="")
THaAnalysisObject * FindModule(const char *name, const char *classname, bool do_error=true)
Bool_t IsOK() const
Bool_t IsInit() const
bool DataValid() const
void PrintInitError(const char *here)
virtual void Clear(Option_t *opt="")
virtual void TransportToLab(Double_t p, Double_t th, Double_t ph, TVector3 &pvect) const
Double_t GetTheta() const
Double_t GetP() const
Double_t GetPhi() const
Bool_t IsOK() const
THaTrackInfo * GetTrackInfo()
A standard Hall C spectrometer apparatus.
Int_t LoadParmValues(const DBRequest *list, const char *prefix="")
Retrieve parameter values from the parameter cache.
Class for the Calculate kinematics of scattering of the primary (beam) particle. These are usually th...
Double_t GetOmega() const
const FourVect * GetA1() const
const FourVect * GetQ() const
const FourVect * GetA() const
const FourVect * GetP1() const
Class for the Calculate kinematics of scattering of the secondary (hadron) particle.
THcPrimaryKine * fPrimary
virtual void Clear(Option_t *opt="")
virtual Int_t DefineVariables(EMode mode=kDefine)
TLorentzVector fB
TLorentzVector fX
virtual Int_t Process(const THaEvData &)
THcSecondaryKine(const char *name, const char *description, const char *secondary_spectro="", const char *primary_kine="", Double_t secondary_mass=0.0)
THcHallCSpectrometer * fSpectro
virtual Int_t ReadDatabase(const TDatime &date)
void SetPrimary(const char *name)
void SetMX(Double_t m)
void SetSpectrometer(const char *name)
Double_t Angle(const TVector3 &v) const
void SetXYZT(Double_t x, Double_t y, Double_t z, Double_t t)
TVector3 Vect() const
Double_t M() const
Double_t P() const
void Boost(const TVector3 &)
Double_t Theta() const
Double_t Phi() const
Double_t E() const
void SetVectM(const TVector3 &spatial, Double_t mass)
const char * GetName() const override
TRotation & SetZAxis(const TVector3 &axis)
TRotation & Invert()
const char * Data() const
Double_t Z() const
Double_t Phi() const
Double_t Y() const
Double_t X() const
Double_t Mag() const
Double_t Theta() const
double beta(double x, double y)
RVec< PromoteType< T > > abs(const RVec< T > &v)
RVec< PromoteTypes< T0, T1 > > pow(const RVec< T0 > &v, const T1 &y)
LVector boost(const LVector &v, const BoostVector &b)
Expr< UnaryOp< Sqrt< T >, Expr< A, T, D, D2, R >, T >, T, D, D2, R > sqrt(const Expr< A, T, D, D2, R > &rhs)
STL namespace.
TMarker m