MatacqProducer.cc

Go to the documentation of this file.

#include "EventFilter/EcalRawToDigi/interface/MatacqProducer.h"
#include "EventFilter/EcalRawToDigi/src/Majority.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "DataFormats/FEDRawData/interface/FEDRawData.h"
#include "DataFormats/FEDRawData/interface/FEDNumbering.h"
#include "DataFormats/EcalDigi/interface/EcalMatacqDigi.h"
#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
#include <boost/algorithm/string.hpp>
#include <boost/format.hpp>

#include <iostream>
#include <memory>

#include <stdio.h>

#include <fstream>
#include <iomanip>

#include "DataFormats/EcalDigi/interface/EcalMatacqDigi.h"

#include <sys/types.h>
#include <unistd.h>
#include <signal.h>
#include <sys/stat.h>

using namespace std;
using namespace boost;
using namespace edm;

// #undef LogInfo
// #define LogInfo(a) cout << "INFO " << a << ": "
// #undef LogWarning
// #define LogWarning(a) cout << "WARN " << a << ": "
// #undef LogDebug
// #define LogDebug(a) cout << "DBG " << a << ": "


//verbose mode for matacq event retrieval debugging:
//static const bool searchDbg = false;

//laser freq is 1 every 112 orbit => >80 orbit
const int MatacqProducer::orbitTolerance_ = 80; 

const MatacqProducer::stats_t MatacqProducer::stats_init = {0,0,0};

static std::string now(){
  struct timeval t;
  gettimeofday(&t, 0);
 
  char buf[256];
  strftime(buf, sizeof(buf), "%F %R %S s", localtime(&t.tv_sec));
  buf[sizeof(buf)-1] = 0;

  stringstream buf2;
  buf2 << buf << " " << ((t.tv_usec+500)/1000)  << " ms";

  return buf2.str();
}


MatacqProducer::MatacqProducer(const edm::ParameterSet& params):
  fileNames_(params.getParameter<std::vector<std::string> >("fileNames")),
  digiInstanceName_(params.getParameter<string>("digiInstanceName")),
  rawInstanceName_(params.getParameter<string>("rawInstanceName")),
  timing_(params.getUntrackedParameter<bool>("timing", false)),
  disabled_(params.getParameter<bool>("disabled")),
  verbosity_(params.getUntrackedParameter<int>("verbosity", 0)),
  produceDigis_(params.getParameter<bool>("produceDigis")),
  produceRaw_(params.getParameter<bool>("produceRaw")),
  inputRawCollection_(params.getParameter<edm::InputTag>("inputRawCollection")),
  mergeRaw_(params.getParameter<bool>("mergeRaw")),
  ignoreTriggerType_(params.getParameter<bool>("ignoreTriggerType")),
  matacq_(0, 0),
  inFile_(0),
  data_(bufferSize),
  openedFileRunNumber_(0),
  lastOrb_(0),
  fastRetrievalThresh_(0),
  orbitOffsetFile_(params.getUntrackedParameter<std::string>("orbitOffsetFile",
                                 "")),
  inFileName_(""),
  stats_(stats_init),
  logFileName_(params.getUntrackedParameter<std::string>("logFileName",
                             "matacqProducer.log")),
  eventSkipCounter_(0),
  onErrorDisablingEvtCnt_(params.getParameter<int>("onErrorDisablingEvtCnt")),
  timeLogFile_(params.getUntrackedParameter<std::string>("timeLogFile", "")),
  runNumber_(0)
{
  if(verbosity_>=4) cout << "[Matacq " << now() << "] in MatacqProducer ctor"  << endl;
  
  gettimeofday(&timer_, 0);

  if(timeLogFile_.size()>0){
    timeLog_.open(timeLogFile_.c_str());
    if(timeLog_.fail()){
      cout << "[LaserSorter " << now() << "] "
           << "Failed to open file " << timeLogFile_ << " to log timing.\n";
      logTiming_ = false;
    } else{
      logTiming_ = true;
    }
  }
  
  posEstim_.verbosity(verbosity_);

  logFile_.open(logFileName_.c_str(), ios::app | ios::out);
  
  if(logFile_.bad()){
    throw cms::Exception("FileOpen") << "Failed to open file "
                     << logFileName_ << " for logging.\n";
  }

  inputRawCollectionToken_ = consumes<FEDRawDataCollection>(params.getParameter<InputTag>("inputRawCollection"));


  if(produceDigis_){
    if(verbosity_>0) cout << "[Matacq " << now() << "] registering new "
                       "EcalMatacqDigiCollection product with instance name '"
                          << digiInstanceName_ << "'\n";
    produces<EcalMatacqDigiCollection>(digiInstanceName_);
  }
  
  if(produceRaw_){
    if(verbosity_>0) cout << "[Matacq " << now() << "] registering new FEDRawDataCollection "
                       "product with instance name '"
                          << rawInstanceName_ << "'\n";
    produces<FEDRawDataCollection>(rawInstanceName_);
  }
  
  startTime_.tv_sec = startTime_.tv_usec = 0;
  if(orbitOffsetFile_.size()>0){
    doOrbitOffset_ = true;
    loadOrbitOffset();
  } else{
    doOrbitOffset_ = false;
  }
  if(verbosity_>=4) cout << "[Matacq " << now() << "] exiting MatacqProducer ctor"  << endl;
}


