From 6a969cf754a46a54a91fc066e3e434b50536c301 Mon Sep 17 00:00:00 2001
From: hallc-online <jlabhallconline@gmail.com>
Date: Fri, 21 Jul 2017 15:08:09 -0400
Subject: [PATCH] Add THcPrimaryKine and modified THcHallCSpectrometer
Add class THcPrimaryKine to calculate physics kinematics quantities for single arm
Modified THcHallCSpectrometer so that it would fill track energy and momentum 3 vector
---
src/THcHallCSpectrometer.cxx | 16 ++-
src/THcHallCSpectrometer.h | 1 +
src/THcPrimaryKine.cxx | 245 +++++++++++++++++++++++++++++++++++
src/THcPrimaryKine.h | 94 ++++++++++++++
4 files changed, 351 insertions(+), 5 deletions(-)
create mode 100644 src/THcPrimaryKine.cxx
create mode 100644 src/THcPrimaryKine.h
diff --git a/src/THcHallCSpectrometer.cxx b/src/THcHallCSpectrometer.cxx
index b912ceb..a4312d5 100644
--- a/src/THcHallCSpectrometer.cxx
+++ b/src/THcHallCSpectrometer.cxx
@@ -165,6 +165,7 @@ Int_t THcHallCSpectrometer::ReadDatabase( const TDatime& date )
{"pcentral", &fPcentral, kDouble },
{"theta_lab", &fTheta_lab, kDouble },
{"partmass", &fPartMass, kDouble },
+ {"phi_lab", &fPhi_lab, kDouble, 0, 1},
{"sel_using_scin", &fSelUsingScin, kInt, 0, 1},
{"sel_using_prune", &fSelUsingPrune, kInt, 0, 1},
{"sel_ndegreesmin", &fSelNDegreesMin, kDouble, 0, 1},
@@ -194,7 +195,7 @@ Int_t THcHallCSpectrometer::ReadDatabase( const TDatime& date )
// Default values
fSelUsingScin = 0;
fSelUsingPrune = 0;
-
+ fPhi_lab = 0.;
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
EnforcePruneLimits();
@@ -217,8 +218,8 @@ Int_t THcHallCSpectrometer::ReadDatabase( const TDatime& date )
fPcentral= fPcentral*(1.+fPCentralOffset/100.);
// Check that these offsets are in radians
fTheta_lab=fTheta_lab + fThetaCentralOffset*TMath::RadToDeg();
- Double_t ph = 0.0+fPhiOffset*TMath::RadToDeg();
-
+ Double_t ph = fPhi_lab+fPhiOffset*TMath::RadToDeg();
+ cout << " Central angles = " << fTheta_lab << endl;
SetCentralAngles(fTheta_lab, ph, false);
Double_t off_x = 0.0, off_y = 0.0, off_z = 0.0;
fPointingOffset.SetXYZ( off_x, off_y, off_z );
@@ -356,8 +357,13 @@ Int_t THcHallCSpectrometer::FindVertices( TClonesArray& tracks )
track->SetDp(sum[3]*100.0+fDeltaOffset); // Percent. (Don't think podd cares if it is % or fraction)
// There is an hpcentral_offset that needs to be applied somewhere.
// (happly_offs)
- track->SetMomentum(fPcentral*(1+track->GetDp()/100.0));
-
+ Double_t ptemp = fPcentral*(1+track->GetDp()/100.0);
+ Double_t mtemp = fPartMass/1000.; // fPartMass in MeV , convert to GeV
+ track->SetMomentum(ptemp);
+ track->SetEnergy(ptemp*ptemp+mtemp*mtemp);
+ TVector3 pvect_temp;
+ TransportToLab(track->GetP(),track->GetTTheta(),track->GetTPhi(),pvect_temp);
+ track->SetPvect(pvect_temp);
}
if (fHodo==0 || (( fSelUsingScin == 0 ) && ( fSelUsingPrune == 0 )) ) {
diff --git a/src/THcHallCSpectrometer.h b/src/THcHallCSpectrometer.h
index 7dfbe25..d0302a6 100644
--- a/src/THcHallCSpectrometer.h
+++ b/src/THcHallCSpectrometer.h
@@ -132,6 +132,7 @@ protected:
Double_t fOopCentralOffset; //Offset of central out-of-plane angle (rad)
Double_t fPCentralOffset; // Offset Central spectrometer momentum (%)
Double_t fTheta_lab; // Central spectrometer angle (deg)
+ Double_t fPhi_lab; // Central spectrometer angle (deg)
// For spectrometer central momentum use fPcentral in THaSpectrometer.h
// THaScintillator *sc_ref; // calculate time track hits this plane
diff --git a/src/THcPrimaryKine.cxx b/src/THcPrimaryKine.cxx
new file mode 100644
index 0000000..c4dd070
--- /dev/null
+++ b/src/THcPrimaryKine.cxx
@@ -0,0 +1,245 @@
+/** \class THcPrimaryKine
+ \ingroup Physics
+
+\brief Class for the Calculate kinematics of scattering of the primary (beam) particle.
+These are usually the electron kinematics.
+*/
+
+#include "THcPrimaryKine.h"
+#include "THaTrackingModule.h"
+#include "THcGlobals.h"
+#include "THcParmList.h"
+#include "THaRunBase.h"
+#include "THaRunParameters.h"
+#include "THaBeam.h"
+#include "VarDef.h"
+#include "TMath.h"
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+
+using namespace std;
+
+ClassImp(THcPrimaryKine)
+
+//_____________________________________________________________________________
+THcPrimaryKine::THcPrimaryKine( const char* name, const char* description,
+ const char* spectro, Double_t particle_mass,
+ Double_t target_mass ) :
+ THaPhysicsModule(name,description), fM(particle_mass),
+ fMA(target_mass), fSpectroName(spectro), fSpectro(NULL), fBeam(NULL)
+{
+ // Standard constructor. Must specify particle mass. Incident particles
+ // are assumed to be along z_lab.
+
+}
+
+//_____________________________________________________________________________
+THcPrimaryKine::THcPrimaryKine( const char* name, const char* description,
+ const char* spectro, const char* beam,
+ Double_t target_mass )
+ : THaPhysicsModule(name,description), fM(-1.0), fMA(target_mass),
+ fSpectroName(spectro), fBeamName(beam), fSpectro(NULL), fBeam(NULL)
+{
+ // Constructor with specification of optional beam module.
+ // Particle mass will normally come from the beam module.
+
+}
+//_____________________________________________________________________________
+
+//_____________________________________________________________________________
+THcPrimaryKine::~THcPrimaryKine()
+{
+ // Destructor
+
+ DefineVariables( kDelete );
+}
+
+//_____________________________________________________________________________
+void THcPrimaryKine::Clear( Option_t* opt )
+{
+ // Clear event-by-event variables.
+
+ THaPhysicsModule::Clear(opt);
+ fQ2 = fOmega = fW = fW2 = fXbj = fScatAngle = fScatAngle_deg = fEpsilon = fQ3mag
+ = fThetaQ = fPhiQ = kBig;
+ // Clear the 4-vectors
+ fP0.SetXYZT(kBig,kBig,kBig,kBig);
+ fP1 = fA = fA1 = fQ = fP0;
+
+}
+
+//_____________________________________________________________________________
+Int_t THcPrimaryKine::DefineVariables( EMode mode )
+{
+ // Define/delete global variables.
+
+ if( mode == kDefine && fIsSetup ) return kOK;
+ fIsSetup = ( mode == kDefine );
+
+ RVarDef vars[] = {
+ { "Q2", "4-momentum transfer squared (GeV^2)", "fQ2" },
+ { "omega", "Energy transfer (GeV)", "fOmega" },
+ { "W2", "Invariant mass (GeV^2) for Mp ", "fW2" },
+ { "W", "Sqrt(Invariant mass) for Mp ", "fW" },
+ { "x_bj", "Bjorken x", "fXbj" },
+ { "scat_ang_rad", "Scattering angle (rad)", "fScatAngle" },
+ { "scat_ang_deg", "Scattering angle (deg)", "fScatAngle_deg" },
+ { "epsilon", "Virtual photon polarization factor", "fEpsilon" },
+ { "q3m", "Magnitude of 3-momentum transfer", "fQ3mag" },
+ { "th_q", "Theta of 3-momentum vector (rad)", "fThetaQ" },
+ { "ph_q", "Phi of 3-momentum vector (rad)", "fPhiQ" },
+ { "nu", "Energy transfer (GeV)", "fOmega" },
+ { "q_x", "x-cmp of Photon vector in the lab", "fQ.X()" },
+ { "q_y", "y-cmp of Photon vector in the lab", "fQ.Y()" },
+ { "q_z", "z-cmp of Photon vector in the lab", "fQ.Z()" },
+ { 0 }
+ };
+ return DefineVarsFromList( vars, mode );
+}
+
+//_____________________________________________________________________________
+THaAnalysisObject::EStatus THcPrimaryKine::Init( const TDatime& run_time )
+{
+ // Initialize the module.
+ // Locate the spectrometer apparatus named in fSpectroName and save
+ // pointer to it.
+
+ fSpectro = dynamic_cast<THaTrackingModule*>
+ ( FindModule( fSpectroName.Data(), "THaTrackingModule"));
+ if( !fSpectro )
+ return fStatus;
+
+ // Optional beam apparatus
+ if( fBeamName.Length() > 0 ) {
+ fBeam = dynamic_cast<THaBeamModule*>
+ ( FindModule( fBeamName.Data(), "THaBeamModule") );
+ if( !fBeam )
+ return fStatus;
+ if( fM <= 0.0 )
+ fM = fBeam->GetBeamInfo()->GetM();
+ }
+
+ // Standard initialization. Calls this object's DefineVariables().
+ if( THaPhysicsModule::Init( run_time ) != kOK )
+ return fStatus;
+
+ if( fM <= 0.0 ) {
+ Error( Here("Init"), "Particle mass not defined. Module "
+ "initialization failed" );
+ fStatus = kInitError;
+ }
+ return fStatus;
+}
+
+//_____________________________________________________________________________
+Int_t THcPrimaryKine::Process( const THaEvData& )
+{
+ // Calculate electron kinematics for the Golden Track of the spectrometer
+ if( !IsOK() || !gHaRun ) return -1;
+
+ THaTrackInfo* trkifo = fSpectro->GetTrackInfo();
+ if( !trkifo || !trkifo->IsOK() ) return 1;
+ // Determine 4-momentum of incident particle.
+ // If a beam module given, use it to get the beam momentum. This
+ // module may apply corrections for beam energy loss, variations, etc.
+ if( fBeam ) {
+ fP0.SetVectM( fBeam->GetBeamInfo()->GetPvect(), fM );
+ } else {
+ // If no beam given, assume beam along z_lab
+ Double_t p_in = gHaRun->GetParameters()->GetBeamP();
+ fP0.SetXYZM( 0.0, 0.0, p_in, fM );
+ }
+
+ fP1.SetVectM( trkifo->GetPvect(), fM );
+ fA.SetXYZM( 0.0, 0.0, 0.0, fMA ); // Assume target at rest
+
+ // proton mass (for x_bj)
+ const Double_t Mp = 0.93827;
+ fMp.SetXYZM( 0.0, 0.0, 0.0, Mp ); // Assume target at rest
+
+
+ // Standard electron kinematics
+ fQ = fP0 - fP1;
+ fQ2 = -fQ.M2();
+ fQ3mag = fQ.P();
+ fOmega = fQ.E();
+ fA1 = fA + fQ;
+ // fW2 = fA1.M2();
+ fMp1 = fMp + fQ;
+ fW2 = fMp1.M2();
+ if (fW2>0) fW = TMath::Sqrt(fW2);
+ fScatAngle = fP0.Angle( fP1.Vect() );
+ fEpsilon = 1.0 / ( 1.0 + 2.0*fQ3mag*fQ3mag/fQ2*
+ TMath::Power( TMath::Tan(fScatAngle/2.0), 2.0 ));
+ fScatAngle_deg = fScatAngle*TMath::RadToDeg();
+ fThetaQ = fQ.Theta();
+ fPhiQ = fQ.Phi();
+ fXbj = fQ2/(2.0*Mp*fOmega);
+
+ fDataValid = true;
+ return 0;
+}
+
+Int_t THcPrimaryKine::ReadDatabase( const TDatime& date )
+{
+
+ static const char* const here = "THcPrimaryKine::ReadDatabase";
+
+#ifdef WITH_DEBUG
+ cout << "In THcPrimaryKine::ReadDatabase()" << endl;
+#endif
+ //
+ char prefix[2];
+
+ cout << " GetName() " << GetName() << endl;
+
+ prefix[0]=tolower(GetName()[0]);
+ prefix[1]='\0';
+
+ DBRequest list[]={
+ {"targmass_amu", &fMA_amu, kDouble, 0, 1},
+ {0}
+ };
+ gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+ //
+ fMA= fMA_amu*931.5/1000.;
+ return kOK;
+}
+
+//_____________________________________________________________________________
+void THcPrimaryKine::SetMass( Double_t m )
+{
+ if( !IsInit())
+ fM = m;
+ else
+ PrintInitError("SetMass()");
+}
+
+//_____________________________________________________________________________
+void THcPrimaryKine::SetTargetMass( Double_t m )
+{
+ if( !IsInit())
+ fMA = m;
+ else
+ PrintInitError("SetTargetMass()");
+}
+
+//_____________________________________________________________________________
+void THcPrimaryKine::SetSpectrometer( const char* name )
+{
+ if( !IsInit())
+ fSpectroName = name;
+ else
+ PrintInitError("SetSpectrometer()");
+}
+
+//_____________________________________________________________________________
+void THcPrimaryKine::SetBeam( const char* name )
+{
+ if( !IsInit())
+ fBeamName = name;
+ else
+ PrintInitError("SetBeam()");
+}
diff --git a/src/THcPrimaryKine.h b/src/THcPrimaryKine.h
new file mode 100644
index 0000000..b0b0ad2
--- /dev/null
+++ b/src/THcPrimaryKine.h
@@ -0,0 +1,94 @@
+#ifndef ROOT_THcPrimaryKine
+#define ROOT_THcPrimaryKine
+
+//////////////////////////////////////////////////////////////////////////
+//
+// THaPrimaryKine
+//
+//////////////////////////////////////////////////////////////////////////
+
+#include "THaPhysicsModule.h"
+#include "TLorentzVector.h"
+#include "TString.h"
+
+class THaTrackingModule;
+class THaBeamModule;
+typedef TLorentzVector FourVect;
+
+class THcPrimaryKine : public THaPhysicsModule {
+
+public:
+ THcPrimaryKine( const char* name, const char* description,
+ const char* spectro = "",
+ Double_t particle_mass = 0.0, /* GeV */
+ Double_t target_mass = 0.0 /* GeV */ );
+ THcPrimaryKine( const char* name, const char* description,
+ const char* spectro, const char* beam,
+ Double_t target_mass = 0.0 /* GeV */ );
+ virtual ~THcPrimaryKine();
+
+ virtual void Clear( Option_t* opt="" );
+
+ Double_t GetQ2() const { return fQ2; }
+ Double_t GetOmega() const { return fOmega; }
+ Double_t GetNu() const { return fOmega; }
+ Double_t GetW2() const { return fW2; }
+ Double_t GetXbj() const { return fXbj; }
+ Double_t GetScatAngle() const { return fScatAngle; }
+ Double_t GetEpsilon() const { return fEpsilon; }
+ Double_t GetQ3mag() const { return fQ3mag; }
+ Double_t GetThetaQ() const { return fThetaQ; }
+ Double_t GetPhiQ() const { return fPhiQ; }
+ Double_t GetMass() const { return fM; }
+ Double_t GetTargetMass() const { return fMA; }
+
+ const FourVect* GetP0() const { return &fP0; }
+ const FourVect* GetP1() const { return &fP1; }
+ const FourVect* GetA() const { return &fA; }
+ const FourVect* GetA1() const { return &fA1; }
+ const FourVect* GetQ() const { return &fQ; }
+
+ virtual EStatus Init( const TDatime& run_time );
+ virtual Int_t Process( const THaEvData& );
+ virtual Int_t ReadDatabase(const TDatime& date);
+ void SetMass( Double_t m );
+ void SetTargetMass( Double_t m );
+ void SetSpectrometer( const char* name );
+ void SetBeam( const char* name );
+
+protected:
+
+ Double_t fQ2; // 4-momentum transfer squared (GeV^2)
+ Double_t fOmega; // Energy transfer (GeV)
+ Double_t fW2; // s = Invariant mass using Mp (GeV^2)
+ Double_t fW; // sqrt(s) (GeV)
+ Double_t fXbj; // x Bjorken
+ Double_t fScatAngle; // Scattering angle (rad)
+ Double_t fScatAngle_deg; // Scattering angle (deg)
+ Double_t fEpsilon; // Virtual photon polarization factor
+ Double_t fQ3mag; // Magnitude of 3-momentum transfer
+ Double_t fThetaQ; // Theta of 3-momentum vector (rad)
+ Double_t fPhiQ; // Phi of 3-momentum transfer (rad)
+ FourVect fP0; // Beam 4-momentum
+ FourVect fP1; // Scattered electron 4-momentum
+ FourVect fA; // Target 4-momentum
+ FourVect fA1; // Recoil system 4-momentum
+ FourVect fQ; // Momentum transfer 4-vector
+ FourVect fMp; // Mp 4-momentum
+ FourVect fMp1; // Recoil Mp 4-momentum
+
+ Double_t fM; // Mass of incident particle (GeV/c^2)
+ Double_t fMA; // Target mass (GeV/c^2)
+ Double_t fMA_amu; // Target mass (amu)
+
+ virtual Int_t DefineVariables( EMode mode = kDefine );
+
+ TString fSpectroName; // Name of spectrometer to consider
+ TString fBeamName; // Name of beam position apparatus
+ THaTrackingModule* fSpectro; // Pointer to spectrometer object
+ THaBeamModule* fBeam; // Pointer to beam position apparatus
+
+ ClassDef(THcPrimaryKine,0) //Single arm kinematics module
+};
+
+#endif
--
GitLab