HarpoReaderRaw.cxx

Go to the documentation of this file.

// 
// File HarpoReaderRaw.cxx
//
#include "HarpoReaderRaw.h"
#include "HarpoDetReader.h"
#include "HarpoConfig.h"
#include "HarpoRunHeader.h"
#include "HarpoEvent.h"
#include "HarpoTimeStamp.h"

#include "HarpoDccReader.h"
#include "HarpoFeminosReader.h"
#include "Pmm2Reader.h"
#include "HarpoSimReader.h"

#include "TSystem.h"

#include <climits>
#include <cstdio>
#include <arpa/inet.h>
#include <cstdlib>
//#include <cstring>
#include <iostream>


ClassImp(HarpoReaderRaw)


HarpoReaderRaw::HarpoReaderRaw() : HarpoReaderBase()
{
  if(gHarpoDebug>0) Info("HarpoReaderRaw()"," ");

  // fNewFile = NULL;
  // fNewFileName = (TString *) NULL;
  // fSave = kFALSE;

  fPmm2Empty = false;
  fSkipPmm2 = false;
  fSync = false;
  fPmm2Sync = 3;
  fMode = hSyncTime;
  fSkipEvents = 0;
  fRunHeader = NULL;
}


HarpoReaderRaw::HarpoReaderRaw(HarpoConfig *cfg) : HarpoReaderBase(cfg)
 {
  if(gHarpoDebug>0) Info("HarpoReaderRaw(cfg)"," ");

  // fNewFile = NULL;
  // fNewFileName = cfg->GetOutFile();
  // fSave = kFALSE;

  fPmm2Empty = false;
  fSkipPmm2 = false;
  fPmm2Sync = 3;

  fSkipEvents = gHConfig->GetSkipEvents();

  fSync = false;
  fMode = hSyncTime;
  fRunHeader = new HarpoRunHeader(cfg);

}

HarpoReaderRaw::~HarpoReaderRaw()
{
  if(gHarpoDebug>0) Info("~HarpoReaderRaw()"," ");
  if( fRunHeader != NULL) delete fRunHeader;
  // ??? Delete Det Constructors
  // if ( fNewFile != NULL) {
  //   fNewFile->Close();
  // }
}

void HarpoReaderRaw::SetDetFile(TString * /* name */ ,Int_t /* plane */)
{
  //todo Do we need this method ?
  //  fInFile[plane] = name; //???
}

 void HarpoReaderRaw::SetRootFile(TString *name)
 {
   //if(gHarpoDebug>0) 
   Info("Dummy SetRootFile","%s",name->Data());
 }
  

