From 075d336dc31a366825ef1dd4ec9595cdfa387c29 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 08:05:03 -0400 Subject: [PATCH 01/19] add routine to check for FCAL insert; update insert-related code to look for new names for insert blocks and rows --- src/libraries/HDGEOMETRY/DGeometry.cc | 27 ++++++++++++++++++++++----- src/libraries/HDGEOMETRY/DGeometry.h | 1 + 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/src/libraries/HDGEOMETRY/DGeometry.cc b/src/libraries/HDGEOMETRY/DGeometry.cc index 9a08bad3f..265ce8df5 100644 --- a/src/libraries/HDGEOMETRY/DGeometry.cc +++ b/src/libraries/HDGEOMETRY/DGeometry.cc @@ -2058,13 +2058,30 @@ bool DGeometry::GetCCALPosition(double &x,double &y,double &z) const } } +// Check for presence of FCAL2 insert +bool DGeometry::HaveInsert() const{ + int ncopy=0; + jgeom->SetVerbose(0); + bool have_insert + =Get("//composition[@name='XTrow0']/mposX[@volume='XTModule']/@ncopy", + ncopy); + jgeom->SetVerbose(1); + return have_insert; +} + //--------------------------------- // GetFCALInsertRowSize //--------------------------------- bool DGeometry::GetFCALInsertRowSize(int &insert_row_size) const { jgeom->SetVerbose(0); // don't print error messages for optional detector elements - bool good = Get("//composition[@name='LeadTungstateFullRow']/mposX[@volume='LTBLwrapped']/@ncopy",insert_row_size); + + bool good = Get("//composition[@name='XTrow0']/mposX[@volume='XTModule']/@ncopy",insert_row_size); + // For backward compatibility with prototype geometry definition: + if (!good){ + good = Get("//composition[@name='LeadTungstateFullRow']/mposX[@volume='LTBLwrapped']/@ncopy",insert_row_size); + } + jgeom->SetVerbose(1); // reenable error messages if(!good){ @@ -2083,7 +2100,7 @@ bool DGeometry::GetFCALInsertRowSize(int &insert_row_size) const bool DGeometry::GetFCALBlockSize(vector &block) const { jgeom->SetVerbose(0); // don't print error messages for optional detector elements - bool good = Get("//box[@name='LGBL']/@X_Y_Z",block); + bool good = Get("//box[@name='LGBU']/@X_Y_Z",block); jgeom->SetVerbose(1); // reenable error messages if(!good){ @@ -2101,12 +2118,12 @@ bool DGeometry::GetFCALBlockSize(vector &block) const bool DGeometry::GetFCALInsertBlockSize(vector &block) const { jgeom->SetVerbose(0); // don't print error messages for optional detector elements - bool good = Get("//box[@name='LTB1']/@X_Y_Z",block); + bool good = Get("//box[@name='XTMD']/@X_Y_Z",block); + // for backward compatiblity + if (!good) good=Get("//box[@name='LTB1']/@X_Y_Z",block); jgeom->SetVerbose(1); // reenable error messages if(!good){ - // NEED TO RETHINK ERROR REPORTING FOR OPTIONAL DETECTOR ELEMENTS - //_DBG_<<"Unable to retrieve ComptonEMcal position."< &block) const; bool GetFCALInsertBlockSize(vector &block) const; + bool HaveInsert() const; bool GetStartCounterGeom(vector >&pos, vector >&norm) const; // < vectors containing positions and norm 3-vectors for start counter From 178f9ee3a0d4d28199b2c298736a066234ca2af5 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 08:06:34 -0400 Subject: [PATCH 02/19] Remove DECALGeometry and move its functionality to the FCAL library where it is needed. --- src/libraries/ECAL/DECALGeometry.cc | 102 --------------------- src/libraries/ECAL/DECALGeometry.h | 73 --------------- src/libraries/ECAL/DECALGeometry_factory.h | 60 ------------ src/libraries/ECAL/DECALHit.h | 4 - src/libraries/ECAL/DECALHit_factory.cc | 95 ++++++++----------- src/libraries/ECAL/DECALHit_factory.h | 73 +++++++++------ src/libraries/ECAL/ECAL_init.cc | 2 - 7 files changed, 83 insertions(+), 326 deletions(-) delete mode 100644 src/libraries/ECAL/DECALGeometry.cc delete mode 100644 src/libraries/ECAL/DECALGeometry.h delete mode 100644 src/libraries/ECAL/DECALGeometry_factory.h diff --git a/src/libraries/ECAL/DECALGeometry.cc b/src/libraries/ECAL/DECALGeometry.cc deleted file mode 100644 index 1ece9034d..000000000 --- a/src/libraries/ECAL/DECALGeometry.cc +++ /dev/null @@ -1,102 +0,0 @@ -// $Id$ -// -/* - * File: DECALHit_factory.h - * - * Created on 01/16/2024 by A.S. - */ - - -#include -#include -using namespace std; - -#include "DECALGeometry.h" -#include "DVector2.h" - -//--------------------------------- -// DECALGeometry (Constructor) -//--------------------------------- -DECALGeometry::DECALGeometry() : - m_numActiveBlocks( 0 ){ - - for( int row = 0; row < kECALBlocksTall; row++ ){ - for( int col = 0; col < kECALBlocksWide; col++ ){ - - // transform to beam axis - m_positionOnFace[row][col] = - DVector2( ( (double)col - kECALMidBlock + 0.5 ) * blockSize(), - ( (double)row - kECALMidBlock + 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++; - } - } -} - -bool -DECALGeometry::isBlockActive( int row, int column ) const -{ - // I'm inserting these lines to effectively disable the - // two assert calls below. They are causing all programs - // (hd_dump, hdview) to exit, even when I'm not interested - // in the FCAL. This does not fix the underlying problem - // of why we're getting invalid row/column values. - // 12/13/05 DL - if( row < 0 || row >= kECALBlocksTall )return false; - if( column < 0 || column >= kECALBlocksWide )return false; - - assert( row >= 0 && row < kECALBlocksTall ); - assert( column >= 0 && column < kECALBlocksWide ); - - return m_activeBlock[row][column]; -} - -int -DECALGeometry::row( float y ) const -{ - return static_cast( y / blockSize() + kECALMidBlock ); -} - -int -DECALGeometry::column( float x ) const -{ - return static_cast( x / blockSize() + kECALMidBlock ); -} - -DVector2 -DECALGeometry::positionOnFace( int row, int column ) const -{ - assert( row >= 0 && row < kECALBlocksTall ); - assert( column >= 0 && column < kECALBlocksWide ); - - return m_positionOnFace[row][column]; -} - -DVector2 -DECALGeometry::positionOnFace( int channel ) const -{ - assert( channel >= 0 && channel < m_numActiveBlocks ); - - return positionOnFace( m_row[channel], m_column[channel] ); -} - -int -DECALGeometry::channel( int row, int column ) const -{ - if( isBlockActive( row, column ) ){ - - return m_channelNumber[row][column]; - } - else{ - - cerr << "ERROR: request for channel number of inactive block!" << endl; - return -1; - } -} diff --git a/src/libraries/ECAL/DECALGeometry.h b/src/libraries/ECAL/DECALGeometry.h deleted file mode 100644 index 8a82c9216..000000000 --- a/src/libraries/ECAL/DECALGeometry.h +++ /dev/null @@ -1,73 +0,0 @@ -// $Id$ -// -// File: DECALGeometry.h -// Created: Tue Nov 30 15:42:41 EST 2010 -// Creator: davidl (on Linux ifarml6 2.6.18-128.el5 x86_64) -// - -#ifndef _DECALGeometry_ -#define _DECALGeometry_ - -#include -#include - -#include "DVector2.h" -#include "units.h" - -class DECALGeometry:public jana::JObject{ - - public: - JOBJECT_PUBLIC(DECALGeometry); - - DECALGeometry(); - ~DECALGeometry(){} - - static const int kECALBlocksWide = 40; - static const int kECALBlocksTall = 40; - static const int kECALMaxChannels = kECALBlocksWide * kECALBlocksTall; - static const int kECALMidBlock = (kECALBlocksWide)/2; - static const int kECALBeamHoleSize = 2; - - static double blockSize() { return 2.09 * k_cm; } - static double blockLength(){ return 20.0 * k_cm; } - static double ecalFaceZ() { return 1279.376 * k_cm; } - - static double ecalMidplane() { return ecalFaceZ() + 0.5 * blockLength() ; } - - bool isBlockActive( int row, int column ) const; - int numActiveBlocks() const { return m_numActiveBlocks; } - - - DVector2 positionOnFace( int row, int column ) const; - DVector2 positionOnFace( int channel ) const; - - int channel( int row, int column ) const; - - int row ( int channel ) const { return m_row[channel]; } - int column( int channel ) const { return m_column[channel]; } - - // get row and column from x and y positions - int row ( float y ) const; - int column( float x ) const; - - void toStrings(vector > &items)const{ - AddString(items, "kECALBlocksWide", "%d", (int) kECALBlocksWide); - AddString(items, "kECALBlocksTall", "%d", (int) kECALBlocksTall); - AddString(items, "kECALMaxChannels", "%d",(int) kECALMaxChannels); - AddString(items, "kECALBeamHoleSize", "%2.3f",(int) kECALBeamHoleSize); - } - - private: - - bool m_activeBlock[kECALBlocksTall][kECALBlocksWide]; - DVector2 m_positionOnFace[kECALBlocksTall][kECALBlocksWide]; - - int m_channelNumber[kECALBlocksTall][kECALBlocksWide]; - int m_row[kECALMaxChannels]; - int m_column[kECALMaxChannels]; - - int m_numActiveBlocks; -}; - -#endif // _DECALGeometry_ - diff --git a/src/libraries/ECAL/DECALGeometry_factory.h b/src/libraries/ECAL/DECALGeometry_factory.h deleted file mode 100644 index 750fb957c..000000000 --- a/src/libraries/ECAL/DECALGeometry_factory.h +++ /dev/null @@ -1,60 +0,0 @@ -// $Id$ -/* - * File: DECALHit_factory.h - * - * Created on 01/16/2024 by A.S. - */ - -#ifndef _DECALGeometry_factory_ -#define _DECALGeometry_factory_ - -#include -#include "DECALGeometry.h" - -class DECALGeometry_factory:public jana::JFactory{ - public: - DECALGeometry_factory(){}; - ~DECALGeometry_factory(){}; - - DECALGeometry *ecalgeometry = nullptr; - - //------------------ - // brun - //------------------ - jerror_t brun(JEventLoop *loop, int32_t runnumber) - { - SetFactoryFlag(NOT_OBJECT_OWNER); - ClearFactoryFlag(WRITE_TO_OUTPUT); - - if( ecalgeometry ) delete ecalgeometry; - - ecalgeometry = new DECALGeometry(); - - return NOERROR; - } - - //------------------ - // evnt - //------------------ - jerror_t evnt(JEventLoop *loop, uint64_t eventnumber) - { - // Reuse existing DBCALGeometry object. - if( ecalgeometry ) _data.push_back( ecalgeometry ); - - return NOERROR; - } - - //------------------ - // erun - //------------------ - jerror_t erun(void) - { - if( ecalgeometry ) delete ecalgeometry; - ecalgeometry = NULL; - - return NOERROR; - } -}; - -#endif // _DECALGeometry_factory_ - diff --git a/src/libraries/ECAL/DECALHit.h b/src/libraries/ECAL/DECALHit.h index 972c86c0d..5f504d135 100644 --- a/src/libraries/ECAL/DECALHit.h +++ b/src/libraries/ECAL/DECALHit.h @@ -21,8 +21,6 @@ class DECALHit:public jana::JObject{ int row; int column; - float x; - float y; float E; float t; float intOverPeak; @@ -31,8 +29,6 @@ class DECALHit:public jana::JObject{ void toStrings(vector > &items)const{ AddString(items, "row", "%4d", row); AddString(items, "column", "%4d", column); - AddString(items, "x(cm)", "%3.1f", x); - AddString(items, "y(cm)", "%3.1f", y); AddString(items, "E(MeV)", "%2.3f", E); AddString(items, "t(ns)", "%2.3f", t); AddString(items, "integral over peak", "%2.3f", intOverPeak); diff --git a/src/libraries/ECAL/DECALHit_factory.cc b/src/libraries/ECAL/DECALHit_factory.cc index cab0ff8cd..2565ccda7 100644 --- a/src/libraries/ECAL/DECALHit_factory.cc +++ b/src/libraries/ECAL/DECALHit_factory.cc @@ -9,7 +9,6 @@ using namespace std; #include "ECAL/DECALDigiHit.h" -#include "ECAL/DECALGeometry.h" #include "ECAL/DECALHit_factory.h" #include "DAQ/Df250PulseIntegral.h" @@ -32,6 +31,26 @@ DECALHit_factory::DECALHit_factory(){ gPARMS->SetDefaultParameter("ECAL:HIT_DEBUG", HIT_DEBUG); gPARMS->SetDefaultParameter("ECAL:DB_PEDESTAL", DB_PEDESTAL); + m_numActiveBlocks=0; + for (int row=0;row=kECALMidBlock-1 && col<=kECALMidBlock + && row>=kECALMidBlock-1 && row<=kECALMidBlock){ + m_activeBlock[row][col]=false; + } + else{ + m_row[m_numActiveBlocks]=row; + m_column[m_numActiveBlocks]=col; + m_channelNumber[row][col]=m_numActiveBlocks; + m_activeBlock[row][col]=true; + + m_numActiveBlocks++; + } + } + + } + + } @@ -41,15 +60,15 @@ DECALHit_factory::DECALHit_factory(){ jerror_t DECALHit_factory::init(void) { // initialize calibration tables - vector< vector > gains_tmp(DECALGeometry::kECALBlocksTall, - vector(DECALGeometry::kECALBlocksWide)); - vector< vector > pedestals_tmp(DECALGeometry::kECALBlocksTall, - vector(DECALGeometry::kECALBlocksWide)); + vector< vector > gains_tmp(kECALBlocksTall, + vector(kECALBlocksWide)); + vector< vector > pedestals_tmp(kECALBlocksTall, + vector(kECALBlocksWide)); - vector< vector > time_offsets_tmp(DECALGeometry::kECALBlocksTall, - vector(DECALGeometry::kECALBlocksWide)); - vector< vector > adc_offsets_tmp(DECALGeometry::kECALBlocksTall, - vector(DECALGeometry::kECALBlocksWide)); + vector< vector > time_offsets_tmp(kECALBlocksTall, + vector(kECALBlocksWide)); + vector< vector > adc_offsets_tmp(kECALBlocksTall, + vector(kECALBlocksWide)); gains = gains_tmp; pedestals = pedestals_tmp; @@ -83,14 +102,6 @@ jerror_t DECALHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber) } pthread_mutex_unlock(&print_mutex); - // extract the ECAL Geometry - vector ecalGeomVect; - eventLoop->Get( ecalGeomVect ); - if (ecalGeomVect.size() < 1) - return OBJECT_NOT_AVAILABLE; - const DECALGeometry& ecalGeom = *(ecalGeomVect[0]); - - vector< double > ecal_gains_ch; vector< double > ecal_pedestals_ch; vector< double > time_offsets_ch; @@ -130,10 +141,10 @@ jerror_t DECALHit_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber) jout << "Error loading /ECAL/adc_offsets !" << endl; - LoadECALConst(gains, ecal_gains_ch, ecalGeom); - LoadECALConst(pedestals, ecal_pedestals_ch, ecalGeom); - LoadECALConst(time_offsets, time_offsets_ch, ecalGeom); - LoadECALConst(adc_offsets, adc_offsets_ch, ecalGeom); + LoadECALConst(gains, ecal_gains_ch); + LoadECALConst(pedestals, ecal_pedestals_ch); + LoadECALConst(time_offsets, time_offsets_ch); + LoadECALConst(adc_offsets, adc_offsets_ch); if(HIT_DEBUG == 1){ @@ -208,15 +219,6 @@ jerror_t DECALHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber) /// data in HDDM format. The HDDM event source will copy /// the precalibrated values directly into the _data vector. - - // extract the ECAL Geometry - vector ecalGeomVect; - eventLoop->Get( ecalGeomVect ); - if (ecalGeomVect.size() < 1) - return OBJECT_NOT_AVAILABLE; - const DECALGeometry& ecalGeom = *(ecalGeomVect[0]); - - vector digihits; loop->Get(digihits); @@ -269,12 +271,6 @@ jerror_t DECALHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber) hit->E = pulse_int_ped_subt; hit->t = pulse_time_correct; - // Get position of blocks on front face. (This should really come from - // hdgeant directly so the poisitions can be shifted in mcsmear.) - DVector2 pos = ecalGeom.positionOnFace(hit->row, hit->column); - hit->x = pos.X(); - hit->y = pos.Y(); - if(pulse_amplitude > 0){ hit->intOverPeak = (pulse_int - integratedPedestal)/pulse_amplitude; } else @@ -307,31 +303,18 @@ jerror_t DECALHit_factory::fini(void) } -void DECALHit_factory::LoadECALConst(ecal_constants_t &table, const vector &ecal_const_ch, - const DECALGeometry &ecalGeom){ - - char str[256]; - - // if (ecalGeom.numActiveBlocks() != ECAL_MAX_CHANNELS) { - // sprintf(str, "ECAL geometry is wrong size! channels=%d (should be %d)", - // ecalGeom.numActiveBlocks(), ECAL_MAX_CHANNELS); - // throw JException(str); - // } - - +void DECALHit_factory::LoadECALConst(ecal_constants_t &table, const vector &ecal_const_ch){ for (int ch = 0; ch < static_cast(ecal_const_ch.size()); ch++) { - - // make sure that we don't try to load info for channels that don't exist - if (ch == ecalGeom.numActiveBlocks()) - break; - - int row = ecalGeom.row(ch); - int col = ecalGeom.column(ch); + + int row = m_row[ch]; + int col = m_column[ch]; // results from DECALGeometry should be self consistent, but add in some // sanity checking just to be sure - if (ecalGeom.isBlockActive(row,col) == false) { + if (isBlockActive(row,col) == false) { + char str[256]; + sprintf(str, "DECALHit: Loading ECAL constant for inactive channel! " "row=%d, col=%d", row, col); throw JException(str); diff --git a/src/libraries/ECAL/DECALHit_factory.h b/src/libraries/ECAL/DECALHit_factory.h index ca6714504..63190550a 100644 --- a/src/libraries/ECAL/DECALHit_factory.h +++ b/src/libraries/ECAL/DECALHit_factory.h @@ -20,35 +20,50 @@ using namespace std; typedef vector< vector > ecal_constants_t; class DECALHit_factory:public jana::JFactory{ - public: - DECALHit_factory(); - ~DECALHit_factory(){}; - - ecal_constants_t gains; - ecal_constants_t pedestals; - ecal_constants_t time_offsets; - ecal_constants_t adc_offsets; - - double adc_en_scale; - double adc_time_scale; - - double base_time; - - private: - jerror_t init(void); ///< Called once at program start.2 - jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber); ///< Called everytime a new run number is detected. - jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber); ///< Called every event. - jerror_t erun(void); ///< Called everytime run number changes, provided brun has been called. - jerror_t fini(void); ///< Called after last event of last event source has been processed. - - void LoadECALConst( ecal_constants_t &table, - const vector &ecal_const_ch, - const DECALGeometry &ecalGeom); - - - unsigned int DB_PEDESTAL; - unsigned int HIT_DEBUG; - +public: + DECALHit_factory(); + ~DECALHit_factory(){}; + + ecal_constants_t gains; + ecal_constants_t pedestals; + ecal_constants_t time_offsets; + ecal_constants_t adc_offsets; + + double adc_en_scale; + double adc_time_scale; + + double base_time; + + bool isBlockActive( int row, int column ) const { + if( row < 0 || row >= kECALBlocksTall )return false; + if( column < 0 || column >= kECALBlocksWide )return false; + + return m_activeBlock[row][column]; + } + +private: + jerror_t init(void); ///< Called once at program start.2 + jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber); ///< Called everytime a new run number is detected. + jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber); ///< Called every event. + jerror_t erun(void); ///< Called everytime run number changes, provided brun has been called. + jerror_t fini(void); ///< Called after last event of last event source has been processed. + + void LoadECALConst(ecal_constants_t &table, + const vector &ecal_const_ch); + + static const int kECALBlocksWide = 40; + static const int kECALBlocksTall = 40; + static const int kECALMidBlock = 20; + static const int kECALMaxChannels = kECALBlocksWide * kECALBlocksTall; + + unsigned int DB_PEDESTAL; + unsigned int HIT_DEBUG; + + bool m_activeBlock[kECALBlocksTall][kECALBlocksWide]; + int m_channelNumber[kECALBlocksTall][kECALBlocksWide]; + int m_row[kECALMaxChannels]; + int m_column[kECALMaxChannels]; + int m_numActiveBlocks; }; #endif // _DECALHit_factory_ diff --git a/src/libraries/ECAL/ECAL_init.cc b/src/libraries/ECAL/ECAL_init.cc index 6bfbe5d8e..cb455f323 100644 --- a/src/libraries/ECAL/ECAL_init.cc +++ b/src/libraries/ECAL/ECAL_init.cc @@ -4,7 +4,6 @@ using namespace jana; #include -#include #include "DECALDigiHit.h" #include "DECALHit_factory.h" @@ -20,7 +19,6 @@ jerror_t ECAL_init(JEventLoop *loop) loop->AddFactory(new DECALHit_factory()); loop->AddFactory(new JFactory("TRUTH")); loop->AddFactory(new DECALTruthShower_factory()); - loop->AddFactory(new DECALGeometry_factory()); return NOERROR; } From c675410f3e2399704a8e024c45abb4f9f689e394 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 08:11:46 -0400 Subject: [PATCH 03/19] Add routine to get FCAL geometry (including insert) from the xml; the old method using a uniform grid is kept for backward compatibility. --- src/libraries/FCAL/DFCALGeometry.cc | 146 +++++++++++++++++++++++++--- src/libraries/FCAL/DFCALGeometry.h | 11 ++- 2 files changed, 141 insertions(+), 16 deletions(-) diff --git a/src/libraries/FCAL/DFCALGeometry.cc b/src/libraries/FCAL/DFCALGeometry.cc index 9a8767cc8..384b26228 100644 --- a/src/libraries/FCAL/DFCALGeometry.cc +++ b/src/libraries/FCAL/DFCALGeometry.cc @@ -16,24 +16,144 @@ using namespace std; // DFCALGeometry (Constructor) //--------------------------------- DFCALGeometry::DFCALGeometry(const DGeometry *geom){ + // Find position of upstream face of FCAL + geom->GetFCALPosition(m_FCALdX,m_FCALdY,m_FCALfront); + + // Find the size of the sensitive volume of each lead glass block + vectorblock; + geom->Get("//box[@name='LGBL']/@X_Y_Z",block); + m_sensitiveBlockSize=block[0]; + + // Check for presence of PbWO4 insert + if (geom->HaveInsert()){ // Geometry based on survey data for positions + // Get size of the insert + int insert_row_size=0; + geom->GetFCALInsertRowSize(insert_row_size); + m_insertSize=insertBlockSize()*double(insert_row_size/2); + + // Find the size of the sensitive volume of each PWO crystal + geom->Get("//box[@name='XTBL']/@X_Y_Z",block); + m_insertSensitiveBlockSize=block[0]; + + // Get the z-position of the upstream face of the insert + geom->GetECALZ(m_insertFront); + + // Get the x-y positions of the crystals + GetSurveyGeometry(geom); + } + else{ // Geometry based on grid(s) of fixed dimensions and block sizes + GetGridGeometry(geom); + } +} + +// The following routine looks up the position of each row in the xml +// specification for the geometry +void DFCALGeometry::GetSurveyGeometry(const DGeometry *geom){ + unsigned int ch=0; + + // Initialize active block map + for( int row = 0; row < 2*kBlocksTall; row++ ){ + for( int col = 0; col < 2*kBlocksWide; col++ ){ + m_activeBlock[row][col]=false; + } + } + + // Extract block positions for each section of the lead glass from xml + string section_names[4]={"LGLowerRow","LGUpperRow","LGNorthRow","LGSouthRow"}; + int num_rows=19; + for (int i=0;i<4;i++){ + if (i>1) num_rows=38; + for (int j=0;jGet(my_mpos_string+"@ncopy",ncopy); + geom->Get(my_mpos_string+"column/@value",col0); + geom->Get(my_mpos_string+"row/@value",row); + geom->Get(my_mpos_string+"@X0",x0); + geom->Get(my_mpos_string+"@dX",dx); + string my_pos_string="//posXYZ[@volume='"+my_row_string+"']/"; + vectorpos; + geom->Get(my_pos_string+"@X_Y_Z",pos); + //vectorrot; + //geom->Get(my_pos_string+"@rot",rot); + //double phi=rot[2]*M_PI/180.; + for (int col=col0;colGet(my_mpos_string+"@ncopy",ncopy); + geom->Get(my_mpos_string+"@X0",x0); + geom->Get(my_mpos_string+"@dX",dx); + string my_pos_string="//posXYZ[@volume='"+my_row_string+"']/"; + vectorpos; + geom->Get(my_pos_string+"@X_Y_Z",pos); + vectorrot; + geom->Get(my_pos_string+"@rot",rot); + double phi=rot[2]*M_PI/180.; + int col0=0,my_row=i; + if (i>=40){ //handle right side of beam hole + col0=21; + my_row=i-21; + } + for (int col=col0;colGetFCALInsertRowSize(insert_row_size); m_insertSize=insertBlockSize()*double(insert_row_size/2); - - geom->GetFCALPosition(m_FCALdX,m_FCALdY,m_FCALfront); + DVector2 XY0(m_FCALdX,m_FCALdY); - vectorblock; - geom->GetFCALBlockSize(block); - double back=m_FCALfront+block[2]; - geom->GetFCALInsertBlockSize(block); - m_insertFront=0.5*(back+m_FCALfront-block[2]); + // The following is for backward compatibility with an older model for the + // FCAL insert, now superceded by a more realistic model of the geometry + // (SJT 4/30/24) + double back=m_FCALfront+blockLength(); + m_insertFront=0.5*(back+m_FCALfront-insertBlockLength()); // Initilize the list of active blocks to false, to be adjusted for the // actual geometry below. @@ -73,6 +193,7 @@ DFCALGeometry::DFCALGeometry(const DGeometry *geom){ } } } + m_numFcalChannels=ch; if (insert_row_size>0){ m_insertMidBlock=(insert_row_size-1)/2; @@ -198,11 +319,8 @@ bool DFCALGeometry::isFiducial(double x,double y) const{ bool DFCALGeometry::inInsert(int channel) const{ - if (fabs(positionOnFace(channel).X()-m_FCALdX) Date: Thu, 22 Aug 2024 08:30:41 -0400 Subject: [PATCH 04/19] Get ECAL hits and and them to the userhit array --- src/libraries/FCAL/DFCALCluster_factory.cc | 34 ++++++++++++++++++++-- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/src/libraries/FCAL/DFCALCluster_factory.cc b/src/libraries/FCAL/DFCALCluster_factory.cc index 98a723cd1..c920bcf41 100644 --- a/src/libraries/FCAL/DFCALCluster_factory.cc +++ b/src/libraries/FCAL/DFCALCluster_factory.cc @@ -17,8 +17,9 @@ using namespace jana; #include "FCAL/DFCALCluster.h" #include "FCAL/DFCALCluster_factory.h" #include "FCAL/DFCALHit.h" +#include #include "FCAL/DFCALGeometry.h" -#include "HDGEOMETRY/DGeometry.h" +#include // Used to sort hits by Energy @@ -91,10 +92,13 @@ jerror_t DFCALCluster_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber) vector fcalhits; eventLoop->Get(fcalhits); + // Hits in PWO insert + vector ecalhits; + eventLoop->Get(ecalhits); // LED events will have hits in nearly every channel. Do NOT // try clusterizing if more than 250 hits in FCAL - if(fcalhits.size() > MAX_HITS_FOR_CLUSTERING) return NOERROR; + if(fcalhits.size()+ecalhits.size() > MAX_HITS_FOR_CLUSTERING) return NOERROR; const DFCALGeometry* fcalGeom=NULL; eventLoop->GetSingle(fcalGeom); @@ -125,8 +129,32 @@ jerror_t DFCALCluster_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber) << nhits << " larger than " << FCAL_USER_HITS_MAX << endl; break; } - } + // Now add the ECAL hits + if (nhits < (int) FCAL_USER_HITS_MAX){ + for (vector::const_iterator hit = ecalhits.begin(); + hit != ecalhits.end(); hit++ ) { + if ( (**hit).E < 1e-6 ) continue; + + hits->hit[nhits].id = (**hit).id; + hits->hit[nhits].E = (**hit).E; + hits->hit[nhits].t = (**hit).t; + hits->hit[nhits].intOverPeak=(**hit).intOverPeak; + + int ch=fcalGeom->channel(100+(**hit).row,100+(**hit).column); + DVector2 positionOnFace=fcalGeom->positionOnFace(ch); + hits->hit[nhits].x = positionOnFace.X(); + hits->hit[nhits].y = positionOnFace.Y(); + + nhits++; + + if (nhits >= (int) FCAL_USER_HITS_MAX) { + cout << "ERROR: FCALCluster_factory: number of hits " + << nhits << " larger than " << FCAL_USER_HITS_MAX << endl; + break; + } + } + } hits->nhits = nhits; int hitUsed[nhits]; From 6b85087fe0d3a19cd80c27754f61671a32aa8c07 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 17:16:00 -0400 Subject: [PATCH 05/19] Add routine to set the color for ECAL/FCAL/CCAL showers by log(E). --- .../Analysis/hdview2/hdv_mainframe.cc | 47 ++++++++++++++++++- src/programs/Analysis/hdview2/hdv_mainframe.h | 3 +- 2 files changed, 47 insertions(+), 3 deletions(-) diff --git a/src/programs/Analysis/hdview2/hdv_mainframe.cc b/src/programs/Analysis/hdview2/hdv_mainframe.cc index 276c92d98..4c32d1df9 100644 --- a/src/programs/Analysis/hdview2/hdv_mainframe.cc +++ b/src/programs/Analysis/hdview2/hdv_mainframe.cc @@ -1712,7 +1712,7 @@ void hdv_mainframe::DrawDetectorsXY(void) fcalblocks.clear(); if(GetCheckButton("fcal")){ - for(unsigned int chan=0; channumChannels(); chan++){ + for(int chan=0; channumChannels(); chan++){ if (fcalgeom->isBlockActive(chan)==false) continue; DVector2 fcalBlockPos=fcalgeom->positionOnFace(chan); @@ -2092,7 +2092,7 @@ void hdv_mainframe::DrawDetectorsRPhi(void) shift[7].Set(+blocksize/2, -blocksize/2); // define a single enclosed space } fcalblocks.clear(); - for(unsigned int chan=0; channumChannels(); chan++){ + for(int chan=0; channumChannels(); chan++){ if (fcalgeom->isBlockActive(chan)==false) continue; DVector2 fcalBlockPos=fcalgeom->positionOnFace(chan); @@ -2686,3 +2686,46 @@ void hdv_mainframe::RedrawAuxillaryWindows(void) if( fmwpcmf ) fmwpcmf->DoNewEvent(); } + +// Set the color of a hit FCAL/ECAL/CCAL block according to energy. +// The aim is to have a log scale in energy (see BCAL) +void hdv_mainframe::SetCalorimeterEnergyColor(TPolyLine *poly,double E) const{ + E *= 1000; // Change Energy to MeV + double logE = log10(E); + + float r,g,b; + if (logE<0){ + r = 1.; + g = 1.; + b = 1.; + } else { + if (logE<1){ + r = 1.; + g = 1.; + b = 1.-logE; + } else { + if (logE<2){ + r = 1.; + g = 1.-(logE-1); + b = 0.; + } else { + if (logE<3){ + r = 1.; + g = 0.; + b = 1.-(logE-2); + } else { + if (logE<4){ + r = 1.-(logE-3); + g = 0.; + b = 1.; + } else { + r = 0; + g = 0; + b = 0; + } + } + } + } + } + poly->SetFillColor(TColor::GetColor(r,g,b)); +} diff --git a/src/programs/Analysis/hdview2/hdv_mainframe.h b/src/programs/Analysis/hdview2/hdv_mainframe.h index f53ade4e6..09f3f59cd 100644 --- a/src/programs/Analysis/hdview2/hdv_mainframe.h +++ b/src/programs/Analysis/hdview2/hdv_mainframe.h @@ -151,7 +151,8 @@ class hdv_mainframe:public TGMainFrame { TPolyLine* GetCCALPolyLine(int row, int col); TPolyLine* GetBCALPolyLine(int mod, int layer, int sector); TPolyLine* GetTOFPolyLine(int translate_side, int tof_ch); - + void SetCalorimeterEnergyColor(TPolyLine *poly,double E) const; + void AddGraphicsSideA(std::vector &v); void AddGraphicsSideB(std::vector &v); void AddGraphicsEndA(std::vector &v); From 46b31106c094701fb0de3b4867277240be8b8321 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 17:16:51 -0400 Subject: [PATCH 06/19] Get ECAL hits; use new method to set color for ECAL/FCAL/CCAL hits. --- src/programs/Analysis/hdview2/MyProcessor.cc | 111 ++++--------------- src/programs/Analysis/hdview2/MyProcessor.h | 2 +- 2 files changed, 20 insertions(+), 93 deletions(-) diff --git a/src/programs/Analysis/hdview2/MyProcessor.cc b/src/programs/Analysis/hdview2/MyProcessor.cc index 919081313..92d30e2a3 100644 --- a/src/programs/Analysis/hdview2/MyProcessor.cc +++ b/src/programs/Analysis/hdview2/MyProcessor.cc @@ -38,6 +38,7 @@ using namespace std; #include "JANA/JGeometry.h" #include "TRACKING/DMCTrajectoryPoint.h" #include "FCAL/DFCALHit.h" +#include "ECAL/DECALHit.h" #include "CCAL/DCCALHit.h" #include "CCAL/DCCALShower.h" #include "TOF/DTOFGeometry.h" @@ -679,6 +680,20 @@ void MyProcessor::FillGraphics(void) // FCAL hits if(hdvmf->GetCheckButton("fcal")){ + // Get FCAL2 insert hits if present + vector ecalhits; + loop->Get(ecalhits); + + for(unsigned int i=0; iE<0.) continue; + + TPolyLine *poly = hdvmf->GetFCALPolyLine(100+hit->row,100+hit->column); + if(!poly) continue; + + hdvmf->SetCalorimeterEnergyColor(poly,hit->E); + } + vector fcalhits; loop->Get(fcalhits); @@ -686,58 +701,8 @@ void MyProcessor::FillGraphics(void) const DFCALHit *hit = fcalhits[i]; TPolyLine *poly = hdvmf->GetFCALPolyLine(hit->row, hit->column); if(!poly)continue; - -#if 0 - double a = hit->E/0.005; - double f = sqrt(a>1.0 ? 1.0:a<0.0 ? 0.0:a); - double grey = 0.8; - double s = 1.0 - f; - - float r = s*grey; - float g = s*grey; - float b = f*(1.0-grey) + grey; -#endif - - // The aim is to have a log scale in energy (see BCAL) - double E = 1000*hit->E; // Change Energy to MeV - if(E<0.0) continue; - double logE = log10(E); - float r,g,b; - if (logE<0){ - r = 1.; - g = 1.; - b = 1.; - } else { - if (logE<1){ - r = 1.; - g = 1.; - b = 1.-logE; - } else { - if (logE<2){ - r = 1.; - g = 1.-(logE-1); - b = 0.; - } else { - if (logE<3){ - r = 1.; - g = 0.; - b = 1.-(logE-2); - } else { - if (logE<4){ - r = 1.-(logE-3); - g = 0.; - b = 1.; - } else { - r = 0; - g = 0; - b = 0; - } - } - } - } - } - poly->SetFillColor(TColor::GetColor(r,g,b)); + hdvmf->SetCalorimeterEnergyColor(poly,hit->E); } } @@ -751,47 +716,9 @@ void MyProcessor::FillGraphics(void) TPolyLine *poly = hdvmf->GetCCALPolyLine(hit->row, hit->column); if(!poly)continue; - // The aim is to have a log scale in energy (see BCAL) - double E = 1000*hit->E; // Change Energy to MeV - E /= 1000.0; // CCAL is currently not calibrated and appears to be at least 1E3 too large 2018-12-10 DL - if(E<0.0) continue; - double logE = log10(E); - - float r,g,b; - if (logE<0){ - r = 1.; - g = 1.; - b = 1.; - } else { - if (logE<1){ - r = 1.; - g = 1.; - b = 1.-logE; - } else { - if (logE<2){ - r = 1.; - g = 1.-(logE-1); - b = 0.; - } else { - if (logE<3){ - r = 1.; - g = 0.; - b = 1.-(logE-2); - } else { - if (logE<4){ - r = 1.-(logE-3); - g = 0.; - b = 1.; - } else { - r = 0; - g = 0; - b = 0; - } - } - } - } - } - poly->SetFillColor(TColor::GetColor(r,g,b)); + // The following code converts from GeV units to MeV units, but the + // CCAL energy units are already in MeV... + hdvmf->SetCalorimeterEnergyColor(poly,0.001*hit->E); } } diff --git a/src/programs/Analysis/hdview2/MyProcessor.h b/src/programs/Analysis/hdview2/MyProcessor.h index d62565414..be72fad63 100644 --- a/src/programs/Analysis/hdview2/MyProcessor.h +++ b/src/programs/Analysis/hdview2/MyProcessor.h @@ -112,7 +112,7 @@ class MyProcessor:public JEventProcessor private: hdv_mainframe *hdvmf; - hdv_fulllistframe *fulllistmf; + hdv_fulllistframe *fulllistmf=NULL; hdv_debugerframe *debugermf; JEventLoop *loop; JEvent last_jevent; From cbad53d2e5b67e2e0dd62897c77ec5c91fbb4671 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 17:19:27 -0400 Subject: [PATCH 07/19] Remove reference to DECALGeometry --- src/libraries/HDDM/DEventSourceHDDM.cc | 31 -------------------------- 1 file changed, 31 deletions(-) diff --git a/src/libraries/HDDM/DEventSourceHDDM.cc b/src/libraries/HDDM/DEventSourceHDDM.cc index 1bcb11862..bcb492ae5 100644 --- a/src/libraries/HDDM/DEventSourceHDDM.cc +++ b/src/libraries/HDDM/DEventSourceHDDM.cc @@ -40,7 +40,6 @@ using namespace std; #include #include #include -#include #include @@ -1848,15 +1847,7 @@ jerror_t DEventSourceHDDM::Extract_DECALHit(hddm_s::HDDM *record, if (tag != "" && tag != "TRUTH") return OBJECT_NOT_AVAILABLE; - // extract the ECAL Geometry (for isBlockActive() and positionOnFace()) - vector ecalGeomVect; - eventLoop->Get( ecalGeomVect ); - if (ecalGeomVect.size() < 1) - return OBJECT_NOT_AVAILABLE; - const DECALGeometry& ecalGeom = *(ecalGeomVect[0]); - vector data; - int hitId = 0; if (tag == "") { @@ -1865,23 +1856,12 @@ jerror_t DEventSourceHDDM::Extract_DECALHit(hddm_s::HDDM *record, for (iter = hits.begin(); iter != hits.end(); ++iter) { int row = iter->getRow(); int column = iter->getColumn(); - - // Filter out non-physical blocks here - if (!ecalGeom.isBlockActive(row, column)) - continue; - - // Get position of blocks on front face. (This should really come from - // hdgeant directly so the poisitions can be shifted in mcsmear.) - DVector2 pos = ecalGeom.positionOnFace(row, column); DECALHit *mchit = new DECALHit(); mchit->row = row; mchit->column = column; - mchit->x = pos.X(); - mchit->y = pos.Y(); mchit->E = iter->getE(); mchit->t = iter->getT(); - mchit->id = hitId++; mchit->intOverPeak = 5.; data.push_back(mchit); } @@ -1893,23 +1873,12 @@ jerror_t DEventSourceHDDM::Extract_DECALHit(hddm_s::HDDM *record, for (iter = hits.begin(); iter != hits.end(); ++iter) { int row = iter->getRow(); int column = iter->getColumn(); - - // Filter out non-physical blocks here - if (!ecalGeom.isBlockActive(row, column)) - continue; - - // Get position of blocks on front face. (This should really come from - // hdgeant directly so the poisitions can be shifted in mcsmear.) - DVector2 pos = ecalGeom.positionOnFace(row, column); DECALHit *mchit = new DECALHit(); mchit->row = row; mchit->column = column; - mchit->x = pos.X(); - mchit->y = pos.Y(); mchit->E = iter->getE(); mchit->t = iter->getT(); - mchit->id = hitId++; mchit->intOverPeak = 5.; data.push_back(mchit); } From c72bb7417d322046591d9d8c6230a6d1c8325482 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 17:23:50 -0400 Subject: [PATCH 08/19] Convert DFCALClusterHit_t structure into a class and add addHit method for adding to the list of DFCALClusterHit_t objects; the vector is now called fHitLists to be more consistent with the style of the other private members. --- src/libraries/FCAL/DFCALCluster.h | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/src/libraries/FCAL/DFCALCluster.h b/src/libraries/FCAL/DFCALCluster.h index 9fca4eb25..c01c324a0 100644 --- a/src/libraries/FCAL/DFCALCluster.h +++ b/src/libraries/FCAL/DFCALCluster.h @@ -47,15 +47,17 @@ class DFCALCluster : public JObject { userhit_t hit[1]; } userhits_t; - typedef struct { - oid_t id; - int ch; - float x; - float y; - float E; - float t; - float intOverPeak; - } DFCALClusterHit_t; + class DFCALClusterHit_t{ + public: + DFCALClusterHit_t(oid_t id,int ch,double E,double x,double y,double t) + :id(id),ch(ch),E(E),x(x),y(y),t(t){} + oid_t id; + int ch; + double E; + double x; + double y; + double t; + }; void saveHits( const userhits_t* const hit ); @@ -88,10 +90,11 @@ class DFCALCluster : public JObject { void resetClusterHits(); bool update( const userhits_t* const hitList, double fcalFaceZ, const DFCALGeometry *fcalgeom ); + void addHit(oid_t id,int ch,double E,double x,double y,double t); // get hits that form a cluster after clustering is finished - inline const vector GetHits() const { return my_hits; } - inline uint32_t GetNHits(void) const { return my_hits.size(); } + inline const vector GetHits() const { return fHitList; } + inline uint32_t GetNHits(void) const { return fHitList.size(); } void toStrings(vector > &items) const { AddString(items, "x(cm)", "%3.1f", getCentroid().x()); @@ -141,7 +144,7 @@ class DFCALCluster : public JObject { double *fEallowed; // allowed energy of hit by cluster (GeV) // container for hits that form a cluster to be used after clustering is done - vector my_hits; + vector fHitList; }; From 3da9771b8249cfcbcb2691196ba5c90676eaa4f5 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 17:29:09 -0400 Subject: [PATCH 09/19] Make sure fHitList vector is cleared before adding hits info. Implemnent and use new addHit methods for this vector. Small change to check for whether or not a hit is in the insert. --- src/libraries/FCAL/DFCALCluster.cc | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/libraries/FCAL/DFCALCluster.cc b/src/libraries/FCAL/DFCALCluster.cc index e4d58318e..fd773b559 100644 --- a/src/libraries/FCAL/DFCALCluster.cc +++ b/src/libraries/FCAL/DFCALCluster.cc @@ -39,6 +39,7 @@ DFCALCluster::DFCALCluster( const int nhits ) fEallowed = 0; fEexpected = 0; } + fHitList.clear(); } DFCALCluster::~DFCALCluster() @@ -51,6 +52,7 @@ DFCALCluster::~DFCALCluster() delete [] fEallowed; if (fEexpected) delete [] fEexpected; + fHitList.clear(); } @@ -58,17 +60,10 @@ void DFCALCluster::saveHits( const userhits_t* const hits ) { for ( int i=0; i < fNhits; i++) { - DFCALClusterHit_t h; JObject::oid_t id = getHitID( hits, i ) ; if ( id != 0 ) { - h.id = (JObject::oid_t) id; - h.ch = getHitCh( hits, i ); - h.E = getHitE( hits, i ) ; - h.x = getHitX( hits, i ) ; - h.y = getHitY( hits, i ) ; - h.t = getHitT( hits, i ) ; - h.intOverPeak = getHitIntOverPeak( hits, i ); - my_hits.push_back(h); + addHit(id,getHitCh( hits, i ),getHitE( hits, i ),getHitX( hits, i ), + getHitY( hits, i ),getHitT( hits, i )); } else { static uint32_t Nwarns=0; @@ -82,6 +77,9 @@ void DFCALCluster::saveHits( const userhits_t* const hits ) } } +void DFCALCluster::addHit(oid_t id,int ch,double E,double x,double y,double t){ + fHitList.push_back(DFCALClusterHit_t(id,ch,E,x,y,t)); +} int DFCALCluster::addHit(const int ihit, const double frac) { @@ -276,7 +274,7 @@ void DFCALCluster::shower_profile( const userhits_t* const hitList, double y = hitList->hit[ihit].y; double moliere_radius=MOLIERE_RADIUS; double min_dist=fcalgeom->blockSize(); - if (fabs(x)insertSize() && fabs(y)insertSize()){ + if (hitList->hit[ihit].ch>=fcalgeom->numFcalChannels()){ moliere_radius=PbWO4_MOLIERE_RADIUS; min_dist=fcalgeom->insertBlockSize(); } From 9b6406cdbca3930c9722f47e4cf77cb957623dde Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 17:32:06 -0400 Subject: [PATCH 10/19] add channel for insert hits when setting userhit data. --- src/libraries/FCAL/DFCALCluster_factory.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/src/libraries/FCAL/DFCALCluster_factory.cc b/src/libraries/FCAL/DFCALCluster_factory.cc index c920bcf41..9c588fdef 100644 --- a/src/libraries/FCAL/DFCALCluster_factory.cc +++ b/src/libraries/FCAL/DFCALCluster_factory.cc @@ -142,6 +142,7 @@ jerror_t DFCALCluster_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber) hits->hit[nhits].intOverPeak=(**hit).intOverPeak; int ch=fcalGeom->channel(100+(**hit).row,100+(**hit).column); + hits->hit[nhits].ch=ch; DVector2 positionOnFace=fcalGeom->positionOnFace(ch); hits->hit[nhits].x = positionOnFace.X(); hits->hit[nhits].y = positionOnFace.Y(); From 9bf3a7c9ca80e6a7092e243762b66707bbaa13ea Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 17:34:16 -0400 Subject: [PATCH 11/19] Add Hitinfo class to store data from both ECAL and FCAL hits. --- .../FCAL/DFCALCluster_factory_Island.h | 27 ++++++++++++++----- 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/libraries/FCAL/DFCALCluster_factory_Island.h b/src/libraries/FCAL/DFCALCluster_factory_Island.h index 13451359c..957707c37 100644 --- a/src/libraries/FCAL/DFCALCluster_factory_Island.h +++ b/src/libraries/FCAL/DFCALCluster_factory_Island.h @@ -12,6 +12,8 @@ #include "DFCALCluster.h" #include "DFCALGeometry.h" #include "DFCALHit.h" +#include + #include "TMatrixD.h" #include "TH1D.h" #include "TH2D.h" @@ -21,6 +23,19 @@ class DFCALCluster_factory_Island:public jana::JFactory{ DFCALCluster_factory_Island(){}; ~DFCALCluster_factory_Island(){}; const char* Tag(void){return "Island";} + + class HitInfo{ + public: + HitInfo(JObject::oid_t id,int row,int column,double E,double x,double y,double t) + :id(id),row(row),column(column),E(E),x(x),y(y),t(t){} + JObject::oid_t id; + int row; + int column; + double E; + double x; + double y; + double t; + }; class PeakInfo{ public: @@ -40,15 +55,15 @@ class DFCALCluster_factory_Island:public jana::JFactory{ jerror_t erun(void); ///< Called everytime run number changes, provided brun has been called. jerror_t fini(void); ///< Called after last event of last event source has been processed. - void FindClusterCandidates(vector&fcal_hits, - vector>&clusterCandidates) const; - bool FitPeaks(const TMatrixD &W,double b,vector&hitList, + void FindClusterCandidates(vector&fcal_hits, + vector>&clusterCandidates) const; + bool FitPeaks(const TMatrixD &W,double b,vector&hitList, vector&peaks,PeakInfo &myPeak,double &chisq, unsigned int &ndf) const; - double CalcClusterEDeriv(double b,const DFCALHit *hit,const PeakInfo &myPeakInfo) const; - double CalcClusterXYDeriv(bool isXDeriv,double b,const DFCALHit *hit, + double CalcClusterEDeriv(double b,const HitInfo &hit,const PeakInfo &myPeakInfo) const; + double CalcClusterXYDeriv(bool isXDeriv,double b,const HitInfo &hit, const PeakInfo &myPeakInfo) const; - void SplitPeaks(const TMatrixD &W,double b,vector&hits, + void SplitPeaks(const TMatrixD &W,double b,vector&hits, vector&peaks,double &chisq,unsigned int &ndf) const; double TIME_CUT,MIN_CLUSTER_SEED_ENERGY,SHOWER_ENERGY_THRESHOLD; From 9cd44c376301b258c1f0a41d8c63626f6af0031a Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 17:37:25 -0400 Subject: [PATCH 12/19] Get ECAL hits and fill local HitInfo vector for both ECAL and FCAL hits. Use addHit interface to the fHitList vector in the cluster class instead of adding associated objects. --- .../FCAL/DFCALCluster_factory_Island.cc | 209 ++++++++++-------- 1 file changed, 120 insertions(+), 89 deletions(-) diff --git a/src/libraries/FCAL/DFCALCluster_factory_Island.cc b/src/libraries/FCAL/DFCALCluster_factory_Island.cc index 555d88cf4..b4643f7b6 100644 --- a/src/libraries/FCAL/DFCALCluster_factory_Island.cc +++ b/src/libraries/FCAL/DFCALCluster_factory_Island.cc @@ -13,8 +13,9 @@ using namespace std; #include "DFCALCluster_factory_Island.h" using namespace jana; -inline bool FCALHit_E_cmp(const DFCALHit *a,const DFCALHit *b){ - return (a->E>b->E); +inline bool FCALHit_E_cmp(const DFCALCluster_factory_Island::HitInfo a, + const DFCALCluster_factory_Island::HitInfo b){ + return (a.E>b.E); } inline bool peak_E_cmp(DFCALCluster_factory_Island::PeakInfo a, @@ -122,19 +123,39 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe { vectorfcal_hits; loop->Get(fcal_hits); - if (fcal_hits.size()==0) return OBJECT_NOT_AVAILABLE; + vectorecal_hits; + loop->Get(ecal_hits); + + unsigned int total_hits=fcal_hits.size()+ecal_hits.size(); + if (total_hits==0) return OBJECT_NOT_AVAILABLE; // LED events will have hits in nearly every channel. Do NOT // try clusterizing if more than 250 hits in FCAL - if(fcal_hits.size() > MAX_HITS_FOR_CLUSTERING) return VALUE_OUT_OF_RANGE; - + if (total_hits > MAX_HITS_FOR_CLUSTERING) return VALUE_OUT_OF_RANGE; + + // Put hit information into local array + vectorhits; + for (size_t i=0;iid,myhit->row,myhit->column,myhit->E,myhit->x, + myhit->y,myhit->t)); + } + for (size_t i=0;irow; + int col=100+myhit->column; + DVector2 pos=dFCALGeom->positionOnFace(row,col); + hits.push_back(HitInfo(myhit->id,row,col,myhit->E,pos.X(),pos.Y(), + myhit->t)); + } + // Sort the hits according to energy. - stable_sort(fcal_hits.begin(),fcal_hits.end(),FCALHit_E_cmp); + stable_sort(hits.begin(),hits.end(),FCALHit_E_cmp); // Associate groups of adjacent hits into cluster candidates - vector>clusterCandidates; - FindClusterCandidates(fcal_hits,clusterCandidates); - + vector>clusterCandidates; + FindClusterCandidates(hits,clusterCandidates); + for (unsigned int i=0;iclusterCandidate=clusterCandidates[i]; + vectorclusterCandidate=clusterCandidates[i]; double Etot=0.,t=0,x=0,y=0; for (unsigned int k=0;kE; + double E=clusterCandidate[k].E; Etot+=E; - t+=E*clusterCandidate[k]->t; - x+=E*clusterCandidate[k]->x; - y+=E*clusterCandidate[k]->y; - - myCluster->AddAssociatedObject(clusterCandidate[k]); + t+=E*clusterCandidate[k].t; + x+=E*clusterCandidate[k].x; + y+=E*clusterCandidate[k].y; + + int channel=dFCALGeom->channel(clusterCandidate[k].row, + clusterCandidate[k].column); + myCluster->addHit(clusterCandidate[k].id,channel,E, + clusterCandidate[k].x,clusterCandidate[k].y, + clusterCandidate[k].t); } x/=Etot; y/=Etot; @@ -162,8 +187,8 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe myCluster->setEnergy(Etot); myCluster->setTimeEWeight(t); - int channel=dFCALGeom->channel(clusterCandidate[0]->row, - clusterCandidate[0]->column); + int channel=dFCALGeom->channel(clusterCandidate[0].row, + clusterCandidate[0].column); myCluster->setChannelEmax(channel); if (APPLY_S_CURVE_CORRECTION){ @@ -194,7 +219,7 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe //------------------------------------------------------------------------ //cout <<"Finding peaks for cluster " << i << endl; - vectorclusterHits=clusterCandidates[i]; + vectorclusterHits=clusterCandidates[i]; vectorpeaks; // Create weight matrix for hits @@ -203,10 +228,10 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe double Esum=0.; double R=0.; // simple average radius of cluster for (unsigned int j=0;jE; + const HitInfo hit=clusterHits[j]; + double E=hit.E; double varE=0.; - if (dFCALGeom->isInsertBlock(hit->row,hit->column)){ + if (dFCALGeom->isInsertBlock(hit.row,hit.column)){ varE=m_insert_Eres[0]+m_insert_Eres[1]*E+m_insert_Eres[2]*E*E; } else{ @@ -214,7 +239,7 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe } W(j,j)=1./varE; Esum+=E; - R+=sqrt(hit->x*hit->x+hit->y*hit->y); + R+=sqrt(hit.x*hit.x+hit.y*hit.y); } R/=double(num_hits); @@ -223,10 +248,10 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe // Find the minimum and maximum row and column numbers int min_row=1000,min_col=1000,max_row=0,max_col=0; for (unsigned int j=0;jrowrow; - if (clusterHits[j]->columncolumn; - if (clusterHits[j]->row>max_row) max_row=clusterHits[j]->row; - if (clusterHits[j]->column>max_col) max_col=clusterHits[j]->column; + if (clusterHits[j].rowmax_row) max_row=clusterHits[j].row; + if (clusterHits[j].column>max_col) max_col=clusterHits[j].column; } // Create arrays to represent the cluster of hits to aid in peak search int num_rows=max_row-min_row+3; @@ -237,20 +262,20 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe unsigned int jmax=0; double Emax=0.; for (unsigned int j=0;jcolumn-min_col+1; - int ir=clusterHits[j]->row-min_row+1; - Emap[ic][ir]=clusterHits[j]->E; + int ic=clusterHits[j].column-min_col+1; + int ir=clusterHits[j].row-min_row+1; + Emap[ic][ir]=clusterHits[j].E; imap[ic][ir]=j; - if (clusterHits[j]->E>Emax){ - Emax=clusterHits[j]->E; + if (clusterHits[j].E>Emax){ + Emax=clusterHits[j].E; jmax=j; } } // Compute estimate for shower shape parameter b double b=0.; - if (dFCALGeom->isInsertBlock(clusterHits[jmax]->row, - clusterHits[jmax]->column)){ + if (dFCALGeom->isInsertBlock(clusterHits[jmax].row, + clusterHits[jmax].column)){ b=INSERT_SHOWER_WIDTH_PAR0+INSERT_SHOWER_WIDTH_PAR1*R; } else{ @@ -262,11 +287,11 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe // Handle the interface between the insert and the lead glass blocks double Esum=0.; for (unsigned int j=0;jE; + double Ej=clusterHits[j].E; Esum+=Ej; } if (Emax>MIN_CLUSTER_SEED_ENERGY){ - PeakInfo myPeak(Esum,clusterHits[jmax]->x,clusterHits[jmax]->y,0,0, + PeakInfo myPeak(Esum,clusterHits[jmax].x,clusterHits[jmax].y,0,0, num_hits); bool good_fit=FitPeaks(W,b,clusterHits,peaks,myPeak,chisq,ndf); if (good_fit){ @@ -279,9 +304,9 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe myPeak.x=0.; myPeak.y=0.; for (unsigned int j=0;jE; - myPeak.x+=Ej*clusterHits[j]->x; - myPeak.y+=Ej*clusterHits[j]->y; + double Ej=clusterHits[j].E; + myPeak.x+=Ej*clusterHits[j].x; + myPeak.y+=Ej*clusterHits[j].y; } myPeak.x/=Esum; myPeak.y/=Esum; @@ -292,7 +317,7 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe chisq=0.; ndf=num_hits-3; for (unsigned int j=0;jE + double dE=clusterHits[j].E -Esum*CalcClusterEDeriv(b,clusterHits[j],myPeak); chisq+=W(j,j)*dE*dE; } @@ -333,8 +358,8 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe } } if (got_peak){ - double x=clusterHits[imap[ic][ir]]->x; - double y=clusterHits[imap[ic][ir]]->y; + double x=clusterHits[imap[ic][ir]].x; + double y=clusterHits[imap[ic][ir]].y; peak_candidates.push_back(PeakInfo(Esum,x,y,ic,ir,nhits_in_peak)); } }// cut on minimum energy of central block @@ -399,7 +424,7 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe chisq=0.; ndf=num_hits-3; for (unsigned int j=0;jE + double dE=clusterHits[j].E -Esum*CalcClusterEDeriv(b,clusterHits[j],peak_guess); chisq+=W(j,j)*dE*dE; } @@ -419,7 +444,7 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe PeakInfo myPeakInfo=peaks[0]; for (unsigned int k=0;kE; + double Ecalc=0.,E=clusterHits[k].E; Ecalc+=myPeakInfo.E*CalcClusterEDeriv(b,clusterHits[k],myPeakInfo); HistdE->Fill(E,E-Ecalc); } @@ -431,7 +456,7 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe // hit to see if we have excess energy that has not been accounted for vectorElist(clusterHits.size()); for (unsigned int m=0;mE; + Elist[m]=clusterHits[m].E; for (unsigned int k=0;kMIN_EXCESS_SEED_ENERGY){ // Make a peak candidate out of the excess energy in the cluster of hits - int ic=clusterHits[mmax]->column-min_col+1; - int ir=clusterHits[mmax]->row-min_row+1; + int ic=clusterHits[mmax].column-min_col+1; + int ir=clusterHits[mmax].row-min_row+1; int lo_col=ic-1; int hi_col=ic+1; int lo_row=ir-1; @@ -464,8 +489,8 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe } } } - double x=clusterHits[mmax]->x; - double y=clusterHits[mmax]->y; + double x=clusterHits[mmax].x; + double y=clusterHits[mmax].y; PeakInfo myPeak(excessE,x,y,ic,ir,num_peak_hits); // Save the current list of peaks @@ -520,7 +545,7 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe double fsum=0.,fmax=0.,t=0.; for (unsigned int j=0;jt; + t+=f*clusterHits[j].t; fsum+=f; if (f>fmax){ fmax=f; @@ -529,8 +554,8 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe } myCluster->setTimeEWeight(t/fsum); - int channel=dFCALGeom->channel(clusterHits[jmax]->row, - clusterHits[jmax]->column); + int channel=dFCALGeom->channel(clusterHits[jmax].row, + clusterHits[jmax].column); myCluster->setChannelEmax(channel); double xc=peaks[k].x,yc=peaks[k].y; @@ -551,20 +576,26 @@ jerror_t DFCALCluster_factory_Island::evnt(JEventLoop *loop, uint64_t eventnumbe // Add hits to the cluster object as associated objects for (unsigned int j=0;jchannel(clusterHits[j].row, + clusterHits[j].column); if (npeaks==1){ - myCluster->AddAssociatedObject(clusterHits[j]); + myCluster->addHit(clusterHits[j].id,channel,clusterHits[j].E, + clusterHits[j].x,clusterHits[j].y, + clusterHits[j].t); } else{ // Output hits surrounding peak position - double dx=clusterHits[j]->x-clusterHits[jmax]->x; - double dy=clusterHits[j]->y-clusterHits[jmax]->y; + double dx=clusterHits[j].x-clusterHits[jmax].x; + double dy=clusterHits[j].y-clusterHits[jmax].y; double d=dFCALGeom->blockSize(); if (dFCALGeom->inInsert(channel)){ d=dFCALGeom->insertBlockSize(); } double dcut=2.5*d; if (fabs(dx)AddAssociatedObject(clusterHits[j]); + myCluster->addHit(clusterHits[j].id,channel,clusterHits[j].E, + clusterHits[j].x,clusterHits[j].y, + clusterHits[j].t); } } } @@ -597,9 +628,9 @@ jerror_t DFCALCluster_factory_Island::fini(void) // and all adjacent blocks, starting with the block with the highest energy. // Each successive cluster candidate starts with the block with the next // highest energy that has not already been included in another cluster. -void DFCALCluster_factory_Island::FindClusterCandidates(vector&fcal_hits, - vector>&clusterCandidates) const { - vectorclusterCandidate; +void DFCALCluster_factory_Island::FindClusterCandidates(vector&fcal_hits, + vector>&clusterCandidates) const { + vectorclusterCandidate; clusterCandidate.push_back(fcal_hits[0]); vectorused_hits(fcal_hits.size()); @@ -611,11 +642,11 @@ void DFCALCluster_factory_Island::FindClusterCandidates(vector& for (size_t i=begin_i;it-fcal_hits[i]->t; + double dt=clusterCandidate[0].t-fcal_hits[i].t; if (fabs(dt)>TIME_CUT) continue; - int drow=clusterCandidate[0]->row-fcal_hits[i]->row; - int dcol=clusterCandidate[0]->column-fcal_hits[i]->column; + int drow=clusterCandidate[0].row-fcal_hits[i].row; + int dcol=clusterCandidate[0].column-fcal_hits[i].column; if (abs(drow)<=1 && abs(dcol)<=1){ clusterCandidate.push_back(fcal_hits[i]); @@ -644,20 +675,20 @@ void DFCALCluster_factory_Island::FindClusterCandidates(vector& // Merge cluster candidates that appear belong to a single cluster, again // looking for adjacent hits. double borderCut=0.5*dFCALGeom->blockSize()+dFCALGeom->insertBlockSize(); - vector>::iterator iter=clusterCandidates.begin(); + vector>::iterator iter=clusterCandidates.begin(); while (iter!=clusterCandidates.end()){ bool matched=false; - vector>::iterator iter2=iter+1; + vector>::iterator iter2=iter+1; for (;iter2!=clusterCandidates.end();iter2++){ for (unsigned int i=0;i<(*iter).size();i++){ for (unsigned int j=0;j<(*iter2).size();j++){ - double dt=(*iter)[i]->t-(*iter2)[j]->t; + double dt=(*iter)[i].t-(*iter2)[j].t; if (fabs(dt)>TIME_CUT) continue; - int row1=(*iter)[i]->row; - int col1=(*iter)[i]->column; - int row2=(*iter2)[j]->row; - int col2=(*iter2)[j]->column; + int row1=(*iter)[i].row; + int col1=(*iter)[i].column; + int row2=(*iter2)[j].row; + int col2=(*iter2)[j].column; if (abs(row1-row2)<=1 && abs(col1-col2)<=1){ matched=true; @@ -668,8 +699,8 @@ void DFCALCluster_factory_Island::FindClusterCandidates(vector& // if present if (MERGE_HITS_AT_BOUNDARY && dFCALGeom->hitPairHasInsertHit(row1,row2)){ - double dx=(*iter)[i]->x-(*iter2)[j]->x; - double dy=(*iter)[i]->y-(*iter2)[j]->y; + double dx=(*iter)[i].x-(*iter2)[j].x; + double dy=(*iter)[i].y-(*iter2)[j].y; if (fabs(dx)& // Fit peaks within a cluster containing a list of hits contained in hitList bool DFCALCluster_factory_Island::FitPeaks(const TMatrixD &W,double b, - vector&hitList, + vector&hitList, vector&peaks, PeakInfo &myNewPeak,double &chisq, unsigned int &ndf @@ -766,7 +797,7 @@ bool DFCALCluster_factory_Island::FitPeaks(const TMatrixD &W,double b, Ecalc+=myNewPeak.E*df_dE; - double Ediff=hitList[i]->E-Ecalc; + double Ediff=hitList[i].E-Ecalc; dE(i,0)=Ediff; //cout << " Ediff " << Ediff << endl; } @@ -826,7 +857,7 @@ bool DFCALCluster_factory_Island::FitPeaks(const TMatrixD &W,double b, } vectorhit_channels(nhits); for (unsigned int i=0;ichannel(hitList[i]->row,hitList[i]->column); + hit_channels[i]=dFCALGeom->channel(hitList[i].row,hitList[i].column); } for (unsigned int j=0;jblockSize(); - if (dFCALGeom->isInsertBlock(hit->row,hit->column)){ + if (dFCALGeom->isInsertBlock(hit.row,hit.column)){ half_block=0.5*dFCALGeom->insertBlockSize(); } double b2=b*b; double f=0.; - double dxc=hit->x-myPeakInfo.x; - double dyc=hit->y-myPeakInfo.y; + double dxc=hit.x-myPeakInfo.x; + double dyc=hit.y-myPeakInfo.y; double Ec=myPeakInfo.E; for (int i=0;i<4;i++){ double dx=dxc+sign1[i]*half_block; @@ -875,17 +906,17 @@ double DFCALCluster_factory_Island::CalcClusterXYDeriv(bool isXDeriv,double b, // normalization of this function in this equation does not appear to be // correct (it's off by 1/sqrt(2pi)). double DFCALCluster_factory_Island::CalcClusterEDeriv(double b, - const DFCALHit *hit, + const HitInfo &hit, const PeakInfo &myPeakInfo) const { double sign1[4]={1,1,-1,-1}; double sign2[4]={1,-1,-1,1}; double half_block=0.5*dFCALGeom->blockSize(); - if (dFCALGeom->isInsertBlock(hit->row,hit->column)){ + if (dFCALGeom->isInsertBlock(hit.row,hit.column)){ half_block=0.5*dFCALGeom->insertBlockSize(); } double f=0.; - double dxc=hit->x-myPeakInfo.x; - double dyc=hit->y-myPeakInfo.y; + double dxc=hit.x-myPeakInfo.x; + double dyc=hit.y-myPeakInfo.y; for (int i=0;i<4;i++){ double dx=dxc+sign1[i]*half_block; double dy=dyc+sign2[i]*half_block; @@ -897,7 +928,7 @@ double DFCALCluster_factory_Island::CalcClusterEDeriv(double b, // Try to split peaks into two following the algorithm (barely) described in // Lednev, IHEP 93-153. void DFCALCluster_factory_Island::SplitPeaks(const TMatrixD &W,double b, - vector&hits, + vector&hits, vector&peaks, double &chisq,unsigned int &ndf) const{ unsigned int npeaks=peaks.size(),nhits=hits.size(); @@ -915,10 +946,10 @@ void DFCALCluster_factory_Island::SplitPeaks(const TMatrixD &W,double b, if (i==k) continue; Ecalc+=peaks[k].E*CalcClusterEDeriv(b,hits[j],peaks[k]); } - Elist[j]=hits[j]->E-Ecalc; + Elist[j]=hits[j].E-Ecalc; E0+=Elist[j]; - x0+=Elist[j]*hits[j]->x; - y0+=Elist[j]*hits[j]->y; + x0+=Elist[j]*hits[j].x; + y0+=Elist[j]*hits[j].y; } x0/=E0; y0/=E0; @@ -926,8 +957,8 @@ void DFCALCluster_factory_Island::SplitPeaks(const TMatrixD &W,double b, // Find the second moments about the center of the cluster of hits double xx=0.,yy=0.,yx=0.; for (unsigned int j=0;jx-x0; - double dy=hits[j]->y-y0; + double dx=hits[j].x-x0; + double dy=hits[j].y-y0; xx+=Elist[j]*dx*dx; yy+=Elist[j]*dy*dy; yx+=Elist[j]*dx*dy; @@ -945,7 +976,7 @@ void DFCALCluster_factory_Island::SplitPeaks(const TMatrixD &W,double b, double r = sqrt(dxc*dxc + dyc*dyc); double alpha = 0.; for(unsigned int i=0;ix-x0)*dxc/r + (hits[i]->y-y0)*dyc/r; + double u=(hits[i].x-x0)*dxc/r + (hits[i].y-y0)*dyc/r; alpha-=Elist[i]*u*fabs(u); } alpha/=E0*rsq; From 44230939e65be1c21d702d9bcaeba40f5fa9c69b Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 22 Aug 2024 17:40:46 -0400 Subject: [PATCH 13/19] Modify getUVFromHits, getE1925FromHits, and GetLogWeightedPosition to use the DFCALClusterHit_t vector so that they can work for both ECAL and FCAL hits. --- src/libraries/FCAL/DFCALShower_factory.cc | 100 +++++++++++----------- src/libraries/FCAL/DFCALShower_factory.h | 8 +- 2 files changed, 56 insertions(+), 52 deletions(-) diff --git a/src/libraries/FCAL/DFCALShower_factory.cc b/src/libraries/FCAL/DFCALShower_factory.cc index 0f0a86d8c..4468cba88 100644 --- a/src/libraries/FCAL/DFCALShower_factory.cc +++ b/src/libraries/FCAL/DFCALShower_factory.cc @@ -512,14 +512,14 @@ jerror_t DFCALShower_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber) shower->setTimeTrack( timeTr ); // now compute some variables at the hit level + unsigned int num_hits=cluster->GetNHits(); + shower->setNumBlocks(num_hits); - vector< const DFCALHit* > fcalHits; - cluster->Get( fcalHits ); - shower->setNumBlocks( fcalHits.size() ); - - double e9e25, e1e9; - getE1925FromHits( e1e9, e9e25, fcalHits, - getMaxHit(cluster->getChannelEmax(),fcalHits) ); + // Get (E,x,y) for each hit in the cluster + const vectorhits=cluster->GetHits(); + + double e9e25, e1e9; + getE1925FromHits(hits, e1e9, e9e25); shower->setE1E9( e1e9 ); shower->setE9E25( e9e25 ); @@ -528,7 +528,7 @@ jerror_t DFCALShower_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber) // if there is no nearest track, the defaults for xTr and yTr will result // in using the beam axis as the directional axis //if (!SHOWER_POSITION_LOG) - getUVFromHits( sumU, sumV, fcalHits, + getUVFromHits( sumU, sumV, hits, DVector3( shower->getPosition().X(), shower->getPosition().Y(), 0 ), DVector3( xTr, yTr, 0 ) ); //else @@ -961,7 +961,7 @@ DFCALShower_factory::getMaxHit( const vector< const DFCALHit* >& hitVec ) const void DFCALShower_factory::getUVFromHits( double& sumUSh, double& sumVSh, - const vector< const DFCALHit* >& hits, + const vector& hits, const DVector3& showerVec, const DVector3& trackVec ) const { @@ -982,47 +982,61 @@ DFCALShower_factory::getUVFromHits( double& sumUSh, double& sumVSh, double sumE = 0; - for( vector< const DFCALHit* >::const_iterator hit = hits.begin(); + for( vector::const_iterator hit = hits.begin(); hit != hits.end(); ++hit ){ - hitLoc.SetX( (**hit).x - showerVec.X() ); - hitLoc.SetY( (**hit).y - showerVec.Y() ); + hitLoc.SetX( (*hit).x - showerVec.X() ); + hitLoc.SetY( (*hit).y - showerVec.Y() ); - sumUSh += (**hit).E * pow( u.Dot( hitLoc ), 2 ); - sumVSh += (**hit).E * pow( v.Dot( hitLoc ), 2 ); + sumUSh += (*hit).E * pow( u.Dot( hitLoc ), 2 ); + sumVSh += (*hit).E * pow( v.Dot( hitLoc ), 2 ); - sumE += (**hit).E; + sumE += (*hit).E; } sumUSh /= sumE; sumVSh /= sumE; } -void -DFCALShower_factory::getE1925FromHits( double& e1e9Sh, double& e9e25Sh, - const vector< const DFCALHit* >& hits, - unsigned int maxIndex ) const { - +void DFCALShower_factory::getE1925FromHits(const vector&hits, + double& e1e9Sh, double& e9e25Sh) const { + unsigned int maxIndex=0; + double maxE=0; + for (unsigned int i=0;imaxE){ + maxE=hits[i].E; + maxIndex=i; + } + } + double E9 = 0; double E25 = 0; + double E9cut=0.,E25cut=0.; - const DFCALHit* maxHit = hits[maxIndex]; - - for( vector< const DFCALHit* >::const_iterator hit = hits.begin(); - hit != hits.end(); ++hit ){ - - if( fabs( (**hit).x - maxHit->x ) < 4.5 && fabs( (**hit).y - maxHit->y ) < 4.5 ) - E9 += (**hit).E; - - if( fabs( (**hit).x - maxHit->x ) < 8.5 && fabs( (**hit).y - maxHit->y ) < 8.5 ) - E25 += (**hit).E; + const DFCALCluster::DFCALClusterHit_t maxHit = hits[maxIndex]; + bool maxHitInInsert=fcalGeom->inInsert(maxHit.ch); + if (maxHitInInsert){ + E9cut=2.3; + E25cut=4.3; + } + else { + E9cut=4.5; + E25cut=8.5; + } + + for( vector::const_iterator hit = hits.begin(); + hit != hits.end(); ++hit ){ + if(fabs((*hit).x - maxHit.x) < E9cut && fabs((*hit).y - maxHit.y) < E9cut ) + E9 += (*hit).E; + if(fabs((*hit).x - maxHit.x) < E25cut && fabs((*hit).y - maxHit.y) < E25cut) + E25 += (*hit).E; } - e1e9Sh = maxHit->E/E9; + e1e9Sh = maxE/E9; e9e25Sh = E9/E25; + } - vector< const DTrackWireBased* > DFCALShower_factory::filterWireBasedTracks( vector< const DTrackWireBased* >& wbTracks ) const { @@ -1078,15 +1092,8 @@ void DFCALShower_factory::GetLogWeightedPosition( const DFCALCluster* cluster, D DVector3 posInCal = cluster->getCentroid(); - vector locHitVector; - cluster->Get(locHitVector); - - int loc_nhits = (int)locHitVector.size(); - if( loc_nhits < 1 ) { - pos_log = posInCal; - return; - } - + vector locHitVector=cluster->GetHits(); + //------ Loop over hits ------// double sW = 0.0; @@ -1096,13 +1103,11 @@ void DFCALShower_factory::GetLogWeightedPosition( const DFCALCluster* cluster, D double ecluster = cluster->getEnergy(); - for( int ih = 0; ih < loc_nhits; ih++ ) { - - const DFCALHit *locHit = locHitVector[ih]; + for( int ih = 0; ih < (int)locHitVector.size(); ih++ ) { - double xcell = locHit->x; - double ycell = locHit->y; - double ecell = locHit->E; + double xcell = locHitVector[ih].x; + double ycell = locHitVector[ih].y; + double ecell = locHitVector[ih].E; W = log_position_const + log( ecell / ecluster ); if( W > 0. ) { @@ -1110,7 +1115,6 @@ void DFCALShower_factory::GetLogWeightedPosition( const DFCALCluster* cluster, D xpos += xcell * W; ypos += ycell * W; } - } double x1, y1; diff --git a/src/libraries/FCAL/DFCALShower_factory.h b/src/libraries/FCAL/DFCALShower_factory.h index 282d4abbb..8bbb0e483 100644 --- a/src/libraries/FCAL/DFCALShower_factory.h +++ b/src/libraries/FCAL/DFCALShower_factory.h @@ -19,6 +19,7 @@ #include class DFCALHit; +class DECALHit; class DTrackWireBased; class DFCALShower_factory:public JFactory{ @@ -47,13 +48,12 @@ class DFCALShower_factory:public JFactory{ unsigned int getMaxHit(int chan_Emax, const vector< const DFCALHit* >& hitVec ) const; void getUVFromHits( double& sumUSh, double& sumVSh, - const vector< const DFCALHit* >& hits, + const vector& hits, const DVector3& showerVec, const DVector3& trackVec ) const; - void getE1925FromHits( double& e1e9Sh, double& e9e25Sh, - const vector< const DFCALHit* >& hits, - unsigned int maxIndex ) const; + void getE1925FromHits(const vector&hits, + double& e1e9Sh, double& e9e25Sh) const; vector< const DTrackWireBased* > filterWireBasedTracks( vector< const DTrackWireBased* >& wbTracks ) const; From ef9e3c3ad814aeb340892c9205a3a6b8d2a34725 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Mon, 26 Aug 2024 08:59:03 -0400 Subject: [PATCH 14/19] Change how hits are sorted so that ECAL hits can be included in the sort. --- src/libraries/FCAL/DFCALCluster_factory.cc | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/libraries/FCAL/DFCALCluster_factory.cc b/src/libraries/FCAL/DFCALCluster_factory.cc index 9c588fdef..f5a75e541 100644 --- a/src/libraries/FCAL/DFCALCluster_factory.cc +++ b/src/libraries/FCAL/DFCALCluster_factory.cc @@ -23,8 +23,11 @@ using namespace jana; // Used to sort hits by Energy -bool FCALHitsSort_C(const DFCALHit* const &thit1, const DFCALHit* const &thit2) { - return thit1->E>thit2->E; +int FCALHitsSort_C(const void *a,const void *b){ + const DFCALCluster::userhit_t hit1=*(const DFCALCluster::userhit_t *)a; + const DFCALCluster::userhit_t hit2=*(const DFCALCluster::userhit_t *)b; + if (hit1.E &fcalhits) { @@ -103,9 +106,6 @@ jerror_t DFCALCluster_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber) const DFCALGeometry* fcalGeom=NULL; eventLoop->GetSingle(fcalGeom); - // Sort hits by energy - sort(fcalhits.begin(), fcalhits.end(), FCALHitsSort_C); - // fill user's hit list int nhits = 0; DFCALCluster::userhits_t* hits = @@ -162,7 +162,10 @@ jerror_t DFCALCluster_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber) for ( int i = 0; i < nhits; i++ ) { hitUsed[i] = 0; } - + + // Sort the hit array by energy + qsort(&hits->hit[0],nhits,sizeof(DFCALCluster::userhit_t),FCALHitsSort_C); + const unsigned int max = 999; DFCALCluster* clusterList[max]; unsigned int clusterCount = 0; From 571bc7af7ac0880b18b5859f0e08dfc59ab6aae6 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Mon, 26 Aug 2024 09:00:54 -0400 Subject: [PATCH 15/19] Update GetECALZ to be consistent with new xml geometry. --- src/libraries/HDGEOMETRY/DGeometry.cc | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/src/libraries/HDGEOMETRY/DGeometry.cc b/src/libraries/HDGEOMETRY/DGeometry.cc index 265ce8df5..9468bf818 100644 --- a/src/libraries/HDGEOMETRY/DGeometry.cc +++ b/src/libraries/HDGEOMETRY/DGeometry.cc @@ -1768,20 +1768,23 @@ bool DGeometry::GetBCALPhiShift(float &bcal_phi_shift) const //--------------------------------- bool DGeometry::GetECALZ(double &z_ecal) const { - vector CrystalEcalpos; - jgeom->SetVerbose(0); // don't print error messages for optional detector elements - bool good = Get("//section/composition/posXYZ[@volume='CrystalECAL']/@X_Y_Z", CrystalEcalpos); - jgeom->SetVerbose(1); // reenable error messages + z_ecal=0; + if (GetFCALZ(z_ecal)){ + jgeom->SetVerbose(0); // don't print error messages for optional detector elements + vectorCrystalEcalpos; + if (Get("//section/composition/posXYZ[@volume='CrystalECAL']/@X_Y_Z", CrystalEcalpos)){ + vector FCALCenter; + Get("//section/composition/posXYZ[@volume='forwardEMcal]/@X_Y_Z", FCALCenter); + vectorblock; + Get("//box[@name='XTBL']/@X_Y_Z",block); - if(!good){ - // NEED TO RETHINK ERROR REPORTING FOR OPTIONAL DETECTOR ELEMENTS - //_DBG_<<"Unable to retrieve ECAL position."<SetVerbose(1); // reenable error messages + + z_ecal += FCALCenter[2]+CrystalEcalpos[2]-0.5*block[2]; return true; - } + } + } + return false; } From e0c4ba616e0907490e2832d9933551d08e771ca2 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Tue, 27 Aug 2024 16:28:28 -0400 Subject: [PATCH 16/19] Add routine to get the size of the FCAL insert. --- src/libraries/HDGEOMETRY/DGeometry.cc | 18 +++++++++++++++--- src/libraries/HDGEOMETRY/DGeometry.h | 1 + 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/src/libraries/HDGEOMETRY/DGeometry.cc b/src/libraries/HDGEOMETRY/DGeometry.cc index 9468bf818..b8e2cf68c 100644 --- a/src/libraries/HDGEOMETRY/DGeometry.cc +++ b/src/libraries/HDGEOMETRY/DGeometry.cc @@ -1787,7 +1787,6 @@ bool DGeometry::GetECALZ(double &z_ecal) const return false; } - //--------------------------------- // GetCCALZ //--------------------------------- @@ -2088,8 +2087,6 @@ bool DGeometry::GetFCALInsertRowSize(int &insert_row_size) const jgeom->SetVerbose(1); // reenable error messages if(!good){ - // NEED TO RETHINK ERROR REPORTING FOR OPTIONAL DETECTOR ELEMENTS - //_DBG_<<"Unable to retrieve ComptonEMcal position."<insertBlock; + // Use nominal crystal separation to be on the safe side. This is slightly + // larger than the actual measured separations from the survey. + return 0.5*double(insertRowSize)*2.09; + } + + return 0.; +} + //--------------------------------- // GetFCALBlockSize //--------------------------------- diff --git a/src/libraries/HDGEOMETRY/DGeometry.h b/src/libraries/HDGEOMETRY/DGeometry.h index 9ce06fe28..2f073621d 100644 --- a/src/libraries/HDGEOMETRY/DGeometry.h +++ b/src/libraries/HDGEOMETRY/DGeometry.h @@ -194,6 +194,7 @@ class DGeometry{ bool GetFCALBlockSize(vector &block) const; bool GetFCALInsertBlockSize(vector &block) const; bool HaveInsert() const; + double GetFCALInsertSize() const; bool GetStartCounterGeom(vector >&pos, vector >&norm) const; // < vectors containing positions and norm 3-vectors for start counter From 09397247133bfb57abbb60d79a2c3ef1a2a62d71 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Tue, 27 Aug 2024 16:29:03 -0400 Subject: [PATCH 17/19] use new routine to get FCAL insert size --- src/libraries/FCAL/DFCALGeometry.cc | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/libraries/FCAL/DFCALGeometry.cc b/src/libraries/FCAL/DFCALGeometry.cc index 384b26228..55a33c6ab 100644 --- a/src/libraries/FCAL/DFCALGeometry.cc +++ b/src/libraries/FCAL/DFCALGeometry.cc @@ -27,9 +27,7 @@ DFCALGeometry::DFCALGeometry(const DGeometry *geom){ // Check for presence of PbWO4 insert if (geom->HaveInsert()){ // Geometry based on survey data for positions // Get size of the insert - int insert_row_size=0; - geom->GetFCALInsertRowSize(insert_row_size); - m_insertSize=insertBlockSize()*double(insert_row_size/2); + m_insertSize=geom->GetFCALInsertSize(); // Find the size of the sensitive volume of each PWO crystal geom->Get("//box[@name='XTBL']/@X_Y_Z",block); From f4dfb8aea499b811d2123bdb45d9961fef4a3706 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Tue, 27 Aug 2024 16:35:21 -0400 Subject: [PATCH 18/19] Add routine to add the extrapolation point at the entrance of the FCAL and at the exit of the block that can be used for both FCAL-1 and FCAL-2 --- .../TRACKING/DTrackFitterKalmanSIMD.cc | 145 ++++++++++++------ .../TRACKING/DTrackFitterKalmanSIMD.h | 6 + 2 files changed, 101 insertions(+), 50 deletions(-) diff --git a/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc b/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc index 6a3358332..149c810ab 100644 --- a/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc +++ b/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc @@ -330,8 +330,16 @@ DTrackFitterKalmanSIMD::DTrackFitterKalmanSIMD(JEventLoop *loop):DTrackFitter(lo +0.5*cdc_endplate_dim[2]; // Outer detector geometry parameters - geom->GetFCALZ(dFCALz); + geom->GetFCALPosition(dFCALx,dFCALy,dFCALz); dFCALzBack=dFCALz+45.; + // The following initializes ECAL position variables to large values so that + // we don't make extroplations to the FCAL insert face if it is not present + dECALz=dECALzBack=1000.; + dECALsize=0.; + if (geom->GetECALZ(dECALz)){ + dECALzBack=dECALz+20.; + dECALsize=geom->GetFCALInsertSize(); + } if (geom->GetDIRCZ(dDIRCz)==false) dDIRCz=1000.; geom->GetFMWPCZ_vec(dFMWPCz_vec); geom->GetFMWPCSize(dFMWPCsize); @@ -8696,11 +8704,11 @@ jerror_t DTrackFitterKalmanSIMD::ExtrapolateToOuterDetectors(const DMatrix5x1 &S const double z_outer_max=1000.; const double x_max=130.; const double y_max=130.; - const double fcal_radius_sq=120.47*120.47; bool hit_tof=false; bool hit_dirc=false; bool hit_fcal=false; bool got_fmwpc=(dFMWPCz_vec.size()>0)?true:false; + bool got_ecal=(dECALz<1000.)?true:false; unsigned int fmwpc_index=0; unsigned int trd_index=0; while (z>Z_MIN && zdECALz){ + newz=dECALz+EPS; + ds=(newz-z)/dz_ds; + } if (fmwpc_indexdFMWPCz_vec[fmwpc_index]){ newz=dFMWPCz_vec[fmwpc_index]+EPS; ds=(newz-z)/dz_ds; @@ -8819,58 +8831,31 @@ jerror_t DTrackFitterKalmanSIMD::ExtrapolateToOuterDetectors(const DMatrix5x1 &S } if (hit_fcal==false && newz>dFCALz){ double r2=S(state_x)*S(state_x)+S(state_y)*S(state_y); - if (r2>fcal_radius_sq) return NOERROR; + if (r2>dFCALradiusSq) return NOERROR; hit_fcal=true; - AddExtrapolation(SYS_FCAL,z,S,t,s); - - // Propagate the track to the back of the FCAL block - int num=int(45./dz); - int m=0; - for (;mBIG) one_over_beta2=BIG; - ds=dz*sqrt(1.+S(state_tx)*S(state_tx)+S(state_ty)*S(state_ty)); - t+=ds*sqrt(one_over_beta2); // in units where c=1 - s+=ds; - - newz=z+dz; - // Step through field - Step(z,newz,dEdx,S); - z=newz; - - r2=S(state_x)*S(state_x)+S(state_y)*S(state_y); - if (r2>fcal_radius_sq){ - // Break out of the loop if the track exits out of the FCAL before - // reaching the back end of the block - break; - } - } - if (m==num){ - newz=dFCALzBack; - - // Propagate t and s to back of FCAL block - double q_over_p_sq=S(state_q_over_p)*S(state_q_over_p); - double one_over_beta2=1.+mass2*q_over_p_sq; - if (one_over_beta2>BIG) one_over_beta2=BIG; - dz=newz-z; - ds=dz*sqrt(1.+S(state_tx)*S(state_tx)+S(state_ty)*S(state_ty)); - t+=ds*sqrt(one_over_beta2); // in units where c=1 - s+=ds; - - // Step through field - Step(z,newz,dEdx,S); - z=newz; - } - // add another extrapolation point at downstream end of FCAL - AddExtrapolation(SYS_FCAL,z,S,t,s); - // .. and exit if the muon detector is not installed - if (got_fmwpc==false) return NOERROR; + // We've reached the front face of the FCAL. Now we need to determine if + // we need to follow the track further in z to match to the ECAL instead + // if the insert is present. + if (got_ecal==false + || (got_ecal==true && (fabs(S(state_x)-dFCALx)>dECALsize + || fabs(S(state_y)-dFCALy)>dECALsize))){ + PropagateThroughFCAL(SYS_FCAL,dEdx,dz,z,S,t,s); + } + // exit if the muon detector or ECAL is not installed + if (got_ecal==false && got_fmwpc==false) return NOERROR; + } + if (got_ecal==true && newz>dECALz){ + if (fabs(S(state_x)-dFCALx)>dECALsize + || fabs(S(state_y)-dFCALy)>dECALsize) return NOERROR; + PropagateThroughFCAL(SYS_ECAL,dEdx,dz,z,S,t,s); + + // .. and exit + return NOERROR; } + // Deal with muon detector - if (hit_fcal==true + if (hit_fcal==true && got_ecal==false && (fabs(S(state_x))>dFMWPCsize || (fabs(S(state_y))>dFMWPCsize))){ return NOERROR; } @@ -10250,3 +10235,63 @@ void DTrackFitterKalmanSIMD::AddExtrapolation(DetectorSystem_t detector, << endl; } } + +// Propagate the track in the forward direction through the FCAL +void DTrackFitterKalmanSIMD::PropagateThroughFCAL(const DetectorSystem_t detector,const double dEdx,const double dz,double &z,DMatrix5x1 &S,double &t,double &s){ + AddExtrapolation(SYS_FCAL,z,S,t,s); + + // Propagate the track to the back of the FCAL/ECAL block + double length=(detector==SYS_FCAL)?45.:20.; + int num=int(length/dz); + int m=0; + for (;mBIG) break; + double ds=dz*sqrt(1.+S(state_tx)*S(state_tx)+S(state_ty)*S(state_ty)); + t+=ds*sqrt(one_over_beta2); // in units where c=1 + s+=ds; + + double newz=z+dz; + // Step through field + Step(z,newz,dEdx,S); + z=newz; + + if (detector==SYS_FCAL){ + double r2=S(state_x)*S(state_x)+S(state_y)*S(state_y); + if (r2>dFCALradiusSq){ + // Break out of the loop if the track exits out of the FCAL before + // reaching the back end of the block + break; + } + // Break out of the loop if the track would go into the insert region + // if ECAL is present + if (fabs(S(state_x)-dFCALx)dECALsize + || fabs(S(state_y)-dFCALy)>dECALsize) break; + } + } + if (m==num){ + double newz=(detector==SYS_FCAL)?dFCALzBack:dECALzBack; + + // Propagate t and s to back of FCAL/ECAL block + double q_over_p_sq=S(state_q_over_p)*S(state_q_over_p); + double one_over_beta2=1.+mass2*q_over_p_sq; + if (one_over_beta2 Get7x7ErrorMatrix(DMatrixDSym C); shared_ptr Get7x7ErrorMatrixForward(DMatrixDSym C); @@ -481,6 +484,9 @@ class DTrackFitterKalmanSIMD: public DTrackFitter{ vectorcdc_origin; // outer detectors double dTOFz,dFCALz,dFCALzBack,dDIRCz,dFMWPCsize,dCTOFz; + double dFCALx,dFCALy; + double dFCALradiusSq=120.47*120.47; + double dECALz,dECALzBack,dECALsize; vectordFMWPCz_vec; vectordTRDz_vec; From 68c7d77689464d4b9b87b59cf233b0383326e286 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Mon, 2 Sep 2024 08:33:39 -0400 Subject: [PATCH 19/19] Return early from brun if TRIG:BYPASS is set to 1. Get number of FCAL channels from FCALGeometry. --- src/libraries/TRIGGER/DL1MCTrigger_factory.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/libraries/TRIGGER/DL1MCTrigger_factory.cc b/src/libraries/TRIGGER/DL1MCTrigger_factory.cc index cbf3be526..ed91e165d 100644 --- a/src/libraries/TRIGGER/DL1MCTrigger_factory.cc +++ b/src/libraries/TRIGGER/DL1MCTrigger_factory.cc @@ -168,6 +168,7 @@ jerror_t DL1MCTrigger_factory::init(void) //------------------ jerror_t DL1MCTrigger_factory::brun(jana::JEventLoop *eventLoop, int32_t runnumber) { + if (BYPASS) return NOERROR; int use_rcdb = 1; @@ -1479,7 +1480,7 @@ int DL1MCTrigger_factory::FindTriggers(DL1MCTrigger *trigger, vector &fcal_const_ch, const DFCALGeometry &fcalGeom){ - for (int ch = 0; ch < static_cast(fcal_const_ch.size()); ch++) { + for (int ch = 0; ch < fcalGeom.numFcalChannels(); ch++) { int row = fcalGeom.row(ch); int col = fcalGeom.column(ch); table[row][col] = fcal_const_ch[ch];