Merge branch 'master' into helicity_tools

src/BMS/Makefile.DANA

@@ -6,4 +6,4 @@ include $(HALLD_HOME)/src/BMS/Makefile.CURL
 
 ADDITIONAL_MODULES += DANA ANALYSIS PID TAGGER TRACKING START_COUNTER
 ADDITIONAL_MODULES += CERE DIRC CDC TRIGGER
-ADDITIONAL_MODULES += FDC TOF BCAL FCAL CCAL HDGEOMETRY HDDM JANA
+ADDITIONAL_MODULES += FDC TOF BCAL FCAL ECAL CCAL HDGEOMETRY HDDM JANA

src/SBMS/SConstruct.plugin

@@ -24,6 +24,18 @@
 # diligent of building for multiple platforms using the same source.
 #
 
+from __future__ import division
+from __future__ import absolute_import
+from __future__ import print_function
+from __future__ import unicode_literals
+
+from builtins import map
+from builtins import str
+from builtins import open
+from builtins import int
+from future import standard_library
+standard_library.install_aliases()
+
 import os
 import sys
 import subprocess
@@ -32,7 +44,7 @@ import glob
 # Get HALLD_RECON_HOME environment variable, verifying it is set
 halld_home = os.getenv('HALLD_RECON_HOME')
 if(halld_home == None):
-	print 'HALLD_RECON_HOME environment variable not set!'
+	print('HALLD_RECON_HOME environment variable not set!')
 	exit(-1)
 
 # Get HALLD_MY if it exists. Otherwise use HALLD_RECON_HOME
@@ -51,7 +63,8 @@ SHOWBUILD = ARGUMENTS.get('SHOWBUILD', 0)
 osname = os.getenv('BMS_OSNAME', 'build')
 
 # Get architecture name
-arch = subprocess.Popen(["uname"], stdout=subprocess.PIPE).communicate()[0].strip()
+arch_out = ROOT_CFLAGS = subprocess.Popen(["uname"], stdout=subprocess.PIPE).communicate()[0].strip()
+arch = str(arch_out, 'utf-8')
 
 # Setup initial environment
 installdir = "%s/%s" %(halld_my, osname)

src/SBMS/sbms.py

@@ -641,7 +641,7 @@ def AddDANA(env):
 	AddTensorflowLite(env) # optional. Used for CPP pi/mu classification
 	DANA_LIBS  = "DANA ANALYSIS KINFITTER PID TAGGER TRACKING START_COUNTER"
 	DANA_LIBS += " CERE DIRC CDC TRIGGER PAIR_SPECTROMETER RF TRD"
-	DANA_LIBS += " FDC TOF BCAL FCAL CCAL TPOL HDGEOMETRY TTAB FMWPC TAC"
+	DANA_LIBS += " FDC TOF BCAL FCAL ECAL CCAL TPOL HDGEOMETRY TTAB FMWPC TAC"
 	DANA_LIBS += " DAQ JANA EVENTSTORE"
 	DANA_LIBS += " expat gfortran" 
 	env.PrependUnique(LIBS = DANA_LIBS.split())

src/SBMS/sbms_setenv.py

@@ -86,6 +86,13 @@ def mk_setenv_csh(env):
 	str += 'setenv HDDS_HOME %s\n' % os.getenv('HDDS_HOME', '$HOME/hdds')
 	str += '\n'
 
+	# HDDM
+	str += '# HDDM\n'
+	str += 'setenv HDDM_HOME %s\n' % os.getenv('HDDM_HOME', '$HOME/hddm')
+	str += 'setenv PATH ${HDDM_HOME}/bin:${PATH}\n'
+	str += 'setenv HDDM_DIR $HDDM_HOME\n'
+	str += '\n'
+
 	# JANA
 	str += '# JANA\n'
 	str += 'setenv JANA_HOME %s\n' % os.getenv('JANA_HOME', '$HOME/jana')
@@ -259,6 +266,13 @@ def mk_setenv_bash(env):
 	str += 'export HDDS_HOME=%s\n' % os.getenv('HDDS_HOME', '$HOME/hdds')
 	str += '\n'
 