Bool_t HarpoReaderRaw::NextEvent()
{
  if(gHarpoDebug>1) Info("NextEvent"," ");
  Bool_t res = true;
  ULong_t ndet;
  Long_t deltat;
  //  ULong_t ts[2],tso[2];
 
  for(ndet=0;ndet<gkNDetectors;ndet++){
    if(gHarpoDebug>1) Info("NextEvent","%ld => %d",ndet,fDetSet->isExist(ndet));
    // if(fDetSet->isSim()){
    //   res = res && fDetReader[SIMDET]->NextEvent();
    //   break;
    // }
    if (fDetSet->isExist(ndet)){
      if(!(ndet==PMM2 && fPmm2Empty)) // ran out of PMM2 events (possible sync issue)
        res = res && fDetReader[ndet]->NextEvent();
      else{
        if(gHarpoDebug>0)
          Info("NextEvent","Skipping PMM2 Reader (no available data)");
      }
      //      fDetReader[ndet]->GetTimeStamp();
    }
  }

  if ( ! res ) return res; // No more event or errors

  //  if (fDetSet->isExist(0) && fDetSet->isExist(1) && !fDetSet->isSim()){
  if (fDetSet->isExist(0) && fDetSet->isExist(1)){
    ULong_t ts[2] = {0, 0};
    ULong_t tso[2] = {0, 0};
    for(ndet=0; ndet<2; ndet++) // Only 0 - plane X and 1 - plane Y
      {
        ts[ndet] = fDetReader[ndet]->GetTimeStamp();
        tso[ndet] = fDetReader[ndet]->GetPrevTimeStamp();
        if(gHarpoDebug>1) 
          Info("NextEvent","Time stamps for det %li: %li %li",ndet,ts[ndet],tso[ndet]);
      }
    
    deltat = HarpoTimeStamp::CompareTimeStamps(ts,tso);
    if(fNevts==0 && deltat<100000){
      HarpoTimeStamp::TSoffset = deltat;
      deltat = 0;
    }

    fSync = ( deltat == 0 );
    if(gHarpoDebug>1) 
      Info("NextEvent","deltat = %ld, fMode = %d",deltat,fMode);
    
    switch (fMode) {
    case hNoCMPTime :  // Read without timestamp comparision
      break;
    case hFirstTime : // return only plane with smaller time stamp
      if ( deltat == 0 ) {
        //fSync = true;
      } else if ( deltat < 0 ) { // x - plane with small time stamp
        fDetReader[1]->PushBack(); // the same y event on next read
      } else { // y smaler
        fDetReader[0]->PushBack();
      }
      break;
    case hSyncTime :   // skip plane event with smaller timestamp
      if ( deltat == 0 ) {
        // fSync = true;
      } else if ( deltat < 0 ) { // x - plane with small time stamp
        fDetReader[1]->PushBack(); // the same y event on next read
        fDetReader[0]->SkipEvent(); // next x
        res = NextEvent();
        fNevts--; // recursion, not real next event
      } else { // y smaler
        fDetReader[0]->PushBack();
        fDetReader[1]->SkipEvent(); // next y
        res = NextEvent();
        fNevts--; // recursion, not real next event
      }
      break;
    default :
      Warning("NextEvent","*** BAD Read Mode %d", fMode);
      return false;
    }
  }
  fNevts++;

  if ((fDetSet->isExist(0) || fDetSet->isExist(1)) && fDetSet->isExist(2) && !fPmm2Empty){

 
    // ULong_t ts[2] = {0, 0};
    // ULong_t tso[2] = {0, 0};
    // if(fDetSet->isExist(0)){
    //   ts[0] = (fDetReader[0]->GetTimeStamp() - fLatency)/10;
    //   tso[0] = (fDetReader[0]->GetPrevTimeStamp() - fLatency)/10;
    // }else{
    //   ts[0] = (fDetReader[1]->GetTimeStamp() - fLatency)/10;
    //   tso[0] = (fDetReader[1]->GetPrevTimeStamp() - fLatency)/10;
    // }
    // ts[1] = fDetReader[2]->GetTimeStamp();
    // tso[1] = fDetReader[2]->GetPrevTimeStamp();

    // deltat = HarpoTimeStamp::CompareTimeStamps(ts,tso);

    // Info("NextEvent","deltat: %li",deltat);


    ULong_t timestamp ;
    if(fDetSet->isExist(0))
      timestamp = fDetReader[0]->GetTimeStamp();
    else
      timestamp = fDetReader[1]->GetTimeStamp();

    if(timestamp)
      timestamp -= fLatency;
    ULong_t timestampPmm2 = 
      fDetReader[2]->GetEvent()->GetTimeStamp();
    if(Double_t(timestamp)/10 - Double_t(timestampPmm2) >= fPmm2Sync){
      //    if(TMath::Abs(deltat) > 2){
      fSkipPmm2 = true;
     if(gHarpoDebug>1)
       Info("NextEvent","Timestamp: %li > PMM2: %li (%g)",timestamp, timestampPmm2, Double_t(timestamp)/10 - Double_t(timestampPmm2));
      if(fDetSet->isExist(0))
        fDetReader[0]->PushBack();
      if(fDetSet->isExist(1))
        fDetReader[1]->PushBack();
      fDetReader[2]->SkipEvent();
      res = NextEvent();
      fNevts--;
      if(!res){
        Warning("NextEvent","Could not find matching PMM2 event for timestamp %li",timestamp);
        //      gHDetSet->UnSet(PMM2);
        fPmm2Empty = true;
        return true;
      }
    }else if(Double_t(timestamp)/10 - Double_t(timestampPmm2) <= -fPmm2Sync){
      //    if(TMath::Abs(deltat) > 2){
      if(gHarpoDebug>1)
        Info("NextEvent","Timestamp: %li < PMM2: %li  (%g)",timestamp, timestampPmm2, Double_t(timestamp)/10 - Double_t(timestampPmm2));

      fSkipPmm2 = true;

      // if(!res){
      //        Warning("NextEvent","Could not find matching PMM2 event for timestamp %li",timestamp);
      //        //      gHDetSet->UnSet(PMM2);
      //        fPmm2Empty = true;
      //        return true;
      // }
      // fNevts--;
    } else{
      if(gHarpoDebug>1)
        Info("NextEvent","******* Timestamp: %li , PMM2: %li ********",timestamp, timestampPmm2);
      fSkipPmm2 = false;
    }
  }

  FillEvent();

  //fSkipPmm2 = false;
  return res;
}