void
MatacqProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup){
  if(verbosity_>=4) cout << "[Matacq " << now() << "] in MatacqProducer::produce"  << endl;
  if(logTiming_){
    timeval t;
    gettimeofday(&t, 0);

    timeLog_ << t.tv_sec << "."
             << setfill('0') << setw(3) << (t.tv_usec+500)/1000 << setfill(' ')<< "\t"
             << (t.tv_usec - timer_.tv_usec)*1. 
      + (t.tv_sec - timer_.tv_sec)*1.e6 << "\t";
    timer_ = t;
  } 
 
  if(startTime_.tv_sec==0) gettimeofday(&startTime_, 0);
  ++stats_.nEvents;  
  if(disabled_) return;
  const uint32_t runNumber = getRunNumber(event);
  if(runNumber!=runNumber_){
    newRun(runNumber_, runNumber);
  }
  addMatacqData(event);

  if(logTiming_){
    timeval t;
      gettimeofday(&t, 0);
      timeLog_ << (t.tv_usec - timer_.tv_usec)*1. 
        + (t.tv_sec - timer_.tv_sec)*1.e6 << "\n";
      timer_ = t;
  }
}

void
MatacqProducer::addMatacqData(edm::Event& event){

  edm::Handle<FEDRawDataCollection> sourceColl;
  event.getByToken(inputRawCollectionToken_, sourceColl);

  std::auto_ptr<FEDRawDataCollection> rawColl;
  if(produceRaw_){
    if(mergeRaw_){
      rawColl = auto_ptr<FEDRawDataCollection>(new FEDRawDataCollection(*sourceColl));
    } else{
      rawColl = auto_ptr<FEDRawDataCollection>(new FEDRawDataCollection());
    }
  }
  
  std::auto_ptr<EcalMatacqDigiCollection>
    digiColl(new EcalMatacqDigiCollection());
  
  if(eventSkipCounter_==0){
    if(sourceColl->FEDData(matacqFedId_).size()>4 && !produceRaw_){
      //input raw data collection already contains matacqData
      formatter_.interpretRawData(sourceColl->FEDData(matacqFedId_),
                  *digiColl);             
    } else{
      bool isLaserEvent = (getCalibTriggerType(event) == laserType);


      //      cout << "---> " << (ignoreTriggerType_?"yes":"no") << " " << getCalibTriggerType(event) << endl;
      
      if(isLaserEvent || ignoreTriggerType_){
    
    const uint32_t runNumber = getRunNumber(event);
    const uint32_t orbitId   = getOrbitId(event);
      
    LogInfo("Matacq") << "Run " << runNumber << "\t Orbit " << orbitId << "\n";
      
    bool fileChange;
    if(doOrbitOffset_){
      map<uint32_t,uint32_t>::iterator it = orbitOffset_.find(runNumber);
      if(it == orbitOffset_.end()){
        LogWarning("Matacq") << "Orbit offset not found for run "
                 << runNumber
                 << ". No orbit correction will be applied.";
      }
    }  
      
    if(getMatacqFile(runNumber, orbitId, &fileChange)){
      //matacq file retrieval succeeded
      LogInfo("Matacq") << "Matacq data file found for "
                << "run " << runNumber << " orbit " << orbitId;
      if(getMatacqEvent(runNumber, orbitId, fileChange)){
        if(produceDigis_){
          formatter_.interpretRawData(matacq_, *digiColl);
        }
        if(produceRaw_){
          uint32_t dataLen64 = matacq_.getParsedLen();
          if(dataLen64 > bufferSize*8 || matacq_.getDccLen()!= dataLen64){
        LogWarning("Matacq") << " Error in Matacq event fragment length! "
                     << "DCC len: " << matacq_.getDccLen()
                     << "*8 Bytes, Parsed len: "
                     << matacq_.getParsedLen() << "*8 Bytes.  "
                     << "Matacq data will not be included for this event.\n";
          } else{
        rawColl->FEDData(matacqFedId_).resize(dataLen64*8);
        copy(data_.begin(), data_.begin() + dataLen64*8,
             rawColl->FEDData(matacqFedId_).data());
          }
        }
        LogInfo("Matacq") << "Associating matacq data with orbit id "
                  << matacq_.getOrbitId()
                  << " to dcc event with orbit id "
                  << orbitId << std::endl;
        if(isLaserEvent){
          ++stats_.nLaserEventsWithMatacq;
        } else{
          ++stats_.nNonLaserEventsWithMatacq;
        }
      } else{
        if(isLaserEvent){
          LogWarning("Matacq") << "No matacq data found for laser event "
                   << "of run " << runNumber << " orbit "
                   << orbitId;
        }
      }
    } else{
      LogWarning("Matacq") << "No matacq file found for event "
                   << event.id();
    }
      }
    }
    if(eventSkipCounter_>0){ //error occured for this events
      //                       and some events will be skipped following
      //                       to this error.
      LogInfo("Matacq") << " [" << now() << "] "
            << eventSkipCounter_
            << " next events will be skipped, following to an "
            << "error on the last processed event, "
            << "which is expected to be persistant.";
    }
  } else{
    --eventSkipCounter_;
  }
  
  if(produceRaw_){
    if(verbosity_>1) cout << "[Matacq " << now() << "] "
                          << "Adding FEDRawDataCollection collection "
                          << " to event.\n";
    event.put(rawColl, rawInstanceName_);
  }

  if(produceDigis_){
    if(verbosity_>1) cout << "[Matacq " << now() << "] "
                          << "Adding EcalMatacqDigiCollection collection "
                          << " to event.\n";
    event.put(digiColl, digiInstanceName_);
  }
}

