AutoPIDMaker.cxx

Go to the documentation of this file.

// $Id: AutoPIDMaker.cxx,v 1.11 2007/11/11 09:12:29 rhatcher Exp $
//
// FILL_IN: [Document your code!!]
//
// kordosky@hep.utexas.edu

#include <cmath>

#include "CalDetPID/AutoPIDMaker.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "RawData/RawRunCommentBlock.h"
#include "RawData/RawRecord.h"
#include "RawData/RawHeader.h"
#include "RawData/RawDaqHeader.h"
#include "Conventions/Munits.h"
#include "CalDetPID/CalDetParticleType.h"
#include "CalDetPID/PIDFormulas.h"
#include "CalDetPID/CalDetTOFRange.h"
#include "CalDetPID/CalDetCERRange.h"
#include "CalDetPID/CalDetCERTimeWin.h"
#include "CalDetPID/CalDetOverlapWin.h"
#include "CalDetPID/CalDetBeamMomentum.h"
#include "Validity/VldTimeStamp.h"
#include "Validity/VldContext.h"
#include "Validity/VldRange.h"

#include "DatabaseInterface/DbiCascader.h"
#include "DatabaseInterface/DbiTableProxyRegistry.h"
#include "DatabaseInterface/DbiWriter.h"
#include "DatabaseInterface/DbiResultPtr.h"
#include "DatabaseInterface/DbiValidityRec.h"


//#include <string>
//#include <map>
#include <sstream>
#include <cstdlib>
#include <cmath>
#include <vector>

using namespace std;
using namespace PIDFormulas;

JOBMODULE(AutoPIDMaker, "AutoPIDMaker",
          "Constructs caldet pid tables from info in the run start block");
CVSID("$Id: AutoPIDMaker.cxx,v 1.11 2007/11/11 09:12:29 rhatcher Exp $");
//......................................................................

//std::string AutoPIDMaker::ftestcomment="beam momentum (GeV) : 5.0\nUS Cerenkov pressure : 0.2\nDS Cerenkov pressure : 3.0";

std::string AutoPIDMaker::ftestcomment="beam momentum (GeV) : +5.0(geV)\nUS Cerenkov pressure : 0.2\nDS Cerenkov pressure : 3.0";

AutoPIDMaker::AutoPIDMaker() : 
     fdidrun(false), ftl(0), ftp(0), ftw(0),
     fbp(0), fuspr(0), fmidpr(0), fdspr(0)
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
     const Registry& r = DefaultConfig();
     Config(r);

}

//......................................................................

AutoPIDMaker::~AutoPIDMaker()
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
}

//......................................................................

void AutoPIDMaker::BeginJob()
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================


}

//......................................................................

void AutoPIDMaker::EndJob()
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
}

//......................................................................

void AutoPIDMaker::BeginFile()
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
}

//......................................................................

void AutoPIDMaker::EndFile()
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
}

//......................................................................

void AutoPIDMaker::BeginRun()
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================

     fdidrun=false;

}

//......................................................................

void AutoPIDMaker::EndRun()
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
}

//......................................................................