HarpoEvent *HarpoReaderRaw::FillEvent()
{
  if(gHarpoDebug>1) Info("GetEvent","event %ld",fNevts);
  HarpoDetEvent *revt;
  ULong_t ndet;

  /* if(fEvent==NULL) */ fEvent = new HarpoEvent();
  /* if(fEventHeader==NULL) */ fEventHeader = new HarpoEventHeader();

  fEventHeader->SetEvtNo(fNevts);  // Event header
  fEvent->SetHeader(fEventHeader);
  fEvent->SetRunHeader(fRunHeader);

  for(ndet=0;ndet<gkNDetectors;ndet++) {
    if (fDetSet->isExist(ndet)) {
      if(ndet==PMM2 && (fPmm2Empty || fSkipPmm2)){
        if(gHarpoDebug>1) Info("FillEvent","No PMM2 Data");
        revt = new Pmm2Event();
      }
      else
        revt = fDetReader[ndet]->GetEvent();
      fEvent->SetDetEvent(revt,ndet);
    }
  }
  
  return fEvent;
}

HarpoRunHeader *HarpoReaderRaw::GetRunHeader()
{
  if(gHarpoDebug>0) Info("GetRunHeader"," ");

  return fRunHeader;
}



Bool_t HarpoReaderRaw::Init(hReadMode mode)
{
  if(gHarpoDebug>0) Info("Init","Readout mode %d",mode);
  Pmm2Reader *lrdr;
  Bool_t res = true;
  fMode = mode;
  fDetSet->print();

  Long64_t pmm2Sync;
  if ( ! gHConfig->Lookup("Reader.pmm2Sync",pmm2Sync) )
    Info("Constructor","Use default pmm2Sync %i",fPmm2Sync);
  else
    fPmm2Sync = pmm2Sync;

  fLatency = 1268;
  Long64_t latency;
  if ( ! gHConfig->Lookup("Reader.latency",latency) )
    Info("Init","Use default Latency %d",fLatency);
  else
    fLatency = latency;

  ULong_t ndet;
  Long64_t trdr;
  for(ndet=0;ndet<gkNDetectors;ndet++) {
    //std::cout << "Init det " << ndet << std::endl;
    if (fDetSet->isExist(ndet)) {
      trdr = gHConfig->GetDetRdrType(static_cast<HarpoDetType>(ndet)); // Det Reader instance type
      switch(ndet){
      case XDCC:
        if (trdr == 0 ) {
          std::cout << "Det is T2K XDCC" << std::endl;
          fDetReader[ndet] = new HarpoDccReader(gHConfig);
        }
        else if (trdr == 1) {
          std::cout << "Det is MINOS XDCC" << std::endl;
          fDetReader[ndet] = new HarpoFeminosReader(gHConfig);
        }
        else {
          std::cout << "UNKNOWN XDCC" << std::endl;
          fDetReader[ndet] = NULL;
        }
        break;
      case YDCC:
        if (trdr == 0 ) {
          std::cout << "Det Reader is T2K YDCC" << std::endl;
          fDetReader[ndet] = new HarpoDccReader(gHConfig);
        }
        else if (trdr == 1) {
          std::cout << "Det Reader is MINOS YDCC" << std::endl;
          fDetReader[ndet] = new HarpoFeminosReader(gHConfig);
        }
        else {
          std::cout << "UNKNOWN YDCC Reader" << std::endl;
          fDetReader[ndet] = NULL;
        }
        break;
      case PMM2:
        if (trdr == 0 ) {
          std::cout << "Det Reader is PMM2 TEXT" << std::endl;
          // lrdr = new Pmm2Reader(gHConfig);
          // lrdr->SetReadOutMode(pTEXT);
          // fDetReader[ndet] = lrdr; 
        }
        else if (trdr == 1) {
          std::cout << "Det Reader is PMM2 Windows" << std::endl;
          // lrdr = new Pmm2Reader(gHConfig);
          // lrdr->SetReadOutMode(pWIN);
          // fDetReader[ndet] = lrdr; 
        }
        else if (trdr == 2) {
          std::cout << "Det Reader is PMM2 Linux" << std::endl;
          // lrdr = new Pmm2Reader(gHConfig);
          // lrdr->SetReadOutMode(pBIN);
          // fDetReader[ndet] = lrdr; 
        }
        else if (trdr == 3) {
          std::cout << "Det Reader is PMM2 Windows 2" << std::endl;
          // lrdr = new Pmm2Reader(gHConfig);
          // lrdr->SetReadOutMode(pWIN2);
          // fDetReader[ndet] = lrdr; 
        }
        // else if (trdr == 4) {
        //   std::cout << "Det Reader is PMM2 Linux 2" << std::endl;
        //   lrdr = new Pmm2Reader(gHConfig);
        //   lrdr->SetReadOutMode(pBIN2);
        //   fDetReader[ndet] = lrdr; 
        // }
        else {
          std::cout << "UNKNOWN PMM2 Reader" << std::endl;
          return 0;
          //      fDetReader[ndet] = NULL;
        }
        lrdr = new Pmm2Reader(gHConfig);
        lrdr->SetReadOutMode(trdr);
        fDetReader[ndet] = lrdr; 
        break;
      case TEMP_PRESSION:
        std::cout << "Det is TEMP_PRESSION" << std::endl;
        fDetReader[ndet] = NULL;
        break;
      case SIMDET:
        std::cout << "Det is SIM" << std::endl;
        fDetReader[ndet] = new HarpoSimReader(gHConfig);
        //std::cout << "DONE" << std::endl;
        break;
      }
      if(fDetReader[ndet] == NULL) continue;

      fDetReader[ndet]->SetDet(ndet);
      fDetReader[ndet]->SetRunNo(fRunNo);
      fDetReader[ndet]->SetDetFile(gHConfig->GetInFile(ndet));
      res = res && fDetReader[ndet]->Init();
      std::cout << res << std::endl;
    }
  }

  // Save det mask
  fRunHeader->SetSimDataFlag(false);
  fRunHeader->SetDetMask(fDetSet->GetMask());

  //  ULong_t ndet;
  HarpoDetHeader *rawhdr;
  for(ndet=0;ndet<gkNDetectors;ndet++) {
    if (fDetSet->isExist(ndet)) {
      std::cout << "Get Det Header " << ndet << std::endl;
      rawhdr = fDetReader[ndet]->GetDetHeader();
      fRunHeader->SetDetHeader(rawhdr,ndet);
    }
  }

  FillRunConfigInfo(); // Fill from DB

  return res;
}