// #if 0
// bool
// MatacqProducer::getMatacqEvent(std::ifstream& f,
//                 uint32_t runNumber,
//                 uint32_t orbitId,
//                 bool doWrap,
//                 std::streamoff maxPos){
//   bool found = false;
//   streampos startPos = f.tellg();

//   while(!f.eof()
//  && !found
//  && (maxPos<0 || f.tellg()<=maxPos)){
//     const streamsize headerSize = 8*8;
//     f.read((char*)&data_[0], headerSize);
//     if(f.eof()) break;
//     int32_t orb = MatacqRawEvent::getOrbitId(&data_[0], headerSize);
//     uint32_t len = MatacqRawEvent::getDccLen(&data_[0], headerSize);
//     uint32_t run = MatacqRawEvent::getRunNum(&data_[0], headerSize);
//     //     cout << "Matacq: orbit = " << orb
//     //    << " len = " << len
//     //    << " run = " << run << endl;
//     if((abs(orb-(int32_t)orbitId) < orbitTolerance_)
//        && (runNumber==0 || runNumber==run)){
//       found = true;
//       //reads the rest of the event:
//       if(data_.size() < len*8){
//  throw cms::Exception("Matacq") << "Buffer overflow";
//       }
//       f.read((char*)&data_[0]+headerSize, len*8-headerSize);
//       matacq_ = MatacqRawEvent((unsigned char*)&data_[0], len*8);
//     } else{
//       //moves to next event:
//       f.seekg(len*8 - headerSize, ios::cur);
//     }
//   }
  
//   f.clear(); //clears eof error to allow seekg  
//   if(doWrap && !found){
//     f.seekg(0, ios::beg);
//     found =  getMatacqEvent(f, runNumber, orbitId, false, startPos);
//   }
//   return found;
// }
//#endif

