Coincidence Time physics module

1) Created Coincidence Time Physics Module
2) Created method to get coin. trig times in the THcTrigDet.cxx class

src/THcCoinTime.cxx

+/** \class THcCoinTime
+    \ingroup PhysMods
+
+\brief Class for calculating and adding the Coincidence Time in the Tree.
+
+//Author: Carlos Yero
+//Date: April 27, 2018
+*/
+#include "THaEvData.h"
+#include "THaCutList.h"
+#include "VarDef.h"
+#include "VarType.h"
+#include "TClonesArray.h"
+
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+
+#include "THcCoinTime.h"
+#include "THcTrigDet.h"
+#include "THaApparatus.h"
+#include "THcHodoHit.h"
+#include "THcGlobals.h"
+#include "THcParmList.h"
+
+using namespace std;
+
+//_____________________________________________________________________________
+THcCoinTime::THcCoinTime (const char *name, const char* description, const char* hadArmName, 
+			  const char* elecArmName, const char* coinname) :
+  
+  THaPhysicsModule(name, description), 
+  fhadArmName(hadArmName),                 //initialize spectro names
+  felecArmName(elecArmName),
+  fCoinDetName(coinname), 
+  fhadSpectro(NULL),                      //initialize spectro objects
+  felecSpectro(NULL),
+  fCoinDet(NULL)
+{
+
+}
+
+//_____________________________________________________________________________
+THcCoinTime::~THcCoinTime()
+{
+  //Destructor
+
+  RemoveVariables();
+
+}
+//_____________________________________________________________________________
+void THcCoinTime::Clear( Option_t* opt )
+{
+  // Clear all event-by-event variables.
+  
+  //  THaPhysicsModule::Clear(opt);
+  fROC1_epCoinTime=kBig;
+  fROC2_epCoinTime=kBig;
+  fROC1_eKCoinTime=kBig;
+  fROC2_eKCoinTime=kBig;
+  fROC1_ePiCoinTime=kBig;
+  fROC2_ePiCoinTime=kBig;
+  fROC1_ePosCoinTime=kBig;
+  fROC2_ePosCoinTime=kBig;
+  fROC1_RAW_CoinTime=kBig;
+  fROC2_RAW_CoinTime=kBig;
+}
+
+//_____________________________________________________________________________
+void THcCoinTime::Reset( Option_t* opt)
+// Clear event-by-event data
+{
+  Clear(opt);
+}
+
+
+//_____________________________________________________________________________
+THaAnalysisObject::EStatus THcCoinTime::Init( const TDatime& run_time )
+{
+  // Initialize THcCoinTime physics module
+  
+  cout << "*************************************************" << endl;
+  cout << "Initializing THcCointTime Physics Modue" << endl;
+  cout << "Hadron Arm   -------> " << fhadArmName << endl;
+  cout << "Electron Arm -------> " << felecArmName << endl;
+  cout << "TrigDet  -------> " << fCoinDetName << endl;
+  cout << "**************************************************" << endl;
+
+  fStatus = kOK;
+
+  fhadSpectro = dynamic_cast<THcHallCSpectrometer*>
+    ( FindModule( fhadArmName.Data(), "THcHallCSpectrometer"));
+  if( !fhadSpectro ) {
+    cout << "THcCoinTime module  Cannnot find Hadron Arm = " <<  fhadArmName.Data() << endl;
+    fStatus = kInitError;
+    return fStatus;
+  }
+  
+  felecSpectro = dynamic_cast<THcHallCSpectrometer*>
+    ( FindModule( felecArmName.Data(), "THcHallCSpectrometer"));
+  if( !felecSpectro ) {
+    cout << "THcCoinTime module  Cannnot find Electron Arm = " << felecArmName.Data() << endl;
+    fStatus = kInitError;
+    return fStatus;
+  }
+    
+  fCoinDet = dynamic_cast<THcTrigDet*>
+    ( FindModule( fCoinDetName.Data(), "THcTrigDet"));
+  if( !fCoinDet ) {
+    cout << "THcCoinTime module  Cannnot find TrigDet = " << fCoinDetName.Data() << endl;
+    fStatus = kInitError;
+    return fStatus;
+  }
+  
+
+  
+  if( (fStatus=THaPhysicsModule::Init( run_time )) != kOK ) {
+    return fStatus;
+  }
+
+  return fStatus;
+}
+
+//_____________________________________________________________________________
+Int_t THcCoinTime::ReadDatabase( const TDatime& date )
+{
+  // Read database. Gets variable needed for CoinTime calculation
+
+
+  DBRequest list[]={
+    {"eHadCoinTime_Offset",  &eHad_CT_Offset, kDouble, 0, 1},   //coin time offset for ep coincidences
+
+    {"HMS_CentralPathLen",  &HMScentralPathLen, kDouble, 0, 1},
+    {"SHMS_CentralPathLen", &SHMScentralPathLen, kDouble, 0, 1},
+    {0}
+  };
+  
+  //Default values if not read from param file
+  eHad_CT_Offset = 0.0;
+
+  HMScentralPathLen = 22.0*100.;
+  SHMScentralPathLen = 18.1*100.;
+
+  
+  gHcParms->LoadParmValues((DBRequest*)&list, "");
+
+  return kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcCoinTime::DefineVariables( EMode mode )
+{
+
+  if( mode == kDefine && fIsSetup ) return kOK;
+  fIsSetup = ( mode == kDefine );
+
+  const RVarDef vars[] = {
+    {"epCoinTime_ROC1",    "ROC1 Corrected ep Coincidence Time",  "fROC1_epCoinTime"},
+    {"epCoinTime_ROC2",    "ROC2 Corrected ep Coincidence Time",  "fROC2_epCoinTime"},
+  
+    {"eKCoinTime_ROC1",    "ROC1 Corrected eK Coincidence Time",  "fROC1_eKCoinTime"},
+    {"eKCoinTime_ROC2",    "ROC2 Corrected eK Coincidence Time",  "fROC2_eKCoinTime"},
+    
+    {"ePiCoinTime_ROC1",    "ROC1 Corrected ePi Coincidence Time",  "fROC1_ePiCoinTime"},
+    {"ePiCoinTime_ROC2",    "ROC2 Corrected ePi Coincidence Time",  "fROC2_ePiCoinTime"},
+      
+    {"ePositronCoinTime_ROC1",    "ROC1 Corrected e-Positorn Coincidence Time",  "fROC1_ePosCoinTime"},
+    {"ePositronCoinTime_ROC2",    "ROC2 Corrected e-Positron Coincidence Time",  "fROC2_ePosCoinTime"},
+    
+    {"CoinTime_RAW_ROC1",    "ROC1 RAW Coincidence Time",  "fROC1_RAW_CoinTime"},
+    {"CoinTime_RAW_ROC2",    "ROC2 RAW Coincidence Time",  "fROC2_RAW_CoinTime"},
+    { 0 }
+  };
+
+  return DefineVarsFromList( vars, mode );
+
+}
+
+
+
+//_____________________________________________________________________________
+Int_t THcCoinTime::Process( const THaEvData& evdata )
+{
+  
+  if( !IsOK() || !gHaRun ) return -1;
+
+  //Declare track information objects for hadron/electron arm
+  THaTrackInfo* had_trkifo = fhadSpectro->GetTrackInfo();
+  THaTrackInfo* elec_trkifo = felecSpectro->GetTrackInfo();
+
+  if( !had_trkifo) cout << " no hadron track " << endl;
+  if( !elec_trkifo) cout << " no electron track " << endl;
+  //Check if the hadron/electron arm had a track
+  if( !had_trkifo || !had_trkifo->IsOK() ) return 1;
+  if( !elec_trkifo || !elec_trkifo->IsOK() ) return 1;
+
+  //Create THaTrack object for hadron/elec arms to get relevant golden track quantities
+  if (felecArmName=="H") {
+    theSHMSTrack =(fhadSpectro->GetGoldenTrack()); 
+    theHMSTrack = (felecSpectro->GetGoldenTrack());
+  } else{
+    theSHMSTrack =(felecSpectro->GetGoldenTrack()); 
+    theHMSTrack = (fhadSpectro->GetGoldenTrack());
+  }
+  //Gather relevant variables for Coincidence time calculation
+  lightSpeed = 29.9792; // in cm/ns
+
+  //Particle Masses (HardCoded)
+  elecMass =  0.510998/1000.0; // electron mass in GeV/c^2
+  positronMass =  0.510998/1000.0;
+  protonMass = 938.27208/1000.0; // proton mass in GeV/c^2	
+  kaonMass = 493.677/1000.0;    //charged kaon mass in GeV/c^2
+  pionMass = 139.570/1000.0;    //charged pion mass in GeV/c^2
+
+ 
+  //Check if there was a golden track in both arms
+  if (!theSHMSTrack || !theHMSTrack)
+    {
+      return 1;
+    }
+
+  //Check if Database is reading the correct elec-arm particle mass
+  if (felecSpectro->GetParticleMass() > 0.00052) return 1;
+   
+     
+      //SHMS arm
+      Double_t shms_xptar = theSHMSTrack->GetTTheta();     
+      Double_t shms_dP = theSHMSTrack->GetDp();            
+      Double_t SHMS_FPtime = theSHMSTrack->GetFPTime();    
+      
+      //HMS arm
+      Double_t hms_xfp = theHMSTrack->GetX();           
+      Double_t hms_xpfp = theHMSTrack->GetTheta();      
+      Double_t hms_ypfp = theHMSTrack->GetPhi();        
+      Double_t HMS_FPtime = theHMSTrack->GetFPTime();    
+      
+      //Get raw TDC Times for HMS/SHMS (3/4 trigger)
+      pTRIG1_rawTdcTime_ROC1 = fCoinDet->Get_pTRG1_ROC1_rawTdctime();  //TDC Channels (0.1 ns/Ch)
+      pTRIG4_rawTdcTime_ROC1 = fCoinDet->Get_pTRG4_ROC1_rawTdctime();
+      pTRIG1_rawTdcTime_ROC2 = fCoinDet->Get_pTRG1_ROC2_rawTdctime();
+      pTRIG4_rawTdcTime_ROC2 = fCoinDet->Get_pTRG4_ROC2_rawTdctime();
+
+	  DeltaSHMSpathLength = -0.11*atan2(shms_xptar,1)*1000 - 0.057*shms_dP;
+	  DeltaHMSpathLength = 12.462*hms_xpfp + 0.1138*hms_xpfp*hms_xfp - 0.0154*hms_xfp - 72.292*hms_xpfp*hms_xpfp - 0.0000544*hms_xfp*had_xfp - 116.52*hms_ypfp*hms_ypfp;
+          // default assume SHMS is electron arm
+	  Double_t ElecPathLength=SHMScentralPathLen +  DeltaSHMSpathLength;
+	  Double_t HadPathLength=HMScentralPathLen +  DeltaHMSpathLength;
+          elec_P = theSHMSTrack->GetP();              //electron golden track arm momentum
+          had_P = theHMSTrack->GetP();              //hadron golden track arm momentum
+          Int_t sign=-1;
+	  if (felecArmName=="H") {
+             ElecPathLength=HMScentralPathLen +  DeltaHMSpathLength;
+	     HadPathLength=SHMScentralPathLen +  DeltaSHMSpathLength;
+             elec_P = theHMSTrack->GetP();              //electron golden track arm momentum
+             had_P = theSHMSTrack->GetP();              //hadron golden track arm momentum
+             sign=1;
+	  }
+
+	  //beta calculations beta = v/c = p/E
+	  elecArm_BetaCalc = elec_P / sqrt(elec_P*elec_P + elecMass*elecMass);
+	  hadArm_BetaCalc_proton = had_P / sqrt(had_P*had_P + protonMass*protonMass);
+	  hadArm_BetaCalc_Kaon = had_P / sqrt(had_P*had_P + kaonMass*kaonMass);
+	  hadArm_BetaCalc_Pion = had_P / sqrt(had_P*had_P + pionMass*pionMass);	
+	  hadArm_BetaCalc_Positron = had_P / sqrt(had_P*had_P + positronMass*positronMass);
+
+
+	  //Coincidence Corrections
+	  elec_coinCorr = (ElecPathLength) / (lightSpeed * elecArm_BetaCalc );
+	  had_coinCorr_proton = (HadPathLength) / (lightSpeed * hadArm_BetaCalc_proton );
+	  had_coinCorr_Kaon =  (HadPathLength)/ (lightSpeed * hadArm_BetaCalc_Kaon );
+	  had_coinCorr_Pion =  (HadPathLength)/ (lightSpeed * hadArm_BetaCalc_Pion );
+	  had_coinCorr_Positron = (HadPathLength) / (lightSpeed * hadArm_BetaCalc_Positron );
+
+	  //Raw, Uncorrected Coincidence Time
+	  fROC1_RAW_CoinTime =  (pTRIG1_rawTdcTime_ROC1*0.1 + SHMS_FPtime) - (pTRIG4_rawTdcTime_ROC1*0.1 + HMS_FPtime);
+	  fROC2_RAW_CoinTime =  (pTRIG1_rawTdcTime_ROC2*0.1 + SHMS_FPtime) - (pTRIG4_rawTdcTime_ROC2*0.1 + HMS_FPtime);
+	  
+	  
+	  //Corrected Coincidence Time for ROC1/ROC2 (ROC1 Should be identical to ROC2)
+          // 
+	  //PROTON
+	  fROC1_epCoinTime = fROC1_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_proton) - eHad_CT_Offset; 
+	  fROC2_epCoinTime = fROC2_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_proton) - eHad_CT_Offset; 
+
+	  //KAON
+	  fROC1_eKCoinTime = fROC1_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_Kaon) - eHad_CT_Offset;
+	  fROC2_eKCoinTime = fROC2_RAW_CoinTime + sign*( elec_coinCorr-had_coinCorr_Kaon) - eHad_CT_Offset;
+	
+	  //PION
+	  fROC1_ePiCoinTime = fROC1_RAW_CoinTime + sign*( elec_coinCorr - had_coinCorr_Pion) - eHad_CT_Offset;	  
+	  fROC2_ePiCoinTime = fROC2_RAW_CoinTime + sign*( elec_coinCorr - had_coinCorr_Pion) - eHad_CT_Offset;
+
+	  //POSITRON
+	  fROC1_ePosCoinTime = fROC1_RAW_CoinTime + sign*( elec_coinCorr + had_coinCorr_Positron) - eHad_CT_Offset ;	  
+	  fROC2_ePosCoinTime = fROC2_RAW_CoinTime + sign*( elec_coinCorr + had_coinCorr_Positron) - eHad_CT_Offset;
+
+   
+  
+  
+  return 0;
+}
+
+//_____________________________________________________________________________
+
+ClassImp(THcCoinTime)
+////////////////////////////////////////////////////////////////////////////////

