AlgCalDetPID.cxx

Go to the documentation of this file.

// $Id: AlgCalDetPID.cxx,v 1.12 2005/02/14 22:24:43 kordosky Exp $
//
// AlgCalDetPID.cxx
//
// AlgCalDetPID is a concrete algorithm class.
//
// Author:  Mike Kordosky 2003.04.27

#include <cassert>
#include "TAxis.h"
#include "CalDetPID/AlgCalDetPID.h"

#include "CandData/CandRecord.h"
#include "Candidate/CandContext.h"
#include "MessageService/MsgService.h"
#include "Validity/VldContext.h"
#include "CalDetPID/CandCalDetPIDHandle.h"
#include "DatabaseInterface/DbiResultPtr.h"
#include "CalDetPID/CalDetTOFRange.h"
#include "CalDetPID/CalDetCERRange.h"
#include "CalDetPID/CalDetOverlapWin.h"
#include "CalDetPID/CalDetCERTimeWin.h"
#include "CandDigit/CandDigitHandle.h"
#include "Conventions/Munits.h"
#include "Conventions/DetectorType.h"
#include "Conventions/SimFlag.h"
#include "Conventions/PlaneView.h"
#include "Conventions/ElecType.h"
#include "Algorithm/AlgConfig.h"
#include "Record/SimSnarlRecord.h"
#include "MinosObjectMap/MomNavigator.h"
#include "Plex/PlexHandle.h"
#include "Conventions/ReadoutType.h"

#include "TClonesArray.h"
#include "TParticle.h"

ClassImp(AlgCalDetPID)

//______________________________________________________________________
CVSID("$Id: AlgCalDetPID.cxx,v 1.12 2005/02/14 22:24:43 kordosky Exp $");

//______________________________________________________________________
AlgCalDetPID::AlgCalDetPID():
   fFirstEvent(kTRUE),
   fNAllowedBefore(1), 
   fNAllowedAfter(1),
   fNAllowedOut(2), 
   fNRequiredIn(1), 
   fTimeFileNum(0)
{
   fTimeFile=0;
   fThist[0]=0;
   fThist[1]=0;

}

//______________________________________________________________________
AlgCalDetPID::~AlgCalDetPID()
{
   if(fTimeFile){
      if(fTimeFile->IsOpen()){
         fTimeFile->Close();
         delete fTimeFile;
         fTimeFile = 0;
      }
   }

}