+	# HDDM
+	str += '# HDDM\n'
+	str += 'export HDDM_HOME=%s\n' % os.getenv('HDDM_HOME', '$HOME/hddm')
+	str += 'export PATH=${HDDM_HOME}/bin:${PATH}\n'
+	str += 'export HDDM_DIR=$HDDM_HOME\n'
+	str += '\n'
+
 	# JANA
 	str += '# JANA\n'
 	str += 'export JANA_HOME=%s\n' % os.getenv('JANA_HOME', '$HOME/jana')

src/libraries/ANALYSIS/DAnalysisUtilities.cc

@@ -626,6 +626,35 @@ void DAnalysisUtilities::Get_UnusedNeutralParticles(JEventLoop* locEventLoop, co
 	}
 }
 
+void DAnalysisUtilities::Get_UnusedTOFPoints(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, vector<const DTOFPoint*>& locUnusedTOFPoints) const
+{
+	locUnusedTOFPoints.clear();
+	locEventLoop->Get(locUnusedTOFPoints);
+
+	vector<const DParticleComboStep*> locParticleComboSteps = locParticleCombo->Get_ParticleComboSteps();
+	for(size_t loc_icombo = 0; loc_icombo < locParticleComboSteps.size(); ++loc_icombo) {
+		const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_icombo);
+		for(size_t loc_ipart = 0; loc_ipart < locParticleComboStep->Get_NumFinalParticles(); ++loc_ipart) {
+
+			auto locParticle = locParticleComboStep->Get_FinalParticle_Measured(loc_ipart);
+			if(locParticle == nullptr || locParticle->charge() == 0) continue; 
+
+			const DChargedTrack* locChargedTrack = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(loc_ipart));
+			const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_Hypothesis(locParticle->PID());
+
+			if(locChargedTrackHypothesis->Get_TOFHitMatchParams() == nullptr) continue;
+			const DTOFPoint *locTOFPoint = locChargedTrackHypothesis->Get_TOFHitMatchParams()->dTOFPoint;
+
+			for(auto locIterator = locUnusedTOFPoints.begin(); locIterator != locUnusedTOFPoints.end();) {
+				if(locTOFPoint != *locIterator)
+					++locIterator;
+				else
+					locIterator = locUnusedTOFPoints.erase(locIterator);
+			}
+		}
+	}
+}
+
 void DAnalysisUtilities::Get_ThrownParticleSteps(JEventLoop* locEventLoop, deque<pair<const DMCThrown*, deque<const DMCThrown*> > >& locThrownSteps) const
 {
  	vector<const DMCThrown*> locMCThrowns;

src/libraries/ANALYSIS/DAnalysisUtilities.h

@@ -72,6 +72,7 @@ class DAnalysisUtilities : public JObject
 
 		void Get_UnusedNeutralShowers(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, vector<const DNeutralShower*>& locUnusedNeutralShowers) const;
 		void Get_UnusedNeutralParticles(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, vector<const DNeutralParticle*>& locUnusedNeutralParticles) const;
+		void Get_UnusedTOFPoints(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo, vector<const DTOFPoint*>& locUnusedTOFPoints) const;
 
 		void Get_ThrownParticleSteps(JEventLoop* locEventLoop, deque<pair<const DMCThrown*, deque<const DMCThrown*> > >& locThrownSteps) const;
 		bool Are_ThrownPIDsSameAsDesired(JEventLoop* locEventLoop, const deque<Particle_t>& locDesiredPIDs, Particle_t locMissingPID = Unknown) const;

src/libraries/ANALYSIS/DEventWriterROOT.cc

@@ -167,6 +167,8 @@ void DEventWriterROOT::Create_DataTree(const DReaction* locReaction, JEventLoop*
 	locBranchRegister.Register_Single<UInt_t>("RunNumber");
 	locBranchRegister.Register_Single<ULong64_t>("EventNumber");
 	locBranchRegister.Register_Single<UInt_t>("L1TriggerBits");
+	locBranchRegister.Register_Single<Double_t>("L1BCALEnergy");
+	locBranchRegister.Register_Single<Double_t>("L1FCALEnergy");
 
 	//create X4_Production
 	locBranchRegister.Register_Single<TLorentzVector>("X4_Production");
@@ -1344,6 +1346,9 @@ void DEventWriterROOT::Fill_DataTree(JEventLoop* locEventLoop, const DReaction*
 	locTreeFillData->Fill_Single<UInt_t>("RunNumber", locEventLoop->GetJEvent().GetRunNumber());
 	locTreeFillData->Fill_Single<ULong64_t>("EventNumber", locEventLoop->GetJEvent().GetEventNumber());
 	locTreeFillData->Fill_Single<UInt_t>("L1TriggerBits", locTrigger->Get_L1TriggerBits());
+	locTreeFillData->Fill_Single<Double_t>("L1BCALEnergy", locTrigger->Get_GTP_BCALEnergy());
+	locTreeFillData->Fill_Single<Double_t>("L1FCALEnergy", locTrigger->Get_GTP_FCALEnergy());
+
 
 	//PRODUCTION X4
 	DLorentzVector locProductionX4 = locVertex->dSpacetimeVertex;
@@ -2409,7 +2414,7 @@ void DEventWriterROOT::Fill_KinFitData(DTreeFillData* locTreeFillData, JEventLoo
         //Look at the decay products:
         //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
         vector<const DKinematicData*> finalParticlesMeasured = locParticleCombos.at(iPC)->Get_FinalParticles_Measured(locReaction, d_AllCharges);
-        vector<const DKinematicData*> finalParticles = locParticleCombos.at(iPC)->Get_FinalParticles(locReaction, d_AllCharges);
+        vector<const DKinematicData*> finalParticles = locParticleCombos.at(iPC)->Get_FinalParticles(locReaction); //, d_AllCharges);
 
         vector<const JObject*> finalParticleObjects = locParticleCombos.at(iPC)->Get_FinalParticle_SourceObjects(d_AllCharges);
         nMeasuredFinalParticles = finalParticlesMeasured.size();

src/libraries/ANALYSIS/DHistogramActions_Independent.cc

@@ -4032,7 +4032,19 @@ bool DHistogramAction_TriggerStudies::Perform_Action(JEventLoop* locEventLoop, c
 	locEventLoop->GetSingle(locTrigger);
 	if(locTrigger == nullptr)
 		return true;
+
+	// allow for this histogram to be called for a particular reaction, and then only plot it for 
+	// events that pass some reasonable kinematic fit cut
+ 	if(dKinFitCLCut >= 0.) {
+		const DKinFitResults* locKinFitResults = locParticleCombo->Get_KinFitResults();
+		if(locKinFitResults == NULL)
+			return true;
+		double locConfidenceLevel = locKinFitResults->Get_ConfidenceLevel();
+		if(locConfidenceLevel < dKinFitCLCut)
+			return true;
+	}
 
+
 
 	//FILL HISTOGRAMS
 	//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

src/libraries/ANALYSIS/DHistogramActions_Independent.h

@@ -993,21 +993,23 @@ class DHistogramAction_TrackMultiplicity : public DAnalysisAction
 class DHistogramAction_TriggerStudies : public DAnalysisAction
 {
 	public:
-		DHistogramAction_TriggerStudies(const DReaction* locReaction, string locActionUniqueString = "") : 
+		DHistogramAction_TriggerStudies(const DReaction* locReaction, string locActionUniqueString = "", double locKinFitCLCut=-1.) : 
 		DAnalysisAction(locReaction, "Hist_TriggerStudies", false, locActionUniqueString),
-		dBCALBins(240), dFCALBins(240), dMaxBCALEnergy(12.), dMaxFCALEnergy(12.) { }
+		dBCALBins(240), dFCALBins(240), dMaxBCALEnergy(12.), dMaxFCALEnergy(12.), dKinFitCLCut(locKinFitCLCut) { }
 
 		DHistogramAction_TriggerStudies(string locActionUniqueString) : 
 		DAnalysisAction(NULL, "Hist_TriggerStudies", false, ""),
-		dBCALBins(240), dFCALBins(240), dMaxBCALEnergy(12.), dMaxFCALEnergy(12.) { }
+		dBCALBins(240), dFCALBins(240), dMaxBCALEnergy(12.), dMaxFCALEnergy(12.), dKinFitCLCut(-1.)  { }
 
 		DHistogramAction_TriggerStudies(void) : 
 		DAnalysisAction(NULL, "Hist_TriggerStudies", false, ""),
-		dBCALBins(240), dFCALBins(240), dMaxBCALEnergy(12.), dMaxFCALEnergy(12.) { }
+		dBCALBins(240), dFCALBins(240), dMaxBCALEnergy(12.), dMaxFCALEnergy(12.), dKinFitCLCut(-1.)  { }
 
 		int dBCALBins, dFCALBins;
 		double dMaxBCALEnergy;
 		double dMaxFCALEnergy;
+
+		double dKinFitCLCut;
 
 		void Initialize(JEventLoop* locEventLoop);
 		void Run_Update(JEventLoop* locEventLoop){}

src/libraries/ANALYSIS/DSourceComboP4Handler.cc

@@ -134,13 +134,13 @@ void DSourceComboP4Handler::Define_DefaultCuts(void)
 	dInvariantMassCuts.emplace(omega, std::make_pair(0.4, 1.2));
 	dInvariantMassCuts.emplace(EtaPrime, std::make_pair(0.6, 1.3));
 	dInvariantMassCuts.emplace(phiMeson, std::make_pair(0.8, 1.2));
-	dInvariantMassCuts.emplace(D0, std::make_pair(1.8, 1.92));
-	dInvariantMassCuts.emplace(AntiD0, std::make_pair(1.8, 1.92));
-	dInvariantMassCuts.emplace(DPlus, std::make_pair(1.8, 1.92));
-	dInvariantMassCuts.emplace(DMinus, std::make_pair(1.8, 1.92));
-	dInvariantMassCuts.emplace(Dstar0, std::make_pair(1.9, 2.1));
-	dInvariantMassCuts.emplace(DstarPlus, std::make_pair(1.9, 2.1));
-	dInvariantMassCuts.emplace(DstarMinus, std::make_pair(1.9, 2.1));
+	dInvariantMassCuts.emplace(D0, std::make_pair(1.7, 2.0));
+	dInvariantMassCuts.emplace(AntiD0, std::make_pair(1.7, 2.0));
+	dInvariantMassCuts.emplace(DPlus, std::make_pair(1.7, 2.0));
+	dInvariantMassCuts.emplace(DMinus, std::make_pair(1.7, 2.0));
+	dInvariantMassCuts.emplace(Dstar0, std::make_pair(1.8, 2.2));
+	dInvariantMassCuts.emplace(DstarPlus, std::make_pair(1.8, 2.2));
+	dInvariantMassCuts.emplace(DstarMinus, std::make_pair(1.8, 2.2));
 //	dInvariantMassCuts.emplace(Jpsi, std::make_pair(2.7, 3.5));
 
 	//INVARIANT MASS CUTS: BARYONS
@@ -157,7 +157,8 @@ void DSourceComboP4Handler::Define_DefaultCuts(void)
 	dInvariantMassCuts.emplace(AntiXi0, dInvariantMassCuts[XiMinus]);
 	dInvariantMassCuts.emplace(AntiXiPlus, dInvariantMassCuts[XiMinus]);
 	dInvariantMassCuts.emplace(OmegaMinus, std::make_pair(1.32, 2.22));
-	dInvariantMassCuts.emplace(Lambda_c, std::make_pair(2.0, 2.6));
+	dInvariantMassCuts.emplace(Lambda_c, std::make_pair(2.1, 2.5));
+	dInvariantMassCuts.emplace(Sigma_cPlusPlus, std::make_pair(2.2, 2.7));
 
 
 	//DEFAULT MISSING MASS SQUARED CUT FUNCTION

src/libraries/ANALYSIS/DTreeInterface.h

@@ -213,7 +213,9 @@ template <typename DType> inline void DTreeInterface::Create_Branch_TObject(stri
 		return; //already created
 
 	dMemoryMap_TObject[locBranchName] = (TObject*)(new DType());
-	dTree->Branch(locBranchName.c_str(), (DType**)&(dMemoryMap_TObject[locBranchName]), 32000, 0); //0: don't split
+	TObject **objptr = new TObject*;
+	*objptr = dMemoryMap_TObject[locBranchName];
+	dTree->Branch(locBranchName.c_str(), objptr, 32000, 0); //0: don't split
 }
 
 template <typename DType> inline void DTreeInterface::Create_Branch_FundamentalArray(string locBranchName, string locArraySizeString, unsigned int locInitialSize)
@@ -232,7 +234,9 @@ template <typename DType> inline void DTreeInterface::Create_Branch_ClonesArray(
 		return; //already created
 
 	dMemoryMap_ClonesArray[locBranchName] = new TClonesArray(DType::Class()->GetName(), locInitialSize);
-	dTree->Branch(locBranchName.c_str(), &(dMemoryMap_ClonesArray[locBranchName]), 32000, 0); //0: don't split
+	TClonesArray **arrptr = new TClonesArray*;
+	*arrptr = dMemoryMap_ClonesArray[locBranchName];
+	dTree->Branch(locBranchName.c_str(), arrptr, 32000, 0); //0: don't split
 }
 
 /************************************************************** TEMPLATE SPECIALIZATIONS **************************************************************/

src/libraries/BCAL/DBCALHit_factory.cc

@@ -61,6 +61,7 @@ jerror_t DBCALHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber)
    vector<double> raw_ADC_timing_offsets;
    vector<double> raw_channel_global_offset;
    vector<double> raw_tdiff_u_d;
+   vector<double> raw_bad_channels;
 
     if(print_messages) jout << "In DBCALHit_factory, loading constants..." << endl;
 
@@ -105,6 +106,8 @@ jerror_t DBCALHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber)
        jout << "Error loading /BCAL/channel_global_offset !" << endl;
    if(eventLoop->GetCalib("/BCAL/tdiff_u_d", raw_tdiff_u_d))
        jout << "Error loading /BCAL/tdiff_u_d !" << endl;
+   if(eventLoop->GetCalib("/BCAL/bad_channels", raw_bad_channels))
+       jout << "Error loading /BCAL/bad_channels !" << endl;
 
    if (PRINTCALIBRATION) jout << "DBCALHit_factory >> raw_gains" << endl;
    FillCalibTable(gains, raw_gains);
@@ -116,6 +119,8 @@ jerror_t DBCALHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber)
    FillCalibTableShort(channel_global_offset, raw_channel_global_offset);
    if (PRINTCALIBRATION) jout << "DBCALHit_factory >> raw_tdiff_u_d" << endl;
    FillCalibTableShort(tdiff_u_d, raw_tdiff_u_d);
+   if (PRINTCALIBRATION) jout << "DBCALHit_factory >> raw_bad_channels" << endl;
+   FillCalibTable(bad_channels, raw_bad_channels);
 
    std::vector<std::map<string,double> > saturation_ADC_pars;
    if(eventLoop->GetCalib("/BCAL/ADC_saturation", saturation_ADC_pars))
@@ -167,6 +172,10 @@ jerror_t DBCALHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
    for(unsigned int i=0; i<digihits.size(); i++){
       const DBCALDigiHit *digihit = digihits[i];
 
+      // throw away hits from bad or noisy channels
+      int quality = GetConstant(bad_channels,digihit);
+      if ( quality > 0 ) continue;
+
       // Error checking for pre-Fall 2016 firmware
       if(digihit->datasource == 1) {
           // There is a slight difference between Mode 7 and 8 data

src/libraries/BCAL/DBCALHit_factory.h

@@ -58,6 +58,7 @@ class DBCALHit_factory:public jana::JFactory<DBCALHit>{
 		bcal_digi_constants_t ADC_timing_offsets;
         bcal_digi_constants_t channel_global_offset;
         bcal_digi_constants_t tdiff_u_d;
+        bcal_digi_constants_t bad_channels;
 
 		const int GetCalibIndex( int module, int layer, int sector ) const {
 			return BCAL_NUM_LAYERS*BCAL_NUM_SECTORS*(module-1) + BCAL_NUM_SECTORS*(layer-1) + (sector-1);

src/libraries/CCAL/CCAL_init.cc

@@ -15,13 +15,13 @@ using namespace jana;
 jerror_t CCAL_init(JEventLoop *loop)
 {
 	/// Create and register CCAL data factories
-  	loop->AddFactory(new JFactory<DCCALDigiHit>());
-  	loop->AddFactory(new JFactory<DCCALRefDigiHit>());
+	loop->AddFactory(new JFactory<DCCALDigiHit>());
+	loop->AddFactory(new JFactory<DCCALRefDigiHit>());
 	loop->AddFactory(new DCCALHit_factory());
 	loop->AddFactory(new DCCALShower_factory());
 	loop->AddFactory(new JFactory<DCCALHit>("TRUTH"));
 	loop->AddFactory(new DCCALTruthShower_factory());
 	loop->AddFactory(new DCCALGeometry_factory());
-
+	
 	return NOERROR;
 }

src/libraries/CCAL/DCCALGeometry.cc

@@ -18,23 +18,23 @@ using namespace std;
 DCCALGeometry::DCCALGeometry() : 
   m_numActiveBlocks( 0 ){
 
-        for( int row = 0; row < kCCALBlocksTall; row++ ){
-	  for( int col = 0; col < kCCALBlocksWide; col++ ){
-      
-	    // transform to beam axis
-	    m_positionOnFace[row][col] = 
-	      DVector2(  ( (double)col - kCCALMidBlock  + 0.5 ) * blockSize(),
-			 ( (double)row - kCCALMidBlock + 0.5 ) * blockSize() );
-	    
-	    m_activeBlock[row][col] = true;
-      
-	    // build the "channel map"
-	    m_channelNumber[row][col]    =  m_numActiveBlocks;
-	    m_row[m_numActiveBlocks]     =  row;
-	    m_column[m_numActiveBlocks]  =  col;
-	    
-	    m_numActiveBlocks++;
-	  }
+	for( int row = 0; row < kCCALBlocksTall; row++ ){
+		for( int col = 0; col < kCCALBlocksWide; col++ ){
+			
+			// transform to beam axis
+			m_positionOnFace[row][col] = 
+			DVector2(  ( (double)col - kCCALMidBlock  + 0.5 ) * blockSize(),
+				( (double)row - kCCALMidBlock + 0.5 ) * blockSize() );
+			
+			m_activeBlock[row][col] = true;
+			
+			// build the "channel map"
+			m_channelNumber[row][col]    =  m_numActiveBlocks;
+			m_row[m_numActiveBlocks]     =  row;
+			m_column[m_numActiveBlocks]  =  col;
+			
+			m_numActiveBlocks++;
+		}
 	}
 }
 
@@ -49,7 +49,7 @@ DCCALGeometry::isBlockActive( int row, int column ) const
 	// 12/13/05  DL
 	if( row < 0 ||  row >= kCCALBlocksTall )return false;
 	if( column < 0 ||  column >= kCCALBlocksWide )return false;
-
+	
 	assert(    row >= 0 &&    row < kCCALBlocksTall );
 	assert( column >= 0 && column < kCCALBlocksWide );
 

src/libraries/CCAL/DCCALGeometry_factory.h

@@ -37,13 +37,13 @@ class DCCALGeometry_factory:public jana::JFactory<DCCALGeometry>{
 		//------------------
 		// evnt
 		//------------------
-		 jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
-		 {
+		jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
+		{
 			// Reuse existing DBCALGeometry object.
 			if( ccalgeometry ) _data.push_back( ccalgeometry );
 
-			 return NOERROR;
-		 }
+			return NOERROR;
+		}
 
 		//------------------
 		// erun