diff --git a/podd b/podd
index cc0ebfe4b74446e27dc5a7040dfe2b4606b24f9e..1154ca9087e06ede79833e8ce42cc04cbb99f479 160000
--- a/podd
+++ b/podd
@@ -1 +1 @@
-Subproject commit cc0ebfe4b74446e27dc5a7040dfe2b4606b24f9e
+Subproject commit 1154ca9087e06ede79833e8ce42cc04cbb99f479
diff --git a/src/THcAerogel.cxx b/src/THcAerogel.cxx
index da2716ecd4505b8557eb492a0416635cbd58930f..610635526567a8ea260e52bf7ff57b5ac0c839b5 100644
--- a/src/THcAerogel.cxx
+++ b/src/THcAerogel.cxx
@@ -19,6 +19,7 @@
#include "THaCutList.h"
#include "THcParmList.h"
#include "THcHitList.h"
+#include "THaApparatus.h"
#include "VarDef.h"
#include "VarType.h"
#include "THaTrack.h"
diff --git a/src/THcCherenkov.cxx b/src/THcCherenkov.cxx
index 236db7c8be1b39244d752dde54a5c2ffd2a2570d..31a3418027ffead7028a3198e1345e62a380520d 100644
--- a/src/THcCherenkov.cxx
+++ b/src/THcCherenkov.cxx
@@ -32,6 +32,7 @@
#include "THaCutList.h"
#include "THcParmList.h"
#include "THcHitList.h"
+#include "THaApparatus.h"
#include "VarDef.h"
#include "VarType.h"
#include "THaTrack.h"
diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index 09ab985c4182268819df4a25fd4d64ad092aa154..f1e8290bdd9c7ea80e97b3b4fe731ece363a055a 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -23,6 +23,7 @@
#include "TClonesArray.h"
#include "TMath.h"
#include "TVectorD.h"
+#include "THaApparatus.h"
#include <cstring>
#include <cstdio>
diff --git a/src/THcDriftChamber.cxx b/src/THcDriftChamber.cxx
index 61f5c02c76029a2d72c04ee6863822dd009e1fdf..b9ed691797b48154cf905f7ca2d4330b3d7bed22 100644
--- a/src/THcDriftChamber.cxx
+++ b/src/THcDriftChamber.cxx
@@ -21,6 +21,7 @@
#include "TMath.h"
#include "TVectorD.h"
#include "THcSpacePoint.h"
+#include "THaApparatus.h"
#include "THaTrackProj.h"
diff --git a/src/THcDriftChamberPlane.cxx b/src/THcDriftChamberPlane.cxx
index 59516413f569d69fc356f030059fb0b0b19337cf..3df5a7bd50a2b073a5f6f7cbbe76f705cab3cbec 100644
--- a/src/THcDriftChamberPlane.cxx
+++ b/src/THcDriftChamberPlane.cxx
@@ -6,6 +6,7 @@
//
//////////////////////////////////////////////////////////////////////////
+#include "THcDC.h"
#include "THcDriftChamberPlane.h"
#include "THcDCWire.h"
#include "THcDCHit.h"
@@ -14,10 +15,11 @@
#include "THcGlobals.h"
#include "THcParmList.h"
#include "THcHitList.h"
-#include "THcDC.h"
+#include "THaApparatus.h"
#include "THcHodoscope.h"
#include "TClass.h"
+
#include <cstring>
#include <cstdio>
#include <cstdlib>
diff --git a/src/THcFormula.cxx b/src/THcFormula.cxx
index 7740c06948fb0d75bcd6c1924e4d81c985d25c1e..114d8f6010d19e0b4d277ef390e999e033be295d 100644
--- a/src/THcFormula.cxx
+++ b/src/THcFormula.cxx
@@ -16,46 +16,64 @@
#include "THaVarList.h"
#include "THaCutList.h"
#include "THaCut.h"
+#include "TMath.h"
#include <iostream>
+#include <numeric>
using namespace std;
static const Double_t kBig = 1e38; // Error value
+enum EFuncCode { kLength, kSum, kMean, kStdDev, kMax, kMin,
+ kGeoMean, kMedian, kIteration, kNumSetBits };
+
+#define ALL(c) (c).begin(), (c).end()
+
+//_____________________________________________________________________________
+static inline Int_t NumberOfSetBits( UInt_t v )
+{
+ // Count number of bits set in 32-bit integer. From
+ // http://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel
+
+ v = v - ((v >> 1) & 0x55555555);
+ v = (v & 0x33333333) + ((v >> 2) & 0x33333333);
+ return (((v + (v >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24;
+}
+
+//_____________________________________________________________________________
+static inline Int_t NumberOfSetBits( ULong64_t v )
+{
+ // Count number of bits in 64-bit integer
+
+ const ULong64_t mask32 = (1LL<<32)-1;
+ return NumberOfSetBits( static_cast<UInt_t>(mask32 & v) ) +
+ NumberOfSetBits( static_cast<UInt_t>(mask32 & (v>>32)) );
+}
+
//_____________________________________________________________________________
THcFormula::THcFormula(const char* name, const char* expression,
const THcParmList* plst, const THaVarList* vlst,
const THaCutList* clst ) :
THaFormula()
{
+ Bool_t do_register=0;
- // We have to duplicate the TFormula constructor code here because of
- // the calls DefinedVariable. Our version will only get called if
- // to Compile(). Compile() only works if fParmList is set.
fParmList = plst;
fVarList = vlst;
fCutList = clst;
- SetName(name);
-
- //eliminate blanks in expression
- Int_t nch = strlen(expression);
- char *expr = new char[nch+1];
- Int_t j = 0;
- for (Int_t i=0;i<nch;i++) {
- if (expression[i] == ' ') continue;
- if (i > 0 && (expression[i] == '*') && (expression[i-1] == '*')) {
- expr[j-1] = '^';
- continue;
- }
- expr[j] = expression[i]; j++;
- }
- expr[j] = 0;
- if (j) SetTitle(expr);
- delete [] expr;
+ if( Init( name, expression ) != 0 ) {
+ RegisterFormula(false);
+ return;
+ }
+
+ SetBit(kNotGlobal,!do_register);
Compile(); // This calls our own Compile()
+
+ if( do_register )
+ RegisterFormula();
}
//_____________________________________________________________________________
@@ -63,15 +81,10 @@ THcFormula& THcFormula::operator=( const THcFormula& rhs )
{
if( this != &rhs ) {
TFormula::operator=(rhs);
- fNcodes = rhs.fNcodes;
fParmList = rhs.fParmList;
fVarList = rhs.fVarList;
fCutList = rhs.fCutList;
- fError = rhs.fError;
- fRegister = rhs.fRegister;
- delete [] fVarDef;
- fVarDef = new FVarDef_t[ kMAXCODES ];
- memcpy( fVarDef, rhs.fVarDef, kMAXCODES*sizeof(FVarDef_t));
+ fInstance = 0;
}
return *this;
}
@@ -83,7 +96,7 @@ THcFormula::~THcFormula()
}
//_____________________________________________________________________________
-Int_t THcFormula::DefinedCut( const TString& name )
+Int_t THcFormula::DefinedCut( TString& name )
{
// Check if 'name' is a known cut. If so, enter it in the local list of
// variables used in this formula.
@@ -102,26 +115,23 @@ Int_t THcFormula::DefinedCut( const TString& name )
} else if (attribute.CompareTo("ncalled")==0) {
thistype = (EVariableType) kCutNCalled;
} else {
- thistype = kUndefined;
+ thistype = kCut;
}
}
// Cut names are obviously only valid if there is a list of existing cuts
if( fCutList ) {
- const THaCut* pcut = fCutList->FindCut( realname );
+ THaCut* pcut = fCutList->FindCut( realname );
if( pcut ) {
- if( fNcodes >= kMAXCODES ) return -1;
- // See if this cut already used earlier in this new cut
- FVarDef_t* def = fVarDef;
- for( Int_t i=0; i<fNcodes; i++, def++ ) {
- if( def->type == thistype && pcut == def->code )
+ // See if this cut already used earlier in the expression
+ for( vector<FVarDef_t>::size_type i=0; i<fVarDef.size(); ++i ) {
+ FVarDef_t& def = fVarDef[i];
+ if( def.type == thistype && pcut == def.obj )
return i;
}
- def->type = thistype;
- def->code = pcut;
- def->index = 0;
+ fVarDef.push_back( FVarDef_t(thistype,pcut,0) );
fNpar = 0;
- return fNcodes++;
+ return fVarDef.size()-1;
}
}
return -1;
@@ -135,80 +145,201 @@ Double_t THcFormula::DefinedValue( Int_t i )
// (calculated at the last evaluation).
// If the variable is a string, return value of its character value
-#ifdef WITH_DEBUG
- R__ASSERT( i>=0 && i<fNcodes );
-#endif
- FVarDef_t* def = fVarDef+i;
- const void* ptr = def->code;
- if( !ptr ) return kBig;
- switch( (Int_t) def->type ) {
+ typedef vector<Double_t>::size_type vsiz_t;
+ typedef vector<Double_t>::iterator viter_t;
+
+ assert( i>=0 && i<(Int_t)fVarDef.size() );
+
+ if( IsInvalid() )
+ return 1.0;
+
+ FVarDef_t& def = fVarDef[i];
+ switch( def.type ) {
case kVariable:
case kString:
- return reinterpret_cast<const THaVar*>(ptr)->GetValue( def->index );
+ case kArray:
+ {
+ const THaVar* var = static_cast<const THaVar*>(def.obj);
+ assert(var);
+ Int_t index = (def.type == kArray) ? fInstance : def.index;
+ assert(index >= 0);
+ if( index >= var->GetLen() ) {
+ SetBit(kInvalid);
+ return 1.0; // safer than kBig to prevent overflow
+ }
+ return var->GetValue( index );
+ }
break;
case kCut:
- return reinterpret_cast<const THaCut*>(ptr)->GetResult();
+ {
+ const THaCut* cut = static_cast<const THaCut*>(def.obj);
+ assert(cut);
+ return cut->GetResult();
+ }
break;
case kCutScaler:
- return reinterpret_cast<const THaCut*>(ptr)->GetNPassed();
+ {
+ const THaCut* cut = static_cast<const THaCut*>(def.obj);
+ assert(cut);
+ return cut->GetNPassed();
+ }
break;
case kCutNCalled:
- return reinterpret_cast<const THaCut*>(ptr)->GetNCalled();
+ {
+ const THaCut* cut = static_cast<const THaCut*>(def.obj);
+ assert(cut);
+ return cut->GetNCalled();
+ }
+ break;
+ case kFunction:
+ {
+ EFuncCode code = static_cast<EFuncCode>(def.index);
+ switch( code ) {
+ case kIteration:
+ return fInstance;
+ default:
+ assert(false); // not reached
+ break;
+ }
+ }
+ break;
+ case kFormula:
+ case kVarFormula:
+ {
+ EFuncCode code = static_cast<EFuncCode>(def.index);
+ THaFormula* func = static_cast<THaFormula*>(def.obj);
+ assert(func);
+
+ vsiz_t ndata = func->GetNdata();
+ if( code == kLength )
+ return ndata;
+
+ if( ndata == 0 ) {
+ //FIXME: needs thought
+ SetBit(kInvalid);
+ return 1.0;
+ }
+ Double_t y;
+ if( code == kNumSetBits ) {
+ // Number of set bits is intended for unsigned int-type expressions
+ y = func->EvalInstance(fInstance);
+ if( y > kMaxULong64 || y < kMinLong64 ) {
+ return 0;
+ }
+ return NumberOfSetBits( static_cast<ULong64_t>(y) );
+ }
+
+ vector<Double_t> values;
+ values.reserve(ndata);
+ for( vsiz_t i = 0; i < ndata; ++i ) {
+ values.push_back( func->EvalInstance(i) );
+ }
+ if( func->IsInvalid() ) {
+ SetBit(kInvalid);
+ return 1.0;
+ }
+ switch( code ) {
+ case kSum:
+ y = accumulate( ALL(values), static_cast<Double_t>(0.0) );
+ break;
+ case kMean:
+ y = TMath::Mean( ndata, &values[0] );
+ break;
+ case kStdDev:
+ y = TMath::RMS( ndata, &values[0] );
+ break;
+ case kMax:
+ y = *max_element( ALL(values) );
+ break;
+ case kMin:
+ y = *min_element( ALL(values) );
+ break;
+ case kGeoMean:
+ y = TMath::GeomMean( ndata, &values[0] );
+ break;
+ case kMedian:
+ y = TMath::Median( ndata, &values[0] );
+ break;
+ default:
+ assert(false); // not reached
+ break;
+ }
+ return y;
+ }
break;
- default:
- return kBig;
}
+ assert(false); // not reached
+ return kBig;
}
//_____________________________________________________________________________
-Int_t THcFormula::DefinedGlobalVariable( const TString& name )
+Int_t THcFormula::DefinedGlobalVariable( TString& name )
{
// Check if 'name' is a known global variable. If so, enter it in the
// local list of variables used in this formula.
// No list of variables or too many variables in this formula?
- if( (!fVarList && !fParmList) || fNcodes >= kMAXCODES )
- return -2;
-
+ if( !fVarList && !fParmList )
+ return -1;
// Parse name for array syntax
- THaArrayString var(name);
- if( var.IsError() )
- return -1;
+ THaArrayString parsed_name(name);
+ if( parsed_name.IsError() ) return -1;
// First check if this name is a Parameter
- const THaVar* obj = fParmList->Find( var.GetName() );
- if ( !obj) { // If not, find a global variable with this name
- obj = fVarList->Find( var.GetName() );
- if( !obj )
+ THaVar* var = fParmList->Find( parsed_name.GetName() );
+ if ( !var) { // If not, find a global variable with this name
+ var = fVarList->Find( parsed_name.GetName() );
+ if( !var )
return -1;
}
- // Error if array requested but the corresponding variable is not an array
- if( var.IsArray() && !obj->IsArray() )
- return -2;
-
- // Subscript(s) within bounds?
+ EVariableType type = kVariable;
Int_t index = 0;
- if( var.IsArray()
- && (index = obj->Index( var )) <0 ) return -2;
-
+ if( parsed_name.IsArray() ) {
+ // Requesting an array element
+ if( !var->IsArray() )
+ // Element requested but the corresponding variable is not an array
+ return -2;
+ if( var->IsVarArray() ) {
+ // Variable-size arrays are always 1-d
+ assert( var->GetNdim() == 1 );
+ if( parsed_name.GetNdim() != 1 )
+ return -2;
+ // For variable-size arrays, we allow requests for any element and
+ // check at evaluation time whether the element actually exists
+ index = parsed_name[0];
+ } else {
+ // Fixed-size arrays: check if subscript(s) within bounds?
+ //TODO: allow fixing some dimensions
+ index = var->Index( parsed_name );
+ if( index < 0 )
+ return -2;
+ }
+ } else if( var->IsArray() ) {
+ // Asking for an entire array
+ type = kArray;
+ SetBit(kArrayFormula);
+ if( var->IsVarArray() )
+ SetBit(kVarArray);
+ }
+
// Check if this variable already used in this formula
- FVarDef_t* def = fVarDef;
- for( Int_t i=0; i<fNcodes; i++, def++ ) {
- if( obj == def->code && index == def->index )
+ for( vector<FVarDef_t>::size_type i=0; i<fVarDef.size(); ++i ) {
+ FVarDef_t& def = fVarDef[i];
+ if( var == def.obj && index == def.index ) {
+ assert( type == def.type );
return i;
+ }
}
// If this is a new variable, add it to the list
- def->type = kVariable;
- def->code = obj;
- def->index = index;
+ fVarDef.push_back( FVarDef_t(type,var,index) );
// No parameters ever for a THaFormula
fNpar = 0;
- return fNcodes++;
+ return fVarDef.size()-1;
}
diff --git a/src/THcFormula.h b/src/THcFormula.h
index 0294ee6417a21b59cbb21a388ad301d79225eb9a..20c962ac8c19d0ba28666f5ed2b8446f5271d4f2 100644
--- a/src/THcFormula.h
+++ b/src/THcFormula.h
@@ -23,12 +23,12 @@ public:
virtual ~THcFormula();
virtual Double_t DefinedValue( Int_t i);
- virtual Int_t DefinedCut( const TString& variable);
- virtual Int_t DefinedGlobalVariable( const TString& variable);
+ virtual Int_t DefinedCut( TString& variable);
+ virtual Int_t DefinedGlobalVariable( TString& variable);
protected:
- enum {kCutScaler = kString+1};
+ enum {kCutScaler = kVarFormula+1};
enum {kCutNCalled = kCutScaler+1};
const THcParmList* fParmList; // Pointer to list of parameters
ClassDef(THcFormula,0) // Formula with cut scalers
diff --git a/src/THcInterface.cxx b/src/THcInterface.cxx
index e7c1d5b2e7631c6da07be1c5fa9148e655ea5442..53c190eb61d30b0d95fb4457be075ec81b30122f 100644
--- a/src/THcInterface.cxx
+++ b/src/THcInterface.cxx
@@ -22,7 +22,7 @@
#include "THcParmList.h"
#include "THcDetectorMap.h"
#include "THaCutList.h"
-#include "THaCodaDecoder.h"
+#include "CodaDecoder.h"
#include "THaGlobals.h"
#include "THcGlobals.h"
#include "THaAnalyzer.h"
@@ -32,7 +32,7 @@
#include "TTree.h"
-
+#include <iostream>
//#include "TGXW.h"
//#include "TVirtualX.h"
@@ -81,7 +81,7 @@ THcInterface::THcInterface( const char* appClassName, int* argc, char** argv,
gHaScalers = new TList;
gHaPhysics = new TList;
// Use the standard CODA file decoder by default
- gHaDecoder = THaCodaDecoder::Class();
+ gHaDecoder = Decoder::CodaDecoder::Class();
// File-based database by default
// gHaDB = new THaFileDB();
gHaTextvars = new THaTextvars;
diff --git a/src/THcParmList.cxx b/src/THcParmList.cxx
index e6caee667ac2d9aaf70b92790fbd6acdad46678a..23ff5be3e1840bdb24c251e632bc79da77dd6285 100644
--- a/src/THcParmList.cxx
+++ b/src/THcParmList.cxx
@@ -339,8 +339,8 @@ void THcParmList::Load( const char* fname, Int_t RunNumber )
if(valstr.IsFloat()) {
fp[currentindex+i] = valstr.Atof();
} else {
- THaFormula* formula = new THaFormula("temp",valstr.Data()
- ,this, 0);
+ THaFormula* formula = new THaFormula
+ ("temp",valstr.Data(), (Bool_t) 0, this, 0);
fp[currentindex+i] = formula->Eval();
delete formula;
}
@@ -377,8 +377,8 @@ void THcParmList::Load( const char* fname, Int_t RunNumber )
if(valstr.IsFloat()) {
existingp[currentindex+i] = valstr.Atof();
} else {
- THaFormula* formula = new THaFormula("temp",valstr.Data()
- ,this, 0);
+ THaFormula* formula = new THaFormula
+ ("temp",valstr.Data(), (Bool_t) 0, this, 0);
existingp[currentindex+i] = formula->Eval();
delete formula;
}
@@ -403,8 +403,8 @@ void THcParmList::Load( const char* fname, Int_t RunNumber )
if(valstr.IsFloat()) {
fp[i] = valstr.Atof();
} else {
- THaFormula* formula = new THaFormula("temp",valstr.Data()
- ,this, 0);
+ THaFormula* formula = new THaFormula
+ ("temp",valstr.Data(), (Bool_t) 0, this, 0);
fp[i] = formula->Eval();
delete formula;
}
diff --git a/src/THcShowerPlane.cxx b/src/THcShowerPlane.cxx
index 8bd3d32e3174eeb900654d98766965c8a6a3d756..81c9f67adf9a6f1c392f97f5593e2ce6d7e77cba 100644
--- a/src/THcShowerPlane.cxx
+++ b/src/THcShowerPlane.cxx
@@ -93,7 +93,6 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date )
THcShower* fParent;
fParent = (THcShower*) GetParent();
-
// Find the number of elements
fNelem = fParent->GetNBlocks(fLayerNum-1);
@@ -152,7 +151,7 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date )
if (fParent->fdbg_init_cal) {
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShowerPlane::ReadDatabase for "
- << GetApparatus()->GetName() << ":" << endl;
+ << GetParent()->GetPrefix() << ":" << endl;
cout << " Layer #" << fLayerNum << ", number of elements " << fNelem
<< endl;
@@ -213,7 +212,7 @@ void THcShowerPlane::Clear( Option_t* )
if ( ((THcShower*) GetParent())->fdbg_decoded_cal ) {
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShowerPlane::Clear for "
- << GetApparatus()->GetName() << ":" << endl;
+ << GetParent()->GetPrefix() << ":" << endl;
cout << " Cleared ADC hits for plane " << GetName() << endl;
cout << "---------------------------------------------------------------\n";
@@ -229,7 +228,7 @@ Int_t THcShowerPlane::Decode( const THaEvData& evdata )
if ( ((THcShower*) GetParent())->fdbg_decoded_cal ) {
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShowerPlane::Decode for "
- << GetApparatus()->GetName() << ":" << endl;
+ << GetParent()->GetPrefix() << ":" << endl;
cout << " Called for plane " << GetName() << endl;
cout << "---------------------------------------------------------------\n";
@@ -356,7 +355,7 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShowerPlane::ProcessHits for "
- << GetApparatus()->GetName() << ":" << endl;
+ << fParent->GetPrefix() << ":" << endl;
cout << " nrawhits = " << nrawhits << " nexthit = " << nexthit << endl;
cout << " Sparsified hits for HMS calorimeter plane #" << fLayerNum
@@ -442,7 +441,7 @@ Int_t THcShowerPlane::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShowerPlane::AcculatePedestals for "
- << GetApparatus()->GetName() << ":" << endl;
+ << GetParent()->GetPrefix() << ":" << endl;
cout << "Processed hit list for plane " << GetName() << ":\n";
@@ -498,7 +497,7 @@ void THcShowerPlane::CalculatePedestals( )
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShowerPlane::CalculatePedestals for"
- << GetApparatus()->GetName() << ":" << endl;
+ << GetParent()->GetPrefix() << ":" << endl;
cout << " ADC pedestals and thresholds for calorimeter plane "
<< GetName() << endl;