//______________________________________________________________________
void AlgCalDetPID::RunAlg(AlgConfig &  ac , 
                          CandHandle &ch, CandContext &cx)
{
     bool ok=true;

   MSG("CalDetPID", Msg::kDebug) 
        << "Starting AlgCalDetPID::RunAlg()" << endl;

   if(fFirstEvent){
      //get filename from alg config
      Int_t tfn;
      if(ac.Get("TimeFileNum", tfn)){
         fTimeFileNum = tfn;
         MSG("CalDetPID", Msg::kDebug)<<"got TimeFileNum from config "
                                      <<fTimeFileNum<<endl;
      }
      else{
         fTimeFileNum = 0;
         MSG("CalDetPID", Msg::kError)<<"Could not get TimeFileName from config."
                                      <<"  Using default "<<fTimeFileNum
                                      <<" instead"<<endl;
      }
      string fTimeFileName;
      switch(fTimeFileNum){
      case 0: //far-far readout t11 2002
         fTimeFileName = "timinghists-2002-t11.root";
         break;
      case 1: //far-far readout t7 2002
         fTimeFileName = "timinghists-2002-t7.root";
         break;
      case 2: //far-near readout 2003
         fTimeFileName = "timinghists-2003-t7-nf.root";
         break;
      case 3: //near readout t7 2003
         fTimeFileName = "timinghists-2003-t7-n.root";
         break;
      case 4: //near readout t11 2003
         fTimeFileName = "timinghists-2003-t11.root";
         break;
      default:
         fTimeFileName = "timinghists.root";
         break;
      }
      
      // open input root file
      string fname(getenv("SRT_PRIVATE_CONTEXT"));
      fname+=("/CalDetPID/data/"+fTimeFileName);
      
      fTimeFile = new TFile(fname.c_str());
      if(!fTimeFile->IsOpen()){
         MSG("CalDetPID",Msg::kInfo)<<"Could not find "<<fname<<" in the base release. Will try the test release."<<endl;
         string fname_save=fname;
         fname=getenv("SRT_PUBLIC_CONTEXT");
         fname+="/CalDetPID/data/"+fTimeFileName;
         //crashes if it get's here because now fTimeFile goes out of scope
         //tv 6-16-2004
//       TFile *fTimeFile = new TFile(fname.c_str());
         fTimeFile = new TFile(fname.c_str());
         if(!fTimeFile->IsOpen()){
            MSG("CalDetPID",Msg::kError)<<"Could not open timing hist file"
                                        <<fname<<endl
                                        <<"Also tried "<<fname_save<<endl;
            return;
         }
      }
      MSG("CalDetPID",Msg::kInfo)<<"Opened timing hist file "<<fname<<endl;
      //      fTimeFile->Print();
      for(int i=0;i<2;i++){
         char hn[50];
         sprintf(hn,"all_thist_c%d",i);
         fThist[i] = (TH1F *)fTimeFile->Get(hn);
         if(fThist[i]!=0){
            MSG("CalDetPID",Msg::kDebug)<<"Found histogram "<<hn<<endl;
         }
         else{
            MSG("CalDetPID",Msg::kError)<<"Could not find histogram "<<hn<<endl;
         }
      }
      fFirstEvent=kFALSE;
   }


   assert(ch.InheritsFrom("CandCalDetPIDHandle"));
   CandCalDetPIDHandle& cpidh = static_cast<CandCalDetPIDHandle&>(ch);

   assert(cx.GetDataIn());
   const CandRecord* crec = cx.GetCandRecord();
   // Get VldContext
   const VldContext vc = *(crec->GetVldContext());
   // is this data?
   const bool isdata = (vc.GetSimFlag()&SimFlag::kData) ? true : false;
   // is this mc?
   const bool ismc = ((vc.GetSimFlag()&SimFlag::kMC) ||
                      (vc.GetSimFlag()&SimFlag::kReroot)) ? true : false;
   
   if(isdata){
        // for data we must have a caldetsi candidate
        const CandCalDetSIHandle *csih =
             dynamic_cast<const CandCalDetSIHandle *> (cx.GetDataIn());
        assert(csih);
        // we must also have some digits
        const CandDigitListHandle* cdlh=
             dynamic_cast<const CandDigitListHandle*>
             (crec->FindCandHandle("CandDigitListHandle","canddigitlist"));
        if(cdlh == 0) {//failure
             MSG("CalDetPID",Msg::kFatal)<<"cdlh==0"<<endl;
             assert(0);
        }

        const Double_t toftime=
             (csih->GetTofTimeStamp()*1.5625 - cdlh->GetAbsTime()/Munits::ns);

        // talk to the db
   
        DbiResultPtr<CalDetTOFRange> rp_tof(vc,0);
        if(!rp_tof.GetValidityRec()){
             MSG("CalDetPID", Msg::kError)
                  <<"Could not look up selection ranges "
                  <<"from CalDetTOFRange table!\n"
                  <<"Bad DbiResultPtr!?"
                  <<endl;
             ok=false;
        }

        DbiResultPtr<CalDetCERRange> rp_cer(vc,0);
        if(!rp_cer.GetValidityRec()){
             MSG("CalDetPID", Msg::kError)
             <<"Could not look up selection ranges "
             <<"from CalDetCERRange table!\n"
             <<"Bad DbiResultPtr!?"
             <<endl;
             ok=false;
        }


   
        // The algorithm:
        // Look through each entry in the db.  
        // There should not be many of them.
        // One entry may correspond to one or more types of particles.
        // The PARTICLETYPE entry in the table rows will be used to set high
        // the corresponding bits in case the data is in the appropriate range.
        // This can be done seperately for tof and cer with the results anded
        // at the end.

        // do id by tof
        UInt_t tof_pid=0;       
        if(ok){   
             // here is the tof data
             const Int_t tdc0=csih->GetTofTDC0();
             const Int_t tdc1=csih->GetTofTDC1();
             const Int_t tdc2=csih->GetTofTDC2();
             MSG("CalDetPID", Msg::kDebug)
                  <<"TOF 0 1 2 = "<<tdc0<<" "<<tdc1<<" "<<tdc2<<endl;
             // loop over rows
             const int nrows_tof = rp_tof.GetNumRows();
             for(int i=0; i<nrows_tof; i++){
                  // get a row
                  const CalDetTOFRange* tr = rp_tof.GetRow(i);
                  if(!tr){
                       // something is really hosed
                       MSG("CalDetPID", Msg::kFatal)
                            <<"Could not look up row "<<i
                            <<" in CalDetTOFRange table!"
                            <<endl;
                       ok=false;
                       assert(0);
                  }
                  else{
                       // for now, only test on tdc2 and tdc0
                       // tdc1 was never used...
                       // when we ran middle and ds counters the middle
                       // was plugged into tdc2
                       if((tdc0 < tr->GetTDC0High()) 
                          && (tdc0 > tr->GetTDC0Low())
                          && (tdc2 < tr->GetTDC2High()) 
                          && (tdc2 > tr->GetTDC2Low())
                          && (tdc2-tdc0 < tr->GetTDC2MinusTDC0High())
                          && (tdc2-tdc0 > tr->GetTDC2MinusTDC0Low())
                            )
                            tof_pid|=(tr->GetParticleType());
                  }
             }
        }

        // do id by cer
        UInt_t cer_pid=0;
        if(ok){
             // here is the cer data
             const Int_t cer0=csih->GetKovADC2();
             const Int_t cer1=csih->GetKovADC1();
             const Int_t cer2=csih->GetKovADC3();
             MSG("CalDetPID", Msg::kDebug)
                  <<"CER 0 1 2 = "<<cer0<<" "<<cer1<<" "<<cer2<<endl;
             // loop over rows
             const int nrows_cer = rp_cer.GetNumRows();
             for(int i=0; i<nrows_cer; i++){
                  const CalDetCERRange* cr = rp_cer.GetRow(i);
                  if(!cr){
                       // something is really hosed
                       MSG("CalDetPID", Msg::kFatal)
                            <<"Could not look up row "<<i
                            <<" in  CalDetCERRange table!"
                            <<endl;
                       ok=false;
                       assert(0);
                  }
                  else{  
                       if((cer0<cr->GetCER0High())
                          &&(cer0>cr->GetCER0Low())
                          &&(cer1<cr->GetCER1High())
                          &&(cer1>cr->GetCER1Low())
                          &&(cer2<cr->GetCER2High())
                          &&(cer2>cr->GetCER2Low())
                            ) 
                            cer_pid|=(cr->GetParticleType());
                  }
             }    
        }

        MSG("CalDetPID", Msg::kDebug)
             <<"tof_pid = "<<hex<<tof_pid<<dec
             <<",  Identified via TOF as "
             <<CalDetParticleType::AsString(tof_pid)<<endl;
        
        MSG("CalDetPID", Msg::kDebug)
             <<"cer_pid = "<<hex<<cer_pid<<dec
             <<",  Identified via CER as "
             <<CalDetParticleType::AsString(cer_pid)<<endl;
        
        // Now we have the basic pid information.
        // Look to see if this is an overlapper.

        // For things that are not overlappers, set the appropriate bits high
        // otherwise do not set them
        // If we cannot talk to the db then cannot_run_overlap=true;
        UInt_t overlap_pid=0;
        bool cannot_run_overlap=false;
        Float_t chi2 = -1.;
        DbiResultPtr<CalDetOverlapWin> rp_cow(vc,0);
        if(!rp_cow.GetValidityRec()){
             MSG("CalDetPID", Msg::kWarning)
                  <<"Could not look up selection ranges "
                  <<"from CalDetOverlapWin table!\n"
                  <<"Bad DbiResultPtr!?"
                  <<endl;
             // but maybe it could live with this?
             cannot_run_overlap=true;
             ok=false;
        }
        else{
             const int nrows_cow=rp_cow.GetNumRows();
             for(int i=0; i<nrows_cow; i++){
                  const CalDetOverlapWin* cow = rp_cow.GetRow(i);
                  if(!cow){
                       // something is really hosed
                       MSG("CalDetPID", Msg::kFatal)
                            <<"Could not look up row "<<i
                            <<" in  CalDetOverlapWin table!"
                            <<endl;
                       ok=false;
                       assert(0);
                  }
                  else{
                       const Float_t lowwin = cow->GetWinLow();
                       const Float_t highwin = cow->GetWinHigh();
                       Int_t nbefore,nafter,nout,nin;
                       HitTimeAnal(cdlh,toftime,lowwin,highwin,
                                   nbefore,nafter,nout,nin,chi2);
                       bool overlap = IsAnOverlap(ac,nbefore,nafter,
                                                  nout,nin);
                       if(!overlap){
                            overlap_pid|=(cow->GetParticleType());
                       }
                  }
             }
        }
                
        //
        // now work on certimewin
        UInt_t ctw_pid=0;
        bool cannot_run_ctw=false;
        DbiResultPtr<CalDetCERTimeWin> rp_ctw(vc,0);
        if(!rp_ctw.GetValidityRec()){
             MSG("CalDetPID", Msg::kError)
                  <<"Could not look up selection ranges "
                  <<"from CalDetCERTimeWin table!\n"
                  <<"Bad DbiResultPtr!?"
                  <<endl;
             // but maybe it could live with this?
             cannot_run_ctw=true;
             ok=false;
        }
        else{
             const int nrows_ctw=rp_ctw.GetNumRows();
             for(int i=0; i<nrows_ctw; i++){
                  const CalDetCERTimeWin* ctw = rp_ctw.GetRow(i);
                  if(!ctw){
                       // something is really hosed
                       MSG("CalDetPID", Msg::kFatal)
                            <<"Could not look up row "<<i
                            <<" in  CalDetCERTimeWin table!"
                            <<endl;
                       ok=false;
                       assert(0);
                  }
                  else{
                       const Float_t win0low = ctw->GetWin0Low();
                       const Float_t win0high = ctw->GetWin0High();
                       const Float_t win1low = ctw->GetWin1Low();
                       const Float_t win1high = ctw->GetWin1High();
                       const Float_t win2low = ctw->GetWin2Low();
                       const Float_t win2high = ctw->GetWin2High();
                       bool cer0=false;
                       bool cer1=false;
                       bool cer2=false;
                       MSG("CalDetPID", Msg::kDebug)
                           <<(*ctw)<<endl;
                       const bool in_win = CERInTimeWin(ac, csih,
                                                        cdlh->GetAbsTime(),
                                                        toftime,
                                                        win0low,win0high,
                                                        win1low,win1high,
                                                        win2low,win2high,
                                                        cer0, cer1, cer2);
                       if(in_win){
                            ctw_pid|=(ctw->GetParticleType());
                       }
                  }
             }
        }
        
        // final step:  
        // and the tof and cer results together
        const UInt_t pid_result = (tof_pid&cer_pid);
        MSG("CalDetPID", Msg::kDebug)
             <<"result = "<<hex<<pid_result<<dec
             <<",  Identified as "
             <<CalDetParticleType::AsString(pid_result)<<endl;
        cpidh.SetPIDType(pid_result);

        MSG("CalDetPID", Msg::kDebug)
             <<"overlap_pid = "<<hex<<overlap_pid<<dec
             <<",  Non overlapping for "
             <<CalDetParticleType::AsString(overlap_pid)<<endl;
        
        cpidh.SetNoOverlapBits(overlap_pid);

        MSG("CalDetPID", Msg::kDebug)
             <<"ctw_pid = "<<hex<<ctw_pid<<dec
             <<",  Cer in time window for "
             <<CalDetParticleType::AsString(ctw_pid)<<endl;

        cpidh.SetInCERTimeBits(ctw_pid);

        Bool_t nooverlap;
        if((pid_result&overlap_pid)==pid_result) nooverlap=kTRUE;
        else nooverlap=kFALSE;
        MSG("CalDetPID", Msg::kDebug)
             <<"NoOverlap = "<<( (nooverlap) ? "true": "false")<<endl;
        cpidh.SetNoOverlap(nooverlap);

        Bool_t incertime;
        if((pid_result&ctw_pid)==pid_result) incertime=kTRUE;
        else incertime=kFALSE;
        MSG("CalDetPID", Msg::kDebug)
             <<"InCERTime = "<<( (incertime) ? "true": "false")<<endl;
        cpidh.SetInCERTime(incertime);

        cpidh.SetOLChi2(chi2);
   }
   else if(ismc){
        cpidh.SetPIDType(CalDetParticleType::kUnknown);
        cpidh.SetNoOverlapBits(0);
        cpidh.SetInCERTimeBits(0);
        cpidh.SetNoOverlap(kTRUE);
        cpidh.SetInCERTime(kTRUE);
        cpidh.SetOLChi2(0.);
        // look for the primary in stdhep
        const SimSnarlRecord* simrec = dynamic_cast<const SimSnarlRecord*>
            (cx.GetMom()->GetFragment("SimSnarlRecord"));
        if(simrec){
            const TClonesArray* stdhep = dynamic_cast<const TClonesArray*> 
                (simrec->FindComponent("TClonesArray", "StdHep"));
            if(stdhep){
                // take the first in the list
                const TParticle* part = dynamic_cast<const TParticle*>
                    (stdhep->At(0));
                if(part){
                    Int_t ipdg = part->GetPdgCode();
                    if(abs(ipdg)==11) 
                        cpidh.SetPIDType(CalDetParticleType::kElectron);
                    else if(abs(ipdg)==13) 
                        cpidh.SetPIDType(CalDetParticleType::kMuon);
                    else if(abs(ipdg)==211) 
                        cpidh.SetPIDType(CalDetParticleType::kPion);
                    else if(abs(ipdg)==321) 
                        cpidh.SetPIDType(CalDetParticleType::kKaon);
                    else if(abs(ipdg)==2212) 
                        cpidh.SetPIDType(CalDetParticleType::kProton);
                    }
                }
            }
       }


}