bool
MatacqProducer::getMatacqEvent(uint32_t runNumber,
                   int32_t orbitId,
                   bool fileChange){
  filepos_t startPos;
  if(!mtell(startPos)) return false;
  
  int32_t startOrb = -1;
  const size_t headerSize = 8*8;
  if(mread((char*)&data_[0], headerSize, "Reading matacq header", true)){
    startOrb = MatacqRawEvent::getOrbitId(&data_[0], headerSize);
    if(startOrb<0) startOrb = 0;
  } else{
    if(verbosity_>2){
      cout << "[Matacq " << now() << "] Failed to read matacq header. Moved to start of "
    " the file.\n";
    }
    mrewind();
    if(mread((char*)&data_[0], headerSize, "Reading matacq header", true)){
      startPos = 0;
      startOrb = MatacqRawEvent::getOrbitId(&data_[0], headerSize);
    } else{
      if(verbosity_>2) cout << "[Matacq " << now() << "] Looks like matacq file is empty"
                            << "\n";
      return false;
    }
  }
  
  if(verbosity_>2) cout << "[Matacq " << now() << "] Last read orbit: " << lastOrb_
                        << " looking for orbit " << orbitId
                        << ". Current file position: " << startPos
                        << " Orbit at current position: " << startOrb << "\n";
  
  //  f.clear();
  bool didCoarseMove = false;

  //FIXME: case where posEtim_.invalid() is false
  if(!posEstim_.invalid()
     && (abs(lastOrb_-orbitId) > fastRetrievalThresh_)){
    filepos_t pos = posEstim_.pos(orbitId);

    //    struct stat st;
    filepos_t fsize;
    //    if(0==stat(inFileName_.c_str(), &st)){
    if(msize(fsize)){
      //      const int64_t fsize = st.st_size;
      if(0!=posEstim_.eventLength() && pos > fsize){
    //estimated position is beyong end of file
    //-> move to beginning of last event:
    int64_t evtSize = posEstim_.eventLength()*sizeof(uint64_t);
    pos = ((int64_t)fsize/evtSize-1)*evtSize;
    if(verbosity_>2){
      cout << "[Matacq " << now() << "] Estimated position was beyond end of file. "
        "Changed to " << pos << "\n";
    }
      }
    } else{
      LogWarning("Matacq") << "Failed to access file " << inFileName_ << ".";
    }
    if(pos>=0){
      if(verbosity_>2) cout << "[Matacq " << now() << "] jumping to estimated position "
                << pos << "\n";
      mseek(pos, SEEK_SET, "Jumping to estimated event position");
      if(mread((char*)&data_[0], headerSize, "Reading matacq header", true)){
    didCoarseMove = true;
      } else{
    //estimated position might have been beyond the end of the file,
    //try, with original position:
    didCoarseMove = false;
    if(!mread((char*)&data_[0], headerSize, "Reading event header", true)){
      return false;
    }
      }
    } else{
      if(verbosity_) cout << "[Matacq " << now() << "] Event orbit outside of orbit range "
               "of matacq data file events\n";
      return false;
    }
  }
  
  int32_t orb = MatacqRawEvent::getOrbitId(&data_[0], headerSize);

  if(didCoarseMove){
    //autoadjustement of threshold for coarse move:
    if(abs(orb-orbitId) > fastRetrievalThresh_){
      if(verbosity_>2) cout << "[Matacq " << now() << "] Fast retrieval threshold increased from "
                            << fastRetrievalThresh_;
      fastRetrievalThresh_ = 2*abs(orb-orbitId);
      if(verbosity_>2) cout << " to " << fastRetrievalThresh_ << "\n";
    }

    //if coarse move did not improve situation, rolls back:
    if(startOrb > 0
       && (abs(orb-orbitId) > abs(startOrb-orbitId))){
      if(verbosity_>2) cout << "[Matacq " << now() << "] Estimation (-> orbit " << orb << ") "
                         "was worst than original position (-> orbit "
                            << startOrb
                            << "). Restoring position (" << startPos << ").\n";
      mseek(startPos, SEEK_SET);
      mread((char*)&data_[0], headerSize, "Reading event header", true);
      orb = MatacqRawEvent::getOrbitId(&data_[0], headerSize);
    }
  }
  
  bool searchBackward = (orb>orbitId)?true:false;
  //BEWARE: len must be signed, because we are using latter in the code (-len)
  //expression
  int len = (int)MatacqRawEvent::getDccLen(&data_[0], headerSize);

  if(len==0){
    cout << "[Matacq " << now() << "] read DCC length is null! Cancels matacq event search "
     << " and move matacq file pointer to beginning of the file. "
     << "(" << __FILE__ << ":" << __LINE__ << ")."
     << "\n";
    //rewind(f);
    mrewind();
    return false;
  }
  
  enum state_t { searching, found, failed } state = searching;
  
  while(state == searching){
    orb = MatacqRawEvent::getOrbitId(&data_[0], headerSize);
    len = (int)MatacqRawEvent::getDccLen(&data_[0], headerSize);
    uint32_t run = MatacqRawEvent::getRunNum(&data_[0], headerSize);
    if(verbosity_>3){
      filepos_t pos = -1;
      mtell(pos);
      cout << "[Matacq " << now() << "] Header read at file position "
       << pos
       << ":  orbit = " << orb
       << " len = " << len << "x8 Byte"
       << " run = " << run << "\n";
    }
    if((abs(orb-orbitId) < orbitTolerance_)
       && (runNumber==0 || runNumber==run)){
      state = found;
      lastOrb_ = orb;
      //reads the rest of the event:
      if((int)data_.size() < len*8){
    throw cms::Exception("Matacq") << "Buffer overflow";
      }
      if(verbosity_>2) cout << "[Matacq " << now() << "] Event found. Reading "
                         " matacq event." << "\n";
      if(!mread((char*)&data_[0], len*8, "Reading matacq event")){
    if(verbosity_>2) cout << "[Matacq " << now() << "] Failed to read matacq event."
                              << "\n";
    state = failed;
      }
      matacq_ = MatacqRawEvent((unsigned char*)&data_[0], len*8);
    } else {
      if((searchBackward && (orb < orbitId))
     || (!searchBackward && (orb > orbitId))){ //search ended
    lastOrb_ = orb;      
    state = failed;
    if(verbosity_>2) cout << "[Matacq " << now()
                << "] No matacq data found for run " << run
                << ", orbit ID " << orbitId << "." << "\n";
      } else{
    off_t offset = (searchBackward?-len:len)*8;
    lastOrb_ = orb; 
    if(verbosity_>3){
      cout << "[Matacq " << now() << "] In matacq file, moving "
           << abs(offset) << " byte " << (offset>0?"forward":"backward")
           << ".\n";
    }   

    if(mseek(offset, SEEK_CUR,
         (searchBackward?"Moving to previous event":
          "Moving to next event"))
       && mread((char*)&data_[0], headerSize, "Reading event header",
            true)){
    } else{
      if(!searchBackward) mseek(-len*8, SEEK_CUR,
                    "Moving to start of last complete event");
      state = failed;
    }
      }
    }
  }
  
  if(state==found){
    filepos_t pos = -1;
    filepos_t fsize = -1;
    mtell(pos);
    msize(fsize);
    if(pos==fsize-1){ //last byte.
      if(verbosity_>2){
    cout << "[Matacq " << now() << "] Event found was at the end of the file. Moving "
      "stream position to beginning of this event."
         << "\n";
      }
      mseek(-(int)len*8-1, SEEK_CUR,
        "Moving to beginning of last matacq event");
    }
  }
  return (state==found);
}