src/THcCoinTime.h

+#ifndef ROOT_THcCoinTime
+#define ROOT_THcCoinTime
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+// THcCoinTime Physics Module                                                //
+// Author: Carlos Yero                                                       //
+// Date: April 27, 2018                                                      //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THaEvData.h"
+#include "THaCutList.h"
+#include "VarDef.h"
+#include "VarType.h"
+#include "TClonesArray.h"
+
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+
+#include "THaPhysicsModule.h"
+#include "THcTrigDet.h" 
+#include "THcHodoscope.h"
+//#include "THaSpectrometer.h"
+#include "THcHallCSpectrometer.h"
+#include "THaTrack.h"
+
+class THcCoinTime : public THaPhysicsModule {
+public:
+  THcCoinTime( const char* name, const char* description, const char* hadArmName="", 
+	       const char* elecArmName="", const char* coinname="");
+
+  virtual ~THcCoinTime();
+
+
+  virtual EStatus Init( const TDatime& run_time );
+  virtual Int_t   Process( const THaEvData& );
+
+  void            Reset( Option_t* opt="" );
+  void            Clear( Option_t* opt="" );
+  
+
+ protected:
+
+  virtual Int_t ReadDatabase( const TDatime& date);
+  virtual Int_t  DefineVariables( EMode mode = kDefine );
+
+  // Data needed for adding coincidence time as a Leaf Variable
+  
+  TString     fHodName;		// Name of hodoscope
+  TString     fCoinDetName;         // Name of Coin Trigger
+  TString     fhadArmName;       //name of hadron arm
+  TString     felecArmName;     // name of electron arm
+
+
+  THcHallCSpectrometer* fhadSpectro;	// hadron Spectrometer object
+  THcHallCSpectrometer* felecSpectro;	// electron Spectrometer object
+  THcTrigDet* fCoinDet;                 //Coin Trigger detector object
+
+  THaTrack* theSHMSTrack;
+  THaTrack* theHMSTrack;
+
+  THcHodoscope* fHod;	                // Hodscope object
+ 
+  //-----Declare Variables used in HMS/SHMS Coin. time correction-----
+  Double_t lightSpeed;
+  Double_t elecMass;
+  Double_t positronMass;
+
+  //hadron masses (the e- could be in coincidence with any of the hadrons)
+  Double_t protonMass;
+  Double_t kaonMass;
+  Double_t pionMass;
+
+  Double_t eHad_CT_Offset;  //e-Hadron coin time Offset
+
+  Double_t SHMScentralPathLen;  
+  Double_t HMScentralPathLen;   
+
+  Double_t DeltaSHMSpathLength;
+  Double_t DeltaHMSpathLength;
+
+  Double_t fROC1_RAW_CoinTime;
+  Double_t fROC2_RAW_CoinTime;
+  
+  
+  Double_t fROC1_epCoinTime;
+  Double_t fROC2_epCoinTime;
+
+  Double_t fROC1_eKCoinTime;
+  Double_t fROC2_eKCoinTime;
+ 
+  Double_t fROC1_ePiCoinTime;
+  Double_t fROC2_ePiCoinTime;
+ 
+  Double_t fROC1_ePosCoinTime;   //electron-positron coin time 
+  Double_t fROC2_ePosCoinTime;
+  
+  Double_t elec_coinCorr;
+  Double_t elecArm_BetaCalc;
+  Double_t elec_hodFPtime;
+  
+  Double_t had_coinCorr_proton;
+  Double_t hadArm_BetaCalc_proton;
+  
+  Double_t had_coinCorr_Kaon;
+  Double_t hadArm_BetaCalc_Kaon;
+   
+  Double_t had_coinCorr_Pion;
+  Double_t hadArm_BetaCalc_Pion;
+  
+  Double_t had_coinCorr_Positron;
+  Double_t hadArm_BetaCalc_Positron;
+  
+
+  Double_t elec_P;     //electron golden track momentum
+  Double_t elec_dP;     //electron golden track delta-> (P-P0 / P0)
+  Double_t elec_xptar;    //electron golden track theta (xptar, :) 
+  Double_t elec_FPtime;   //electron focal plane time
+
+  Double_t had_P;     //hadron golden track momentum
+  Double_t had_xfp;      //hadron x-focal plane
+  Double_t had_xpfp;     //hadron xp focal plane
+  Double_t had_ypfp;     //hadron yp focal plane
+  Double_t had_FPtime;   //hadron focal plane time
+
+  //Raw trigger times pTrig1 (SHMS 3/4 trig) and pTrig4 (HMS 3/4 trig)
+  Int_t pTRIG1_rawTdcTime_ROC1;
+  Int_t pTRIG4_rawTdcTime_ROC1;
+  Int_t pTRIG1_rawTdcTime_ROC2;
+  Int_t pTRIG4_rawTdcTime_ROC2;
+
+  //--------------------------------------------------------------------
+
+  ClassDef(THcCoinTime,0) 	// Coincidence Time Module
+};
+
+#endif