JobCResult AutoPIDMaker::Reco(MomNavigator* mom )
{
//======================================================================
//
// Calculate PID parameters and write to db
//
//======================================================================

     // if I've already seen this run then do nothing!
     if(fdidrun==true) return JobCResult::kPassed;
     
     // do a check on the data
     
     if(!SanityCheck()) {
          MSG("CalDetPID", Msg::kWarning)<<"SanityCheck() has failed. Will not calculate PID criteria for this run!"<<endl; 
          return JobCResult::kPassed;
     }

     // a list of particles and their masses, in gev
     const int npart = 5;
     CalDetParticleType::CalDetParticleType_t 
          particles[npart]={CalDetParticleType::kElectron,
                        CalDetParticleType::kMuon,
                        CalDetParticleType::kPion,
                        CalDetParticleType::kKaon,
                        CalDetParticleType::kProton};
     float m[npart] = {5E-4, 0.1057, 0.1396, 0.494, 0.938}; 

     // vectors to hold dbi table rows
     vector<CalDetTOFRange*> trv;
     vector<CalDetCERRange*> crv;
     vector<CalDetOverlapWin*> owv;
     vector<CalDetCERTimeWin*> ctv;

     // loop over particle types
     for(int i=0; i<npart; i++){

          // TOF CALCULATION
          const float tof_lim=4096.0; // should be settable?
          
          // ftp= Mean for electrons (by convention)
          const float tof_mean=TimeOfFlight(m[i],fbp,ftl) 
               + ftp-TimeOfFlight(m[0], fbp, ftl);
          float w = ftofsigma;
          float t = 0.0;
          if(particles[i]==CalDetParticleType::kProton && fcorrproton!=0){
             // correct tof for low energy protons
             w = CorrProtonWidth(fbp, w);
             t = CorrProtonTOF(fbp);
          }
          trv.push_back( new CalDetTOFRange( particles[i],
                                         -tof_lim, tof_lim,
                                        -tof_lim, tof_lim,
                                        -tof_lim, tof_lim,
                                        (tof_mean+t)-w*ftw, 
                                        (tof_mean+t)+w*ftw,
                                        -2*tof_lim, 2*tof_lim) );
     
          // CER CALCULATION
          const float n_minus_one_C02 = 4.1E-4;
          // dont care values will pass everything!
          const float dont_care_low=-1.0; 
          const float dont_care_high=1E5;
          // zero values will require the signal to be zero!
          const float zero_low=-1.0;
          const float zero_high=1.0;
          // cer timing window
          const float bigctw=1E5;

          //uscer
          float uslow=dont_care_low;
          float ushigh=dont_care_high;
          float usctwlow=-bigctw;
          float usctwhigh=bigctw;
          if(fuspr==-1.0){// do not use uscer
          }
          else if(CERThreshold(m[i],fbp,n_minus_one_C02)<fuspr){
               // this species over threshold for uscer
               uslow=fuscerlow;
               ushigh=fuscerhigh;
               usctwlow=fusctwlow;
               usctwhigh=fusctwhigh;
          }
          else {
               // this species under threshold for uscer
               // require no signals for this particle
               uslow=zero_low;
               ushigh=zero_high;
          }

          //midcer
          float midlow=dont_care_low;
          float midhigh=dont_care_high;
          float midctwlow=-bigctw;
          float midctwhigh=bigctw;
          if(fmidpr==-1.0){// do not use midcer
          }
          else if(CERThreshold(m[i],fbp,n_minus_one_C02)<fmidpr){
               // this species over threshold for midcer
               midlow=fmidcerlow;
               midhigh=fmidcerhigh;
               midctwlow=fmidctwlow;
               midctwhigh=fmidctwhigh;
          }
          else {
               // this species under threshold for midcer
               // require no signals for this particle
               midlow=zero_low;
               midhigh=zero_high;
          }

          //dscer
          float dslow=dont_care_low;
          float dshigh=dont_care_high;
          float dsctwlow=-bigctw;
          float dsctwhigh=bigctw;
          if(fdspr==-1.0){// do not use dscer
          }
          else if(CERThreshold(m[i],fbp,n_minus_one_C02)<fdspr){
               // this species over threshold for dscer
               dslow=fdscerlow;
               dshigh=fdscerhigh;
               dsctwlow=fdsctwlow;
               dsctwhigh=fdsctwhigh;
          }
          else {
               // this species under threshold for dscer
               // require no signals for this particle
               dslow=zero_low;
               dshigh=zero_high;
          }
          
          // store the cerrange table row
          crv.push_back( new CalDetCERRange(particles[i],
                                            midlow, midhigh,
                                            dslow, dshigh,
                                            uslow, ushigh) );
          
          // store cer timing window table row
          ctv.push_back( new CalDetCERTimeWin(particles[i],
                                              midctwlow, midctwhigh,
                                              dsctwlow, dsctwhigh,
                                              usctwlow, usctwhigh) );

          // construct overlap window
          float ollow=-500.0;
          float olhigh=500.0;
          if(particles[i]==CalDetParticleType::kElectron){
               ollow=felecollow;
               olhigh=felecolhigh;
          }
          else {
               ollow=fhadollow;
               olhigh=fhadolhigh;
          }
          owv.push_back( new CalDetOverlapWin(particles[i],ollow,olhigh) );
                         
     }

     // ok, now can write out to the db
     bool db_ok=true;
     // if we are testing, then require that we write to a temp table
     if(ftest!=0) fwritetemp=1;
     // the job can choose to write to a temp table
     if(fwritetemp!=0){
          MSG("CalDetPID", Msg::kDebug)<<"Creating temp tables."<<endl;
          // then we should write to a temp table
          DbiCascader& cascader 
               = DbiTableProxyRegistry::Instance().GetCascader();
          // beam momentum
          Int_t bp_no = 
               cascader.CreateTemporaryTable("CALDETBEAMMOMENTUM",
                                             CalDetBeamMomentum::GetTableDesc());
          if(bp_no<0){
               MSG("CalDetPID", Msg::kError)
                    <<"Trying to write a temp table but "
                    <<"no database will accept it!"<<endl;
               db_ok=false;
          }
          // tof range
          Int_t tr_no = 
               cascader.CreateTemporaryTable("CALDETTOFRANGE",
                                             CalDetTOFRange::GetTableDesc());
          if(tr_no<0){
               MSG("CalDetPID", Msg::kError)
                    <<"Trying to write a temp table but "
                    <<"no database will accept it!"<<endl;
               db_ok=false;
          }
          // cer range
          Int_t cr_no = 
               cascader.CreateTemporaryTable("CALDETCERRANGE",
                                             CalDetCERRange::GetTableDesc());
          if(cr_no<0){
               MSG("CalDetPID", Msg::kError)
                    <<"Trying to write a temp table but "
                    <<"no database will accept it!"<<endl;
               db_ok=false;
          }
          // cer win
          Int_t cw_no = 
               cascader.CreateTemporaryTable("CALDETCERTIMEWIN",
                                             CalDetCERTimeWin::GetTableDesc());
          if(cw_no<0){
               MSG("CalDetPID", Msg::kError)
                    <<"Trying to write a temp table but "
                    <<"no database will accept it!"<<endl;
               db_ok=false;
          }
          // overlap win
          Int_t ow_no = 
               cascader.CreateTemporaryTable("CALDETOVERLAPWIN",
                                             CalDetOverlapWin::GetTableDesc());
          if(ow_no<0){
               MSG("CalDetPID", Msg::kError)
                    <<"Trying to write a temp table but "
                    <<"no database will accept it!"<<endl;
               db_ok=false;
          }
          
     }
     // do the actual writing
     // get the run start time from the raw header
     MSG("CalDetPID", Msg::kDebug)<<"Record business."<<endl;
     VldTimeStamp start;
     Int_t run = 0;
     const RawRecord* rr = 
          dynamic_cast<const RawRecord*>(mom->GetFragment("RawRecord"));
     if(!rr){
          MSG("CalDetPID", Msg::kError)
               <<"Did not find a RawRecord. Will not write any tables."<<endl;
          db_ok=false;
     }
     else{
          const RawDaqHeader* rh = 
               dynamic_cast<const RawDaqHeader*>(rr->GetRawHeader());
          if(rh) {
               start=rh->GetVldContext().GetTimeStamp();
               run = rh->GetRun();
          }
          else{
               MSG("CalDetPID", Msg::kError)
                    <<"Did not find a RawHeader. "
                    <<"Will not write any tables."<<endl;
               db_ok=false;
          }
     }
     MSG("CalDetPID", Msg::kDebug)<<"DB ok. Trying to write."<<endl;
     if(db_ok){
          // if the db is ok at this point we can try to write
          const Detector::Detector_t det=Detector::kCalDet;
          const SimFlag::SimFlag_t sim=SimFlag::kData;
          const int aggno=-1;
          const Dbi::Task task=0;
          const VldTimeStamp create=start;
          // ending vld time = start + offset
          // the plan is subsequent runs within the offset period
          // will overlay and overide this one
          // Thanks dbi!
          const time_t end_offset_sec=12*3600; // 12 hour validity time
          VldTimeStamp end=start; end.Add(VldTimeStamp(end_offset_sec,0));
          const VldRange range(det,sim,start,end,"AutoPIDMaker");
          
          // beam momentum
          CalDetBeamMomentum cdbm(run, fbp);
          DbiWriter<CalDetBeamMomentum> w_bp(range,aggno,task,create);
          MSG("CalDetPID", Msg::kInfo)<<"Writing row:\n"
                                      << cdbm<<endl;

          w_bp<<cdbm;
          w_bp.Close();

          // the multiple row writes
          DbiWriter<CalDetTOFRange> w_tr(range,aggno,task,create);
          DbiWriter<CalDetCERRange> w_cr(range,aggno,task,create);
          DbiWriter<CalDetCERTimeWin> w_ct(range,aggno,task,create);
          DbiWriter<CalDetOverlapWin> w_ow(range,aggno,task,create);
          MSG("CalDetPID", Msg::kDebug)<<"Looping to write rows."<<endl;
          for(int i=0; i<npart; i++){
               // trv, crv, owv, ctv
               // tof range
               MSG("CalDetPID", Msg::kInfo)<<"Writing row:\n"
                                            << *(trv[i])<<endl;
               w_tr<<*(trv[i]);
               // cer range
               MSG("CalDetPID", Msg::kInfo)<<"Writing row:\n"
                                            << *(crv[i])<<endl;
               w_cr<<*(crv[i]);
               // cer window
               MSG("CalDetPID", Msg::kInfo)<<"Writing row:\n"
                                            << *(ctv[i])<<endl;
               w_ct<<*(ctv[i]);
               // overlap window
               MSG("CalDetPID", Msg::kInfo)<<"Writing row:\n"
                                            << *(owv[i])<<endl;
               w_ow<<*(owv[i]);
          }
          // close writers
          w_tr.Close();
          w_cr.Close();
          w_ct.Close();
          w_ow.Close();
          
     }

     MSG("CalDetPID", Msg::kDebug)<<"Now deleting table rows."<<endl;
     // delete table row entries from vectors
     for(int i=0; i<npart; i++){
          if(trv[i]!=0) {delete trv[i]; trv[i]=0;}
          if(crv[i]!=0) {delete crv[i]; crv[i]=0;}
          if(ctv[i]!=0) {delete ctv[i]; ctv[i]=0;}
          if(owv[i]!=0) {delete owv[i]; owv[i]=0;}
     }

     // all done
     fdidrun=true;

     return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

//......................................................................

JobCResult AutoPIDMaker::Get(MomNavigator* mom)
{
//======================================================================
// Looks for a run comment block
// Reads the block and extracts member variables
//======================================================================
     
     // if I've already seen this run then do nothing!
     if(fdidrun==true) return JobCResult::kPassed;
  
     const RawRecord* rr = 
          dynamic_cast<const RawRecord*>(mom->GetFragment("RawRecord"));
     
     if(!rr){
          MSG("CalDetPID", Msg::kInfo)
               <<"Did not find a RawRecord. Doing nothing."<<endl;
          return JobCResult::kFailed;
     }
     
     const RawRunCommentBlock* rcb = 
          dynamic_cast<const RawRunCommentBlock*>(
               rr->FindRawBlock("RawRunCommentBlock"));
     if(!rcb){
          MSG("CalDetPID", Msg::kDebug)
               <<"Could not find a RawRunCommentBlock!\n"
               <<"Cannot compute PID parameters until I do!"
               <<endl;
          return JobCResult::kFailed;
     }
     
     // got a good runcommentblock
     MSG("CalDetPID", Msg::kInfo)
          <<"Got the following RawRunCommentBlock:\n"<<*rcb<<endl;
     string s;
     if(ftest==0){
          // copy run comment into string
          s=rcb->GetRunComment();
          MSG("CalDetPID", Msg::kInfo)<<"The comment was: \n"<<s<<endl;
     }
     else {
          s=ftestcomment;
          MSG("CalDetPID", Msg::kInfo)<<"Testing with the comment:\n"
                                      <<s<<endl;
     }
     // convert run comment into a map
     const int slen = s.length();
     // if the run comment doesn't end with a newline then add one
     if(s[slen-1]!='\n') s+='\n';
     map<string,string> ms;
     RCtoMap(s,ms);

     // dig through the map to pull out relevant values
     // try to reject crappy entries
     
     bool ok=true;
     double dval=0.0;
     
     if(FindDouble(ms,"beam momentum (GeV) ", dval)) fbp=dval;
     else ok&=false;

     if(FindDouble(ms,"US Cerenkov pressure ",dval)) fuspr=dval;
     else ok&=false;

     if(FindDouble(ms,"DS Cerenkov pressure ",dval)) fdspr=dval;
     else ok&=false;
     
//     if(FindDouble(ms,"MIDCER pressure (atm) ",dval)) fmidpr=dval;
//     else ok&=false;

     // all done!
//     fdidrun=true;
     return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

//......................................................................

JobCResult AutoPIDMaker::Ana(const MomNavigator* /* mom */)
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

//......................................................................

const Registry& AutoPIDMaker::DefaultConfig() const
{
//======================================================================
// Supply the default configuration for the module
//======================================================================
  static Registry r; // Default configuration for module

  // Set name of config
  std::string name = this->GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  // Set values in configuration
  r.UnLockValues();
  r.Set("writetemp", 1);
  r.Set("toflength", 9.11);
  r.Set("tofped",    417.2);
  r.Set("tofwidth",  3.0);
  r.Set("beamp",     0.0);
  r.Set("usceratm",  -1.0);
  r.Set("midceratm", -1.0);
  r.Set("dsceratm",  -1.0);

  // tof sigma
  r.Set("tofsigma", 2.5);

  // cer selection ranges
  r.Set("uscerlow", 50.0);
  r.Set("uscerhigh", 1.6E4);
  r.Set("midcerlow", 50.0);
  r.Set("midcerhigh", 1.6E4);
  r.Set("dscerlow", 50.0);
  r.Set("dscerhigh", 1.6E4);

  // cer timing window
  r.Set("usctwhigh", -63.0);
  r.Set("usctwlow", -103.0);
  r.Set("midctwhigh", -1E5);
  r.Set("midctwlow", 1E5);
  r.Set("dsctwhigh", 14.0);
  r.Set("dsctwlow", -26.0);

  // overlap ranges
  r.Set("elecollow", -80.0);
  r.Set("elecolhigh", 0.0);
  r.Set("hadollow", -80.0);
  r.Set("hadolhigh", 160.0);

  // for testing
  r.Set("test", 0);

  // low energy proton tof correction
  r.Set("corrproton", 1);

//  r.Set("testcomment","");
  r.LockValues();

  return r;
}

//......................................................................

void AutoPIDMaker::Config(const Registry& r)
{
//======================================================================
// Configure the module given the Registry r
//======================================================================
  double tmpd;
  int tmpi;

  if (r.Get("writetemp",tmpi)) { fwritetemp = tmpi; }

  if (r.Get("toflength",tmpd)) { ftl = tmpd; }
  if (r.Get("tofped",tmpd)) { ftp = tmpd; }
  if (r.Get("tofwidth",tmpd)) { ftw = tmpd; }
  if (r.Get("beamp",tmpd)) { fbp = tmpd; }
  if (r.Get("usceratm",tmpd)) { fuspr = tmpd; }
  if (r.Get("midceratm",tmpd)) { fmidpr = tmpd; }
  if (r.Get("dsceratm",tmpd)) { fdspr = tmpd; }

  // tof sigma
  if (r.Get("tofsigma",tmpd)) { ftofsigma = tmpd; }

  // cer selection ranges
  if (r.Get("uscerlow",tmpd)) { fuscerlow = tmpd; }
  if (r.Get("uscerhigh",tmpd)) { fuscerhigh = tmpd; }
  if (r.Get("midcerlow",tmpd)) { fmidcerlow = tmpd; }
  if (r.Get("midcerhigh",tmpd)) { fmidcerhigh = tmpd; }
  if (r.Get("dscerlow",tmpd)) { fdscerlow = tmpd; }
  if (r.Get("dscerhigh",tmpd)) { fdscerhigh = tmpd; }

  // cer timing window
  if (r.Get("usctwlow",tmpd)) { fusctwlow = tmpd; }
  if (r.Get("usctwhigh",tmpd)) { fusctwhigh = tmpd; }
  if (r.Get("midctwlow",tmpd)) { fmidctwlow = tmpd; }
  if (r.Get("midctwhigh",tmpd)) { fmidctwhigh = tmpd; }
  if (r.Get("dsctwlow",tmpd)) { fdsctwlow = tmpd; }
  if (r.Get("dsctwhigh",tmpd)) { fdsctwhigh = tmpd; }

  // overlap window
  if (r.Get("elecollow",tmpd)) { felecollow = tmpd; }
  if (r.Get("elecolhigh",tmpd)) { felecolhigh = tmpd; }
  if (r.Get("hadollow",tmpd)) { fhadollow = tmpd; }
  if (r.Get("hadolhigh",tmpd)) { fhadolhigh = tmpd; }

  // for testing
  if (r.Get("test", tmpi)) { ftest=tmpi;}

  // low energy proton tof correction
  if (r.Get("corrproton", tmpi)) {fcorrproton=tmpi;}

  

/*
  const char* tmpc;
  if (r.Get("testcomment", tmpc)) {
       MSG("CalDetPID", Msg::kDebug)<<"Got testcomment: \n"<<tmpc<<endl;
       ftestcomment=tmpc;
  }
*/

}

//......................................................................

void AutoPIDMaker::Help()
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
}


void AutoPIDMaker::RCtoMap(std::string& s, 
                           std::map<std::string, std::string>& ms) const
{
     // converts the run comment string into a map
     // the key is the bit left of the : and the value is the bit to the right
     std::istringstream ss(s);
     while(!ss.eof()){
          char lbuf[256]={'0'};
          char rbuf[256]={'0'};
          ss.getline(lbuf, 256, ':');
          ss.getline(rbuf, 256);
          if(ss.eof()) break;
          ms.insert(pair<std::string,std::string>(std::string(lbuf), 
                                                  std::string(rbuf)));
     }
     return;
}

bool AutoPIDMaker::ExtractDouble(std::string& v, double& val) const
{
   // this algorithm was bad when people typed "+1 GEV" or some such
   // try to fix it by stripping out anything that is not a number
   // or .+-
   int found_dig=0;
   string vstrip;
   for(unsigned int i=0; i<v.length(); i++){
      if( (isdigit(v[i]))||(v[i]=='.')||(v[i]=='+')||(v[i]=='-') ){
         // accept these
         vstrip+=v[i];
         if(isdigit(v[i])) found_dig++;
      } 

   }
   // checks if s can be converted to a number and extracts it if so

   // the checking bit
   if(vstrip.length()==0) return false;// nothing at all
   if(found_dig==0) return false; // no numbers
   
/*
   bool ok=true;
   for(unsigned int i=0; i<vstrip.length(); i++){
      if(!( (isdigit(vstrip[i]))
            ||(isspace(vstrip[i]))
            ||(vstrip[i]=='.')
            ||(vstrip[i]=='+')
            ||(vstrip[i]=='-')
            )){ok=false; break;}
   }
     
   if(!ok) return false;
*/

   // looks like a number so extract it
   istringstream in(vstrip);
   in>>val;

   return true;
}

bool AutoPIDMaker::FindDouble(std::map<std::string, std::string>& ms,
                              const char* k, double& value) const 
{

     double dval;
     typedef map<string,string>::iterator MIT;
     MIT it=ms.begin();
     string key(k);
     if(key.length()==0) return false;
     
     it = ms.find(key);
     if(it!=ms.end()){ 
          if(ExtractDouble(it->second, dval)){
               value=dval;
               return true;
          }
          else{
               MSG("CalDetPID", Msg::kError)
                    <<"Failed to extract a double from \""<<it->second
                    <<"\"!"<<endl;
               return false;
          }
     }
     else{
          MSG("CalDetPID", Msg::kError)
               <<"Could not find a value for the key \""<<key<<"\" !"<<endl;
          return false;
     }
}

bool AutoPIDMaker::SanityCheck() const
{
     // does a simple sanity check on the member data to decide if 
     // it is reasonable to try and make pid tables from it
     bool ok=true;
     // check for reasonable beam momentum
     if( !(abs(fbp)>0.05 && abs(fbp)<11.0) ) ok&=false;
     // check tof baseline
     if( !(ftl>0.0 && ftl<13.0) ) ok&=false;
     // check tof width... just needs to be > 0
     if( !(ftw>0.0) ) ok&=false;
     // check tof pedestal
     if( !(abs(ftp)<4096) ) ok&=false;
     // check if cer pressure is in the reasonable range
     // or, if pressure==-1 then that channel will be considered unused
     // it will not be used in selections and any signal will be valid
     // for any type of particle
     if( !((fuspr>0.0 && fuspr<5.0) || fuspr==-1.0) ) ok&=false;
     if( !((fmidpr>0.0 && fmidpr<5.0) || fmidpr==-1.0) ) ok&=false;
     if( !((fdspr>0.0 && fdspr<5.0) || fdspr==-1.0) ) ok&=false;

     // all done
     return ok;
}

float AutoPIDMaker::CorrProtonWidth(float p, float w) const
{
   if(p<0.8) p=0.8; // apply correction for 0.8 below 0.8
   float corw = 2.91 + 13.56/pow(p,4);
   if(corw<w) corw=w; // use the width w as a minimum
   return corw;
}

float AutoPIDMaker::CorrProtonTOF(float p) const
{
   float tc=0;
   if(! (fabs(p)>1.2 ) ) tc = -3.6 + 11.15/pow(p,4);
   return tc;
}

---