bool
MatacqProducer::getMatacqFile(uint32_t runNumber, uint32_t orbitId,
                  bool* fileChange){
  if(openedFileRunNumber_!=0
     && openedFileRunNumber_==runNumber){
    if(fileChange!=0) *fileChange = false;
    return misOpened();
  }

  if(fileNames_.size()==0) return 0;

  const string runNumberFormat = "%08d";
  string sRunNumber = str(boost::format(runNumberFormat) % runNumber);
  //cout << "Run number string: " << sRunNumber << "\n";
  bool found = false;
  string fname;
  for(unsigned i=0; i < fileNames_.size() && !found; ++i){
    fname = fileNames_[i];
    boost::algorithm::replace_all(fname, "%run_subdir%",
                  runSubDir(runNumber));
    boost::algorithm::replace_all(fname, "%run_number%", sRunNumber);

    if(verbosity_>0) cout << "[Matacq " << now() << "] "
              << "Looking for a file with path "
              << fname << "\n";
    
    if(mcheck(fname)){
      LogInfo("Matacq") << "Uses matacq data file: '" << fname << "'\n";
      found = true;
    }
  }
  if(!found){
    if(verbosity_>=0) cout << "[Matacq " << now() << "] no matacq file found "
            "for run " << runNumber << "\n";
    eventSkipCounter_ = onErrorDisablingEvtCnt_;
    openedFileRunNumber_ = 0;
    if(fileChange!=0) *fileChange = false;
    return 0;
  }
  
  if(!mopen(fname)){
    LogWarning("Matacq") << "Failed to open file " << fname << "\n";
    eventSkipCounter_ = onErrorDisablingEvtCnt_;
    openedFileRunNumber_ = 0;
    if(fileChange!=0) *fileChange = false;
    return false;
  } else{
    openedFileRunNumber_ = runNumber;
    lastOrb_ = 0;
    posEstim_.init(this);
    if(fileChange!=0) *fileChange = true;
    return true;
  }
}
 

uint32_t MatacqProducer::getRunNumber(edm::Event& ev) const{
  return ev.run();
}

uint32_t MatacqProducer::getOrbitId(edm::Event& ev) const{
  //on CVS HEAD (June 4, 08), class Event has a method orbitNumber()
  //we could use here. The code would be shorten to:
  //return ev.orbitNumber();
  //we have to deal with what we have in current CMSSW releases:
  edm::Handle<FEDRawDataCollection> rawdata;
  ev.getByToken(inputRawCollectionToken_, rawdata);
  if(!(rawdata.isValid())){
    throw cms::Exception("NotFound")
      << "No FED raw data collection found. ECAL raw data are "
      "required to retrieve the orbit ID";
  }
  
  int orbit = 0;
  for(int id=601; id<=654; ++id){
    if(!FEDNumbering::inRange(id)) continue;
    const FEDRawData& data = rawdata->FEDData(id);
    const int orbitIdOffset64 = 3;
    if(data.size()>=8*(orbitIdOffset64+1)){//orbit id is in 4th 64-bit word
      const unsigned char* pOrbit = data.data() + orbitIdOffset64*8;
      int thisOrbit = pOrbit[0]
    | (pOrbit[1] <<8)
    | (pOrbit[2] <<16)
    | (pOrbit[3] <<24);
      if(orbit!=0 && thisOrbit!=0 && abs(orbit-thisOrbit)>orbitTolerance_){
    //throw cms::Exception("EventCorruption")
    //  << "Orbit ID inconsitency in DCC headers";
    LogWarning("EventCorruption")
      << "Orbit ID inconsitency in DCC headers";
    orbit = 0;
    break;
      }
      if(thisOrbit!=0) orbit = thisOrbit;
    }
  }
  
  if(orbit==0){
    //    throw cms::Exception("NotFound")
    //  << "Failed to retrieve orbit ID of event "<< ev.id();
    LogWarning("NotFound") << "Failed to retrieve orbit ID of event "
                << ev.id();
  }
  return orbit;
}