//______________________________________________________________________
void AlgCalDetPID::Trace(const char * /* c */) const
{
}

bool AlgCalDetPID::HitTimeAnal(const CandDigitListHandle *cdlh, 
                               Double_t toftime, 
                               Float_t lowwin, Float_t highwin,
                               Int_t& nbefore, Int_t& nafter,
                               Int_t& nout, Int_t& nin, Float_t& chi2) const
{
          // returns false if there is an error
          
          // zero counters
          nbefore=0; nafter=0; nout=0; nin=0;
          
          //make histograms to hold timing of current event
          int NBFAR = fThist[0]->GetNbinsX();
          Axis_t LOWFAR = fThist[0]->GetBinLowEdge(1);
          Axis_t HIFAR = fThist[0]->GetBinLowEdge(NBFAR)+
             fThist[0]->GetBinWidth(NBFAR);

          int NBNEAR = fThist[1]->GetNbinsX();
          Axis_t LOWNEAR = fThist[1]->GetBinLowEdge(1);
          Axis_t HINEAR = fThist[1]->GetBinLowEdge(NBNEAR)+
             fThist[1]->GetBinWidth(NBNEAR);

          MSG("CalDetPID",Msg::kDebug)<<"making event thists "
                                      <<" NBFAR "<<NBFAR
                                      <<" LOWFAR "<<LOWFAR
                                      <<" HIFAR "<<HIFAR
                                      <<" NBNEAR "<<NBNEAR
                                      <<" LOWNEAR "<<LOWNEAR
                                      <<" HINEAR "<<HINEAR<<endl;

          TH1F ehfar("ehfar","ehfar",NBFAR,LOWFAR,HIFAR);
          TH1F ehnear("ehnear","ehnear",NBNEAR,LOWNEAR,HINEAR);
  
          // check digit list
          CandDigitHandleItr hi(cdlh->GetDaughterIterator());
          if(!hi.IsValid()){
               MSG("CalDetPID",Msg::kError)<<"DaugherIter not valid"<<endl;
               return false;    
          }
          //begin looping over hits
          for(hi.Reset();hi.IsValid();hi.Next()){
              // to check that it isn't a cosmic counter or unknown
              const PlaneView::PlaneView_t pv=
                   (*hi)->GetPlexSEIdAltL().GetPlaneView();
              // to check that it IS connected to scint
              const PlexHandle ph(*((*hi)->GetVldContext()));
              const ReadoutType::Readout_t rt= 
                   ph.GetReadoutType((*hi)->GetChannelId());
              const ElecType::Elec_t et = (*hi)->GetChannelId().GetElecType();
/*
  // the logic used to be
              if((pv==PlaneView::kA)||(pv==PlaneView::kB)
                 ||(pv==PlaneView::kUnknown)) {
*/
/*  changing the logic so that it will consider near detector hits
              if( !(pv==PlaneView::kU || pv==PlaneView::kV) 
                  || rt!=ReadoutType::kScintStrip || (et!=ElecType::kVA)){
                  // don't do anything with this hit
                   continue;
              }
*/
              if( !(pv==PlaneView::kU || pv==PlaneView::kV) 
                  || rt!=ReadoutType::kScintStrip){
                  // don't do anything with this hit
                   continue;
              }

              Float_t time = (*hi)->GetSubtractedTime(CalTimeType::kT0);
              time/=Munits::nanosecond;
              MSG("CalDetPID",Msg::kDebug)<<"hittime - toftime= "
                                             <<time<<"-"<<toftime
                                             <<" = "<<time-toftime<<endl;
              time-=toftime;
              if(et==ElecType::kVA){
                 ehfar.Fill(time);
                 if(time<lowwin){nbefore++; nout++;}
                 else if(time>highwin){nafter++; nout++;}
                 else nin++;
              }
              else if(et==ElecType::kQIE){
                 ehnear.Fill(time);
              }
          }
          chi2 = OverlapLikelihood(&ehfar,&ehnear);
          return true;
}

