various photon fixes

src/libraries/HDDM/DEventSourceREST.h

@@ -16,6 +16,8 @@
 #include <JANA/JEventSource.h>
 #include <JANA/jerror.h>
 #include <JANA/JCalibration.h>
+#include <JANA/JCalibrationCCDB.h>
+#include <JANA/JCalibrationGeneratorCCDB.h>
 
 #include "hddm_r.hpp"
 
@@ -90,8 +92,6 @@ class DEventSourceREST:public JEventSource
    jerror_t Extract_DRFTime(hddm_r::HDDM *record,
                     JFactory<DRFTime>* factory);
 #endif
-   jerror_t Extract_DDIRCPmtHit(hddm_r::HDDM *record,
-                    JFactory<DDIRCPmtHit>* factory, JEventLoop* locEventLoop);
 
    void Get7x7ErrorMatrix(double mass, const double vec[5], const TMatrixFSym* C5x5, TMatrixFSym* loc7x7ErrorMatrix);
  private:
@@ -108,6 +108,9 @@ class DEventSourceREST:public JEventSource
 	int dDIRCMaxChannels;
 	enum dirc_status_state {GOOD, BAD, NOISY};
 	map<unsigned int, vector<vector<int>>> dDIRCChannelStatusMap; //unsigned int is run number
+
+	map<unsigned int, DVector2> dBeamCenterMap,dBeamDirMap;
+	map<unsigned int, double> dBeamZ0Map;
 
 	DFCALShower_factory *dFCALShowerFactory;
 	DBCALShower_factory_IU *dBCALShowerFactory;
@@ -118,6 +121,14 @@ class DEventSourceREST:public JEventSource
 
    std::ifstream *ifs;		// input hddm file ifstream
    hddm_r::istream *fin;	// provides hddm layer on top of ifstream
+
+   string REST_JANA_CALIB_CONTEXT = "";
+   JCalibrationGeneratorCCDB *calib_generator;
+
+   	map<unsigned int, JCalibration *> dJCalib_olds; //unsigned int is run number
+   	map<unsigned int, DTAGHGeometry *> dTAGHGeoms; //unsigned int is run number
+   	map<unsigned int, DTAGMGeometry *> dTAGMGeoms; //unsigned int is run number
+
 };
 
 #endif //_JEVENT_SOURCEREST_H_

src/libraries/PID/DBeamPhoton_factory.cc