Long_t      HarpoReaderRaw::Loop(Long_t maxevents)
{
  Long_t cnt=0;
  if(gHarpoDebug>0) Info("Loop", "maxevents = %ld",maxevents);

  std::cout << "HarpoReaderRaw::Loop" << std::endl;

  while (NextEvent()) {
    //GetEvent();
    // Pre
    Int_t st = hAnaOK;
    fEvent->GetHeader()->SetEvAnaStatus(st);

    if(fNevts<fFirstEvent || fNevts%(fSkipEvents+1)!=0){
      fEvent->Delete();
      continue;
    }

    for(Int_t i = 0; i<fNanalyses; i++) {
      fAnalyse[i]->ProcessEvent(fEvent);
      st = fEvent->GetHeader()->GetEvAnaStatus();
      if ( st < hAnaNextEvent ) break;
    }

    //    std::cout << "-- Raw Loop " << fNevts << " " << st << std::endl;
    fEvent->Delete();

    if (st == hAnaNextFile )
      {
        //how to go to next file ?
        break;
      }
    else if (st == hAnaStop)
      {
        break;
      }
    else if (st <= hAnaAbort)
      {
        Error(__FUNCTION__,
              "Fatal Analysis Error %d aborting prgram.\n",st);
        std::abort();
      }
    
    if (++cnt >= maxevents) break;
  }
  return cnt;  
}