bool AlgCalDetPID::IsAnOverlap(AlgConfig &ac,
                               Int_t nbefore, Int_t nafter,
                               Int_t nout, Int_t nin)
{
     // returns false if there is no overlap
     // returns true if there is an overlap
     // get the number of hits allowed out of the time interval from
     // AlgConfig
     
     Int_t tmpi;
     if(ac.Get("NAllowedBefore", tmpi)) fNAllowedBefore=tmpi;
     if(ac.Get("NAllowedAfter", tmpi)) fNAllowedAfter=tmpi;
     if(ac.Get("NAllowedOut", tmpi)) fNAllowedOut=tmpi;
     if(ac.Get("NRequiredIn", tmpi)) fNRequiredIn=tmpi;
   
     if(nout<=fNAllowedOut&&
        nbefore<=fNAllowedBefore&&
        nafter<=fNAllowedAfter&&
        nin>=fNRequiredIn){
          return false;
     }
     else{
          return true;
     }  

}

bool AlgCalDetPID::CERInTimeWin(AlgConfig& /* ac */,
                                const CandCalDetSIHandle* csih,
                                Double_t abstime, Double_t toftime,
                                Float_t win0low, Float_t win0high,
                                Float_t win1low, Float_t win1high,
                                Float_t win2low, Float_t win2high,
                                bool& cer0, bool& cer1, bool& cer2){
     // returns in cerX true if counter X passes the selection
     // returns cer0&cer1&cer2

     cer0=cer1=cer2=false;
     if(csih==0){
          MSG("CalDetPID",Msg::kError)
               <<"CERInTimeWin got a null CandCalDetSIHandle!"
               <<"Returning false!"<<endl;
          return false;
     }
     else{
          // in principle could use algorithm config to avoid considering
          // one of the counters for some particles

          // candcaldetsihandle starts counting cer counters at 1
          const Double_t tdc_convert=1.5625;


          Int_t cer0adc=csih->GetKovADC2(); // yes, I really mean this
          Float_t cer0time= toftime-
              (csih->GetKovTimeStamp2()*tdc_convert -abstime/Munits::ns); // too
          Int_t cer1adc=csih->GetKovADC1(); // too
          Float_t cer1time= toftime - 
              (csih->GetKovTimeStamp1()*tdc_convert -abstime/Munits::ns);//too
          Int_t cer2adc=csih->GetKovADC3();
          Float_t cer2time= toftime - 
               (csih->GetKovTimeStamp3()*tdc_convert-abstime/Munits::ns);

          MSG("CalDetPID",Msg::kDebug)<<"toftime - (kov2*tdc_convert -abstime/Munits::ns)= "
                                         <<toftime<<"- ("
                                         <<csih->GetKovTimeStamp2()*tdc_convert
                                         <<" - "<<abstime/Munits::ns<<")"
                                         <<" = "<<toftime -(csih->GetKovTimeStamp2()*tdc_convert -abstime/Munits::ns)
                                         <<endl;
          
          MSG("CalDetPID",Msg::kDebug)<<"toftime - (kov1*tdc_convert -abstime/Munits::ns)= "
                                         <<toftime<<"- ("
                                         <<csih->GetKovTimeStamp1()*tdc_convert
                                         <<" - "<<abstime/Munits::ns<<")"
                                         <<" = "<<toftime -(csih->GetKovTimeStamp1()*tdc_convert -abstime/Munits::ns)
                                         <<endl;

          MSG("CalDetPID",Msg::kDebug)<<"toftime - (kov3*tdc_convert -abstime/Munits::ns)= "
                                         <<toftime<<"- ("
                                         <<csih->GetKovTimeStamp3()*tdc_convert
                                         <<" - "<<abstime/Munits::ns<<")"
                                         <<" = "<<toftime -(csih->GetKovTimeStamp3()*tdc_convert -abstime/Munits::ns)
                                         <<endl;
          const float nullwin=100000;
          // set cerX true if counter X passes the selection
          if(cer0adc>0){
              if(win0low==-nullwin && win0high==nullwin) cer0=true;
              else if((cer0time<win0low) || (cer0time>win0high)) cer0=false;
              else cer0=true;
          }
          else cer0=true;

          if(cer1adc>0){
              if(win1low==-nullwin && win1high==nullwin) cer1=true;
              else if((cer1time<win1low) || (cer1time>win1high)) cer1=false;
              else cer1=true;
          }
          else cer1=true;

          if(cer2adc>0){
              if(win2low==-nullwin && win2high==nullwin) cer2=true;
              else if((cer2time<win2low) || (cer2time>win2high)) cer2=false;
              else cer2=true;
          }
          else cer2=true;
          
          return (cer0&cer1&cer2);
     }    

}