src/THcTrigDet.cxx

@@ -221,8 +221,19 @@ void THcTrigDet::Clear(Option_t* opt) {
   };
 }
 
+//Added function to SET coincidence trigger times
+void THcTrigDet::SetCoinTrigTimes() 
+{ 
+  pTrig1_ROC1 = fTdcTimeRaw[fidx0];
+  pTrig4_ROC1 = fTdcTimeRaw[fidx1];
+  pTrig1_ROC2 = fTdcTimeRaw[fidx2];
+  pTrig4_ROC2 = fTdcTimeRaw[fidx3];
+
+}
+
 
 Int_t THcTrigDet::Decode(const THaEvData& evData) {
+    
   // Decode raw data for this event.
   Bool_t present = kTRUE;	// Don't suppress reference time warnings
   if(HaveIgnoreList()) {
@@ -289,10 +300,15 @@ Int_t THcTrigDet::Decode(const THaEvData& evData) {
     ++iHit;
   }
 
+  //Set raw Tdc coin. trigger times for pTRIG1/4
+  SetCoinTrigTimes();
+
   return 0;
 }
 
 
+
+
 void THcTrigDet::Setup(const char* name, const char* description) {
   // Prefix for parameters in `param` file.
   string kwPrefix = string(GetApparatus()->GetName()) + "_" + name;
@@ -348,11 +364,44 @@ Int_t THcTrigDet::ReadDatabase(const TDatime& date) {
   fAdcNames = vsplit(adcNames);
   fTdcNames = vsplit(tdcNames);
 
+  //default index values
+  fidx0 = 27;
+  fidx1 = 30;
+  fidx2 = 58;
+  fidx3 = 61;	
+
+  //Assign an index to coincidence trigger times strings
+  for (int i = 0; i <fNumTdc; i++)
+    {
+     
+      if(fTdcNames.at(i)=="pTRIG1_ROC1")
+	{
+	  fidx0 = i;
+	}
+      else if(fTdcNames.at(i)=="pTRIG4_ROC1")
+	{
+	  fidx1 = i;
+	}
+      else if(fTdcNames.at(i)=="pTRIG1_ROC2")
+	{
+	  fidx2 = i;
+	}
+      else if(fTdcNames.at(i)=="pTRIG4_ROC2")
+	{
+	  fidx3 = i;
+	}
+
+    }
+  
+  
+
   return kOK;
 }
 
 
 Int_t THcTrigDet::DefineVariables(THaAnalysisObject::EMode mode) {
+
+
   if (mode == kDefine && fIsSetup) return kOK;
   fIsSetup = (mode == kDefine);
 
@@ -461,6 +510,8 @@ Int_t THcTrigDet::DefineVariables(THaAnalysisObject::EMode mode) {
   for (int i=0; i<fNumTdc; ++i) {
     tdcTimeRawTitle.at(i) = fTdcNames.at(i) + "_tdcTimeRaw";
     tdcTimeRawVar.at(i) = TString::Format("fTdcTimeRaw[%d]", i);
+    
+
     RVarDef entry1 {
       tdcTimeRawTitle.at(i).Data(),
       tdcTimeRawTitle.at(i).Data(),

src/THcTrigDet.h

@@ -32,11 +32,32 @@ class THcTrigDet : public THaDetector, public THcHitList {
     virtual Bool_t IsIgnoreType(Int_t evtype) const;
     virtual Bool_t HaveIgnoreList() const;
     Int_t          End(THaRunBase* run);
+    
+    //Funtions to get raw TdcTimes 
+    Int_t Get_pTRG1_ROC1_rawTdctime() {return pTrig1_ROC1;}
+    Int_t Get_pTRG4_ROC1_rawTdctime() {return pTrig4_ROC1;}
+    Int_t Get_pTRG1_ROC2_rawTdctime() {return pTrig1_ROC2;}
+    Int_t Get_pTRG4_ROC2_rawTdctime() {return pTrig4_ROC2;}
 
   protected:
     void Setup(const char* name, const char* description);
     virtual Int_t ReadDatabase(const TDatime& date);
     virtual Int_t DefineVariables(EMode mode=kDefine);
+    
+    //Added for coin. time calculation
+    void SetCoinTrigTimes();
+
+    //indices to store pTrig string locations
+    Int_t fidx0;
+    Int_t fidx1;
+    Int_t fidx2;
+    Int_t fidx3;
+
+
+    Int_t pTrig1_ROC1;
+    Int_t pTrig4_ROC1;
+    Int_t pTrig1_ROC2;
+    Int_t pTrig4_ROC2;
 
     std::string fKwPrefix;