Long_t      HarpoReaderRaw::Find(Long_t eventNo)
{
  if(gHarpoDebug>0) Info("Find","Event %ld",eventNo);
  Long_t cnt=0;
  Long_t maxevents = fMaxEvents;
  Bool_t found = false;
  HarpoEvent *evt;

  while (NextEvent()) 
    {
      evt = GetEvent();  // We need this as we assemble full event in GetEvent
      Info("Find","%ld => %ld",eventNo,(evt->GetHeader())->GetEvtNo());
      if ((evt->GetHeader())->GetEvtNo() == eventNo) {
        found = true;
        break;
      }
      SkipEvent();
      if (++cnt >= maxevents) break;
    }
  return found;

}

void      HarpoReaderRaw::SkipEvent()
{
  if(gHarpoDebug>0) Info("SkipEvent","event %ld",fNevts);
  ULong_t ndet;
  for(ndet=0;ndet<gkNDetectors;ndet++) {
    if (fDetSet->isExist(ndet)) {
      fDetReader[ndet]->SkipEvent();
    }
  }
}

void HarpoReaderRaw::FillRunConfigInfo()
{
  HarpoRunHeader * rh = GetRunHeader();
  Long_t run = rh->GetRun();
  HarpoRunConfig *info = new HarpoRunConfig();

  if ( ! info->GetRunConfig(run,gHConfig->GetDbUrl()) )
    {
      Info("FillRunConfig","No RunConfig in DB ,run %ld",run);
      return;
    }

  if (info->Lookup("Angle")) {
    TString val,type;
    if (info->GetVal("Angle",val)) {
      //  std::cout << "Angle val " << val << std::endl;
      rh->SetAngle(std::atoi(val.Data()));
    } else {
      std::cout << "Angle value not found" << std::endl;
    }
      
    // if (info->GetType("Angle",type)) {
    //   std::cout << "Angle type " << type << std::endl;
    // } else {
    //   std::cout<< "Angle type not found" << std::endl;
    // }
  }
  else
    std::cout << "Angle not found" << std::endl;

  if (info->Lookup("Ebeam")) {
    TString val,type;
    if (info->GetVal("Ebeam",val)) {
      //  std::cout << "Ebeam val " << val << std::endl;
      rh->SetEbeam(std::atof(val.Data()));
    } else {
      std::cout << "Ebeam value not found" << std::endl;
    }
      
    // if (info->GetType("Ebeam",type)) {
    //   std::cout << "Ebeam type " << type << std::endl;
    // } else {
    //   std::cout<< "Ebeam type not found" << std::endl;
    // }
  }
  else
    std::cout << "Ebeam not found" << std::endl;

  if (info->Lookup("Ephoton")) {
    TString val,type;
    if (info->GetVal("Ephoton",val)) {
      //  std::cout << "Ephoton val " << val << std::endl;
      rh->SetEgamma(std::atof(val.Data()));
    } else {
      std::cout << "Ephoton value not found" << std::endl;
    }
      
    // if (info->GetType("Ephoton",type)) {
    //   std::cout << "Ephoton type " << type << std::endl;
    // } else {
    //   std::cout<< "Ephoton type not found" << std::endl;
    // }
  }
  else
    std::cout << "Ephoton not found" << std::endl;

  //void SetLaser(Int_t value)    {fLaser = value;}
  // void SetTrigger(Int_t value)  {fTrigger = value;}
  // void SetIlaser(Float_t value) {fIlaser = value;}
  // void SetPola(Int_t value)     {fPola = value;}

  rh->SetRunConfig(info);

}