Float_t AlgCalDetPID::OverlapLikelihood(TH1F *ehfar, TH1F *ehnear) const
{
   //now loop over the histogram bins to calculate the chi2
   Int_t nhits=0;
   Float_t flsum=0.;
   
   TH1F *eventhist[2];
   eventhist[0] = ehfar;
   eventhist[1] = ehnear;

   for(int i=0;i<2;i++){
      const int nevents = (const int)eventhist[i]->Integral();
      nhits+=nevents;
      for(int j=1;j<=eventhist[i]->GetNbinsX();j++){
         if(eventhist[i]->GetBinContent(j)>0){
            int x = (int)eventhist[i]->GetBinContent(j);
            float mean = fThist[i]->GetBinContent(j)*nevents;
            MSG("CalDetPID",Msg::kDebug)<<" crate "<<i
                                   <<" at time "<<eventhist[i]->GetBinCenter(j)
                                   <<" x "<<x<<" mean "<<mean<<endl;

            if(mean==0){
               continue;
            }

            double lnprob =  1.*x*log(mean)-mean;
            double lnnorm = 1.*x*log(1.*x)-1.*x;

            flsum+=(lnprob-lnnorm);
            
            MSG("CalDetPID",Msg::kDebug)<<"i "<<i
                                   <<" ln prob "<<lnprob
                                   <<" ln norm "<<lnnorm
                                   <<" fl sum "<<flsum<<endl;
         }
      }
      MSG("CalDetPID",Msg::kDebug)<<"i "<<i
                             <<" flsum "<<flsum<<endl;
   }
   if(nhits>0){
      flsum = -2*(flsum/(1.*nhits));  
   }
   else{
      flsum = 0;
   }
   return flsum;
}

---