int MatacqProducer::getCalibTriggerType(edm::Event& ev) const{  
  edm::Handle<FEDRawDataCollection> rawdata;
  ev.getByToken(inputRawCollectionToken_, rawdata);
  if(!(rawdata.isValid())){
    throw cms::Exception("NotFound")
      << "No FED raw data collection found. ECAL raw data are "
      "required to retrieve the trigger type";
  }
  
  Majority<int> stat;
  for(int id=601; id<=654; ++id){
    if(!FEDNumbering::inRange(id)) continue;
    const FEDRawData& data = rawdata->FEDData(id);
    const int detailedTrigger32 = 5;
    if(data.size()>=4*(detailedTrigger32+1)){
      const unsigned char* pTType = data.data() + detailedTrigger32*4;
      int tType = pTType[1] & 0x7;
      stat.add(tType);
    }
  }
  double p;
  int tType = stat.result(&p);
  if(p<0){
    //throw cms::Exception("NotFound") << "No ECAL DCC data found\n";
    LogWarning("NotFound")  << "No ECAL DCC data found\n";
    tType = -1;
  }
  if(p<.8){
    //throw cms::Exception("EventCorruption") << "Inconsitency in detailed trigger type indicated in ECAL DCC data headers\n";
    LogWarning("EventCorruption") << "Inconsitency in detailed trigger type indicated in ECAL DCC data headers\n";
    tType = -1;
  }
  return tType;
}
  
void MatacqProducer::PosEstimator::init(MatacqProducer* mp){
  mp->mrewind();

  const size_t headerSize = 8*8;
  unsigned char data[headerSize];
  if(!mp->mread((char*)data, headerSize)){
    if(verbosity_) cout << "[Matacq " << now() << "] reached end of file!\n"; 
    firstOrbit_ = eventLength_ = orbitStepMean_ = 0;
    return;
  } else{
    firstOrbit_ = MatacqRawEvent::getOrbitId(data, headerSize);
    eventLength_ = MatacqRawEvent::getDccLen(data, headerSize);
    if(verbosity_>1) cout << "[Matacq " << now() << "] First event orbit: " << firstOrbit_
                          << " event length: " << eventLength_
                          << "*8 byte\n";
  }

  mp->mrewind();
  
  if(eventLength_==0){    
    if(verbosity_) cout << "[Matacq " << now() << "] event length is null!" << endl; 
    return;
  }

  filepos_t s;
  mp->msize(s);

  //number of complete events:
  const unsigned nEvents = s/eventLength_/8;

  if(nEvents==0){
    if(verbosity_) cout << "[Matacq " << now() << "] File is empty!" << endl;
    orbitStepMean_ = 0;
    return;
  }

  if(verbosity_>1) cout << "[Matacq " << now() << "] File size: " << s
                        << " Number of events: " << nEvents << endl;
  
  //position of last complete events:
  off_t last = (nEvents-1)*(off_t)eventLength_*8;
  mp->mseek(last, SEEK_SET, "Moving to beginning of last complete "
        "matacq event");
  if(!mp->mread((char*) data, headerSize, "Reading matacq header", true)){
    LogWarning("Matacq") << "Fast matacq event retrieval failure. "
      "Falling back to safe retrieval mode.";
    orbitStepMean_ = 0;
  }
  
  int32_t lastOrb = MatacqRawEvent::getOrbitId(data, headerSize);
  int32_t lastLen = MatacqRawEvent::getDccLen(data, headerSize);

  if(verbosity_>1) cout << "[Matacq " << now() << "] Last event orbit: " << lastOrb
                        << " last event length: " << lastLen << endl;
  
  //some consistency check
  if(lastLen!=eventLength_){
    LogWarning("Matacq")
      //throw cms::Exception("Matacq")
      << "Fast matacq event retrieval failure: it looks like "
      "the matacq file contains events of different sizes.";
      //      " Falling back to safe retrieval mode.";
    invalid_ = false; //true;
    orbitStepMean_ = 112; //0;
    return;
  }

  orbitStepMean_ = (lastOrb - firstOrbit_)/nEvents;
  
  if(verbosity_>1) cout << "[Matacq " << now() << "] Orbit step mean: " << orbitStepMean_
                        << "\n";

  invalid_ = false;
}