@@ -88,6 +88,7 @@ jerror_t DBeamPhoton_factory::evnt(jana::JEventLoop *locEventLoop, uint64_t locE
     for (unsigned int ih=0; ih < tagh_hits.size(); ++ih)
     {
         if (!tagh_hits[ih]->has_fADC) continue; // Skip TDC-only hits (i.e. hits with no ADC info.)
+        if (!tagh_hits[ih]->has_TDC) continue;  // Skip fADC-only hits (i.e. hits with no TDC info.)
         DBeamPhoton *gamma = nullptr;
         for (unsigned int jh=0; jh < _data.size(); ++jh)
         {
@@ -125,7 +126,10 @@ void DBeamPhoton_factory::Set_BeamPhoton(DBeamPhoton* gamma, const DTAGMHit* hit
     gamma->setPosition(pos);
     gamma->setTime(hit->t);
     gamma->dCounter = hit->column;
-    gamma->dSystem = SYS_TAGM;
+    if(gamma->dCounter == 0)   // handle photons from simulation that miss tagger counters
+    	gamma->dSystem = SYS_NULL;
+    else
+    	gamma->dSystem = SYS_TAGM;
     gamma->AddAssociatedObject(hit);
 
 	auto locCovarianceMatrix = dResourcePool_TMatrixFSym->Get_SharedResource();
@@ -143,7 +147,10 @@ void DBeamPhoton_factory::Set_BeamPhoton(DBeamPhoton* gamma, const DTAGHHit* hit
     gamma->setPosition(pos);
     gamma->setTime(hit->t);
     gamma->dCounter = hit->counter_id;
-    gamma->dSystem = SYS_TAGH;
+    if(gamma->dCounter == 0)   // handle photons from simulation that miss tagger counters
+    	gamma->dSystem = SYS_NULL;
+    else
+	    gamma->dSystem = SYS_TAGH;
     gamma->AddAssociatedObject(hit);
 
 	auto locCovarianceMatrix = dResourcePool_TMatrixFSym->Get_SharedResource();

src/libraries/PID/DBeamPhoton_factory_MCGEN.cc

@@ -3,6 +3,9 @@
 using namespace std;
 
 #include "DBeamPhoton_factory_MCGEN.h"
+#include "TAGGER/DTAGHGeometry.h"
+#include "TAGGER/DTAGMGeometry.h"
+
 using namespace jana;
 
 //------------------
@@ -62,9 +65,30 @@ jerror_t DBeamPhoton_factory_MCGEN::evnt(jana::JEventLoop *locEventLoop, uint64_
 		}
 	}
 
+	// extract the TAGH geometry
+   vector<const DTAGHGeometry*> taghGeomVect;
+   eventLoop->Get(taghGeomVect);
+   if (taghGeomVect.empty())
+      return NOERROR;
+   const DTAGHGeometry* taghGeom = taghGeomVect[0];
+
+   // extract the TAGM geometry
+   vector<const DTAGMGeometry*> tagmGeomVect;
+   eventLoop->Get(tagmGeomVect);
+   if (tagmGeomVect.empty())
+      return NOERROR;
+   const DTAGMGeometry* tagmGeom = tagmGeomVect[0];
+
+
 	//Photon is NOT TAGGED //Create a beam object from the DMCReaction
 	auto *locBeamPhoton = new DBeamPhoton;
 	*(DKinematicData*)locBeamPhoton = locMCReactions[0]->beam;
+	if(tagmGeom->E_to_column(locBeamPhoton->energy(), locBeamPhoton->dCounter))
+		locBeamPhoton->dSystem = SYS_TAGM;
+	else if(taghGeom->E_to_counter(locBeamPhoton->energy(), locBeamPhoton->dCounter))
+		locBeamPhoton->dSystem = SYS_TAGH;
+	else
+		locBeamPhoton->dSystem = SYS_NULL;
 	_data.push_back(locBeamPhoton);
 
 	return NOERROR;

src/libraries/PID/DBeamPhoton_factory_TRUTH.cc

@@ -51,7 +51,10 @@ jerror_t DBeamPhoton_factory_TRUTH::evnt(jana::JEventLoop *locEventLoop, uint64_
 		gamma->setMomentum(mom);
 		gamma->setPosition(pos);
 		gamma->setTime(tagm_hits[ih]->t);
-	    gamma->dSystem = SYS_TAGM;
+		if(gamma->dCounter == 0)   // handle photons from simulation that miss tagger counters
+			gamma->dSystem = SYS_NULL;
+		else
+		    gamma->dSystem = SYS_TAGM;
 		gamma->dCounter = tagm_hits[ih]->column;
 		gamma->AddAssociatedObject(tagm_hits[ih]);
 		_data.push_back(gamma);
@@ -67,7 +70,10 @@ jerror_t DBeamPhoton_factory_TRUTH::evnt(jana::JEventLoop *locEventLoop, uint64_
 		gamma->setMomentum(mom);
 		gamma->setPosition(pos);
 		gamma->setTime(tagh_hits[ih]->t);
-	    gamma->dSystem = SYS_TAGH;
+		if(gamma->dCounter == 0)   // handle photons from simulation that miss tagger counters
+			gamma->dSystem = SYS_NULL;
+		else
+		    gamma->dSystem = SYS_TAGH;
 		gamma->dCounter = tagh_hits[ih]->counter_id;
 		gamma->AddAssociatedObject(tagh_hits[ih]);
 		_data.push_back(gamma);

src/libraries/TAGGER/DTAGHGeometry.cc

@@ -3,25 +3,77 @@
 // Created: Sat Jul 5 10:18:56 EST 2014
 // Creator: jonesrt on gluey.phys.uconn.edu
 //
+//  10/19/2019 A.S 
+//
+//  Modify calculation of the photon beam energy to account 
+//  for the fact that the energy of bremsstrahlung electrons detected
+//  by each tagger counter does not depend on the electron beam energy.
+//  The photon beam energy E_gamma has to be computed as
+//  
+//     E_gamma = R * E_endpoint_calib  +  DE,  where
+//     DE = E_endpoint - E_endpoint_calib
+//
+// 
 
 #include <stdlib.h>
 #include <iostream>
 #include <map>
 
-#include <JANA/JApplication.h>
-#include <JANA/JEvent.h>
 #include "DTAGHGeometry.h"
 
 const unsigned int DTAGHGeometry::kCounterCount = 274;
 
+// Only print messages for one thread whenever run number change
+static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER;
+static set<int> runs_announced;
+
+
 //---------------------------------
 // DTAGHGeometry    (Constructor)
 //---------------------------------
 DTAGHGeometry::DTAGHGeometry(JEventLoop *loop)
+{
+	// keep track of which runs we print out messages for
+	int32_t runnumber = loop->GetJEvent().GetRunNumber();
+	pthread_mutex_lock(&print_mutex);
+	bool print_messages = false;
+	if(runs_announced.find(runnumber) == runs_announced.end()){
+	  print_messages = true;
+	  runs_announced.insert(runnumber);
+	}
+	pthread_mutex_unlock(&print_mutex);
+
+	JEvent &event = loop->GetJEvent();
+	DApplication* dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
+	JCalibrationCCDB *jcalib =  dynamic_cast<JCalibrationCCDB*>( dapp->GetJCalibration(event.GetRunNumber()) );
+
+	Initialize(jcalib, print_messages);
+}
+
+//---------------------------------
+// DTAGHGeometry    (Constructor)
+//---------------------------------
+DTAGHGeometry::DTAGHGeometry(JCalibration *jcalib, int32_t runnumber)
+{
+	pthread_mutex_lock(&print_mutex);
+	bool print_messages = false;
+	if(runs_announced.find(runnumber) == runs_announced.end()){
+	  print_messages = true;
+	  runs_announced.insert(runnumber);
+	}
+	pthread_mutex_unlock(&print_mutex);
+
+	Initialize(jcalib, print_messages);
+}
+
+//---------------------------------
+// Initialize
+//---------------------------------
+void DTAGHGeometry::Initialize(JCalibration *jcalib, bool print_messages)
 {
    /* read tagger set endpoint energy from calibdb */
    std::map<string,double> result1;
-   loop->GetCalib("/PHOTON_BEAM/endpoint_energy", result1);
+   jcalib->Get("/PHOTON_BEAM/endpoint_energy", result1);
    if (result1.find("PHOTON_BEAM_ENDPOINT_ENERGY") == result1.end()) {
       std::cerr << "Error in DTAGHGeometry constructor: "
                 << "failed to read photon beam endpoint energy "
@@ -34,6 +86,63 @@ DTAGHGeometry::DTAGHGeometry(JEventLoop *loop)
 
    /* read hodoscope counter energy bounds from calibdb */
    std::vector<std::map<string,double> > result2;
+   jcalib->Get("/PHOTON_BEAM/hodoscope/scaled_energy_range", result2);
+   if (result2.size() != kCounterCount) {
+      jerr << "Error in DTAGHGeometry constructor: "
+           << "failed to read photon beam scaled_energy_range table "
+           << "from calibdb at /PHOTON_BEAM/hodoscope/scaled_energy_range" << std::endl;
+      for (unsigned int i=0; i <= TAGH_MAX_COUNTER; ++i) {
+         m_counter_xlow[i] = 0;
+         m_counter_xhigh[i] = 0;
+      }
+   }
+   else {
+      for (unsigned int i=0; i < result2.size(); ++i) {
+	int ctr = (result2[i])["counter"];
+         m_counter_xlow[ctr] = (result2[i])["xlow"];
+         m_counter_xhigh[ctr] = (result2[i])["xhigh"];
+      }
+   }
+
+   int status = 0;
+   m_endpoint_energy_calib_GeV = 0.;
+
+   std::map<string,double> result3;
+   status = jcalib->Get("/PHOTON_BEAM/hodoscope/endpoint_calib",result3);
+
+
+   if(!status){
+     if (result3.find("TAGGER_CALIB_ENERGY") == result3.end()) {
+       std::cerr << "Error in DTAGHGeometry constructor: "
+		 <<  "failed to read  endpoint_calib field "
+		 <<  "from /PHOTON_BEAM/hodoscope/endpoint_calib table" << std::endl;
+
+     } else {
+       m_endpoint_energy_calib_GeV  = result3["TAGGER_CALIB_ENERGY"];
+
+	   if(print_messages)
+       	  jout << " Correct Beam Photon Energy (TAGH) = " << m_endpoint_energy_calib_GeV << " (GeV)" << std::endl;
+
+     }
+   }
+
+
+   /*
+   // read tagger set endpoint energy from calibdb 
+   std::map<string,double> result1;
+   loop->GetCalib("/PHOTON_BEAM/endpoint_energy", result1);
+   if (result1.find("PHOTON_BEAM_ENDPOINT_ENERGY") == result1.end()) {
+      std::cerr << "Error in DTAGHGeometry constructor: "
+                << "failed to read photon beam endpoint energy "
+                << "from calibdb at /PHOTON_BEAM/endpoint_energy" << std::endl;
+      m_endpoint_energy_GeV = 0;
+   }
+   else {
+      m_endpoint_energy_GeV = result1["PHOTON_BEAM_ENDPOINT_ENERGY"];
+   }
+
+   // read hodoscope counter energy bounds from calibdb 
+   std::vector<std::map<string,double> > result2;
    loop->GetCalib("/PHOTON_BEAM/hodoscope/scaled_energy_range", result2);
    if (result2.size() != kCounterCount) {
       jerr << "Error in DTAGHGeometry constructor: "
@@ -46,40 +155,87 @@ DTAGHGeometry::DTAGHGeometry(JEventLoop *loop)
    }
    else {
       for (unsigned int i=0; i < result2.size(); ++i) {
-         int ctr = (result2[i])["counter"];
+	int ctr = (result2[i])["counter"];
          m_counter_xlow[ctr] = (result2[i])["xlow"];
          m_counter_xhigh[ctr] = (result2[i])["xhigh"];
       }
    }
+
+   int status = 0;
+   m_endpoint_energy_calib_GeV = 0.;
+   
+   std::map<string,double> result3;
+   status = loop->GetCalib("/PHOTON_BEAM/hodoscope/endpoint_calib",result3);
+   
+   
+   if(!status){
+     if (result3.find("TAGGER_CALIB_ENERGY") == result3.end()) {
+       std::cerr << "Error in DTAGHGeometry constructor: "
+		 <<  "failed to read  endpoint_calib field "
+		 <<  "from /PHOTON_BEAM/hodoscope/endpoint_calib table" << std::endl;
+       
+     } else {
+       m_endpoint_energy_calib_GeV  = result3["TAGGER_CALIB_ENERGY"];
+
+	   if(print_messages)
+       	  jout << " Correct Beam Photon Energy (TAGH) = " << m_endpoint_energy_calib_GeV << " (GeV)" << std::endl;
+
+     }
+   }
+   */
+
 }
 
 DTAGHGeometry::~DTAGHGeometry() { }
 
+
 bool DTAGHGeometry::E_to_counter(double E, unsigned int &counter) const
-{
-   double x = E/m_endpoint_energy_GeV;
-   for (counter=1; counter <= kCounterCount; ++counter) {
-      if ( x >= m_counter_xlow[counter] &&
-           x <= m_counter_xhigh[counter] )
+{  
+  for (counter = 1; counter <= kCounterCount; ++counter) {     
+
+    double Emin = getElow(counter);
+    double Emax = getEhigh(counter);
+
+    if ( E >= Emin &&  E <= Emax )
       {
-         return true;
+	return true;
       }
-   }
-   return false;
+  }
+  return false;
 }
 
+
 double DTAGHGeometry::getElow(unsigned int counter) const
 {
-   if (counter > 0 && counter <= kCounterCount)
+  if (counter > 0 && counter <= kCounterCount){
+    if(m_endpoint_energy_calib_GeV  > 0){
+
+      double delta_E  =  m_endpoint_energy_GeV  -  m_endpoint_energy_calib_GeV;
+      double Emin     =  m_counter_xlow[counter]*m_endpoint_energy_calib_GeV  +  delta_E;
+
+      return Emin;
+
+    } else
       return m_endpoint_energy_GeV * m_counter_xlow[counter];
-   else
-      return 0;
+
+  }  else
+    return 0;
 }
 
+
 double DTAGHGeometry::getEhigh(unsigned int counter) const
 {
-   if (counter > 0 && counter <= kCounterCount)
+  if (counter > 0 && counter <= kCounterCount){
+    if(m_endpoint_energy_calib_GeV  > 0){
+
+      double delta_E  =  m_endpoint_energy_GeV  -  m_endpoint_energy_calib_GeV;
+      double Emax     =  m_counter_xhigh[counter]*m_endpoint_energy_calib_GeV  +  delta_E;
+
+      return Emax;
+
+    } else    
       return m_endpoint_energy_GeV * m_counter_xhigh[counter];
-   else
-      return 0;
+
+  } else
+    return 0;
 }