Merge pull request #219 from MarkKJones/primarykine

Add THcPrimaryKine and modified THcHallCSpectrometer

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 )) ) {

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
 

src/THcPrimaryKine.cxx

+/** \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()");
+}

src/THcPrimaryKine.h

+#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