int64_t MatacqProducer::PosEstimator::pos(int orb) const{
  if(orb<firstOrbit_) return -1;
  uint64_t r = orbitStepMean_!=0?
    (((uint64_t)(orb-firstOrbit_))/orbitStepMean_)*eventLength_*8
    :0;
  if(verbosity_>2) cout << "[Matacq " << now() << "] Estimated Position for orbit  " << orb
                        << ": " << r << endl;
  return r;
}

MatacqProducer::~MatacqProducer(){
  mclose();
  timeval t;
  gettimeofday(&t, 0);
  if(logTiming_ && startTime_.tv_sec!=0){
    //not using logger, to allow timing with different logging options
    cout << "[Matacq " << now() << "] Time elapsed between first event and "
      "destruction of MatacqProducer: "
     << ((t.tv_sec-startTime_.tv_sec)*1.
         + (t.tv_usec-startTime_.tv_usec)*1.e-6) << "s\n";
  }
}

void MatacqProducer::loadOrbitOffset(){
  std::ifstream f(orbitOffsetFile_.c_str());
  if(f.bad()){
    throw cms::Exception("Matacq")
      << "Failed to open orbit ID correction file '"
      << orbitOffsetFile_ << "'\n";
  }

  cout << "[Matacq " << now() << "] "
       << "Offset to substract to Matacq events Orbit ID: \n"
       << "#Run Number\t Offset\n";

  int iline = 0;
  string s;
  stringstream buf;
  while(f.eof()){
    getline(f, s);
    ++iline;
    if(s[0]=='#'){//comment
      //skip line:
      f.ignore(numeric_limits<streamsize>::max(), '\n');
      continue;
    }
    buf.str("");
    buf << s;
    int run;
    int orbit;
    buf >> run;
    buf >> orbit;
    if(buf.bad()){
      throw cms::Exception("Matacq")
    << "Syntax error in Orbit offset file '"
    << orbitOffsetFile_ << "'";
    }
    cout << run << "\t" << orbit << "\n";
    orbitOffset_.insert(pair<int, int>(run, orbit));
  }
}

#ifdef USE_STORAGE_MANAGER
bool MatacqProducer::mseek(filepos_t offset, int whence, const char* mess){
  if(0==inFile_.get()) return false;
  try{
    Storage::Relative wh;
    if(whence==SEEK_SET) wh = Storage::SET;
    else if(whence==SEEK_CUR) wh = Storage::CURRENT;
    else if(whence==SEEK_END) wh = Storage::END;
    else throw cms::Exception("Bug") << "Bug found in "
                                     << __FILE__ << ": "<< __LINE__ << "\n";
                   
    inFile_->position(offset, wh);
  } catch(cms::Exception& e){
    if(verbosity_){
      cout << "[Matacq " << now() << "] ";
      if(mess) cout << mess << ". ";
      cout << "Random access error on input matacq file. ";
      if(whence==SEEK_SET) cout << "Failed to seek absolute position " << offset;
      else if(whence==SEEK_CUR) cout << "Failed to move " << offset << " bytes forward";
      else if(whence==SEEK_END) cout << "Failed to seek position at " << offset << " bytes before end of file";
    cout << ". Reopening file. " << e.what() << "\n";
      mopen(inFileName_);
      return false;
    }
  }
  return true;
}

bool MatacqProducer::mtell(filepos_t& pos){
  if(0==inFile_.get()) return false;
  pos = inFile_->position();
  return true;
}

bool MatacqProducer::mread(char* buf, size_t n, const char* mess, bool peek){
  if(0==inFile_.get()) return false;
  
  filepos_t pos = -1;
  if(!mtell(pos)) return false;

  bool rc = false;
  try{
    rc =  (n==inFile_->xread(buf, n));
  } catch(cms::Exception& e){
    if(verbosity_){
      cout << "[Matacq " << now() << "] ";
      if(mess) cout << mess << ". ";
      cout << "Read failure from input matacq file: "
           << e.what() << "\n";
    }
    //recovering from error:
    mopen(inFileName_);
    mseek(pos);
    return false;
  }
  if(peek){//asked to restore original file position
    mseek(pos);
  }
  return rc;
}

bool MatacqProducer::msize(filepos_t& s){
  if(inFile_.get()==0) return false;
  s = inFile_.get()->size();
  return true;
}

bool MatacqProducer::mrewind(){
  Storage* file = inFile_.get();
  if(file==0) return false;
  try{
    file->rewind();
  } catch(cms::Exception e){
    if(verbosity_) cout << "Exception cautgh while rewinding file "
                        << inFileName_ << ": " << e.what() << ". "
                        << "File will be reopened.";
    return mopen(inFileName_);
  }
  return true;
}

bool MatacqProducer::mcheck(const std::string& name){
  return StorageFactory::get()->check(name);
}

