From 2680e82a126e7cf4ea43bf0ecf0d12f241d08b13 Mon Sep 17 00:00:00 2001
From: "Stephen A. Wood" <saw@jlab.org>
Date: Wed, 18 Feb 2015 14:32:11 -0500
Subject: [PATCH] Modify hcana for compatibility with analyzer(podd) 1.6
Modify THcFormula and THcParmList to match changes to THaFormula (needed
for reports) Small changes in THcInterface and detector classes for
compatibility with new OO decoder and podd 1.6
---
podd | 2 +-
src/THcAerogel.cxx | 1 +
src/THcCherenkov.cxx | 1 +
src/THcDC.cxx | 1 +
src/THcDriftChamber.cxx | 1 +
src/THcDriftChamberPlane.cxx | 4 +-
src/THcFormula.cxx | 285 +++++++++++++++++++++++++----------
src/THcFormula.h | 6 +-
src/THcInterface.cxx | 6 +-
src/THcParmList.cxx | 12 +-
src/THcShowerPlane.cxx | 13 +-
11 files changed, 234 insertions(+), 98 deletions(-)
diff --git a/podd b/podd
index cc0ebfe..1154ca9 160000
--- a/podd
+++ b/podd
@@ -1 +1 @@
-Subproject commit cc0ebfe4b74446e27dc5a7040dfe2b4606b24f9e
+Subproject commit 1154ca9087e06ede79833e8ce42cc04cbb99f479
diff --git a/src/THcAerogel.cxx b/src/THcAerogel.cxx
index da2716e..6106355 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 236db7c..31a3418 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 09ab985..f1e8290 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 61f5c02..b9ed691 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 5951641..3df5a7b 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 7740c06..114d8f6 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 0294ee6..20c962a 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 e7c1d5b..53c190e 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 e6caee6..23ff5be 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 8bd3d32..81c9f67 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;
--
GitLab