bool MatacqProducer::mopen(const std::string& name){
  //close already opened file if any:
  mclose();
  
  try{
    inFile_
      = auto_ptr<Storage>(StorageFactory::get()->open(name,
                                                      IOFlags::OpenRead));
    inFileName_ = name;
  } catch(cms::Exception& e){
    LogWarning("Matacq") << e.what();
    inFile_.reset();
    inFileName_ = "";
    return false;
  }
  return true;
}

void MatacqProducer::mclose(){
  if(inFile_.get()!=0){
    inFile_->close();
    inFile_.reset();
  }
}

bool MatacqProducer::misOpened(){
  return inFile_.get()!=0;
}

bool MatacqProducer::meof(){
  if(inFile_.get()==0) return true;
  return inFile_->eof();
}

#else //USE_STORAGE_MANAGER not defined
bool MatacqProducer::mseek(off_t offset, int whence, const char* mess){
  if(0==inFile_) return false;
  const int rc = fseeko(inFile_, offset, whence);
  if(rc!=0 && verbosity_){
    cout << "[Matacq " << now() << "] ";
    if(mess) cout << mess << ". ";
    cout << "Random access error on input matacq file. "
      "Rewind file.\n";
    mrewind();
  }
  return rc==0;
}

bool MatacqProducer::mtell(filepos_t& pos){
  if(0==inFile_) return false;
  pos = ftello(inFile_);
  return pos != -1;
    
}

bool MatacqProducer::mread(char* buf, size_t n, const char* mess, bool peek){
  if(0==inFile_) return false;
  off_t pos = ftello(inFile_);
  bool rc = (pos!=-1) && (1==fread(buf, n, 1, inFile_));
  if(!rc){
    if(verbosity_){
      cout << "[Matacq " << now() << "] ";
      if(mess) cout << mess << ". ";
      cout << "Read failure from input matacq file.\n";
    }
    clearerr(inFile_);
  }
  if(peek || !rc){//need to restore file position
    if(0!=fseeko(inFile_, pos, SEEK_SET)){
      if(verbosity_){
    cout << "[Matacq " << now() << "] ";
    if(mess) cout << mess << ". ";
    cout << "Failed to restore file position of "
      "before read error. Rewind file.\n";
      }
      //rewind(inFile_.get());
      mrewind();
      lastOrb_ = 0;
    }
  }
  return rc;
}

bool MatacqProducer::msize(filepos_t& s){
  if(0==inFile_) return false;
  struct stat buf;
  if(0!=fstat(fileno(inFile_), &buf)){
    s = 0;
    return false;
  } else{
    s = buf.st_size;
    return true;
  }
}

bool MatacqProducer::mrewind(){
  if(0==inFile_) return false;
  clearerr(inFile_);
  return fseeko(inFile_, 0, SEEK_SET)!=0; 
}

bool MatacqProducer::mcheck(const std::string& name){
  struct stat dummy;
  return 0==stat(name.c_str(), &dummy);
//   if(stat(name.c_str(), &dummy)==0){
//     return true;
//   } else{
//     cout << "[Matacq " << now() << "] Failed to stat file '"
//   << name.c_str() << "'. " 
//   << "Error " << errno << ": " << strerror(errno) << "\n";
//     return false;
//   }
}

bool MatacqProducer::mopen(const std::string& name){
  if(inFile_!=0) mclose();
  inFile_ = fopen(name.c_str(), "r");
  if(inFile_!=0){
    inFileName_ = name;
    return true;
  } else{
    inFileName_ = "";
    return false;
  }
}

void MatacqProducer::mclose(){
  if(inFile_!=0) fclose(inFile_);
  inFile_ = 0;
}

bool MatacqProducer::misOpened(){
  return inFile_!=0;
}

bool MatacqProducer::meof(){
  if(0==inFile_) return true;
  return feof(inFile_)==0;
}

#endif //USE_STORAGE_MANAGER defined

std::string MatacqProducer::runSubDir(uint32_t runNumber){
  int millions = runNumber / (1000*1000);
  int thousands = (runNumber-millions*1000*1000) / 1000;
  int units = runNumber-millions*1000*1000 - thousands*1000;
  return str(boost::format("%03d/%03d/%03d") % millions % thousands % units);
}

void MatacqProducer::newRun(int prevRun, int newRun){
    runNumber_ = newRun;
    eventSkipCounter_ = 0;
    logFile_ << "[" << now() << "] Event count for run "
         << runNumber_ << ": "
         << "total: " << stats_.nEvents << ", "
         << "Laser event with Matacq data: "
         << stats_.nLaserEventsWithMatacq << ", "
         << "Non laser event (according to DCC header) with Matacq data: "
         << stats_.nNonLaserEventsWithMatacq << "\n" << flush;

    stats_.nEvents = 0;
    stats_.nLaserEventsWithMatacq = 0;
    stats_.nNonLaserEventsWithMatacq = 0;

    
}