HarpoReaderRoot.cxx

Go to the documentation of this file.

// 
// File HarpoReaderRoot.cxx
//
#include "HarpoReaderRoot.h"
#include "HarpoDetReader.h"
#include "HarpoConfig.h"
#include "HarpoRunHeader.h"
#include "HarpoDccHeader.h"
#include "HarpoFeminosHeader.h"
#include "HarpoEvent.h"
#include "HarpoRecoEvent.h"
#include "HarpoDccEvent.h"
#include "HarpoFeminosEvent.h"
#include "Pmm2Event.h"
#include "HarpoSimEvent.h"
#include "HarpoTimeStamp.h"

#include "HarpoDccReader.h"
#include "HarpoSimReader.h"

#include "TClonesArray.h"
#include "TBranch.h"
#include "TSystemDirectory.h"

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


ClassImp(HarpoReaderRoot)


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

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

  fSync = false;
  fMode = hSyncTime;
}


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

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

  fSync = false;
  fMode = hSyncTime;
}

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

void HarpoReaderRoot::SetDetFile(TString * /* name */ , Int_t /* det */)
{
  //  if(det==0)
  //  fInputFileName = name; //???
}

 void HarpoReaderRoot::SetRootFile(TString *name)
 {
   //if(gHarpoDebug>0) 
   Info("Dummy SetRootFile","%s",name->Data());
 }
  
Bool_t HarpoReaderRoot::NextEvent()
{
  if(gHarpoDebug>1) Info("NextEvent","%ld",fNevts);
  Bool_t res = true;

  if(gHarpoDebug>1){
    for (UInt_t ndet = 0; ndet < gkNDetectors; ndet++)
      Info("NextEvent","*Det %d => %p",ndet,fEvents[ndet]);
    Info("NextEvent","*Reco => %p",fRecoEvent);
  }

  if(gHarpoDebug>1)
    Info("NextEvent","GetEvent %ld",fNevts);
  if(!fInputTree->GetEvent(fNevts))
    return false;

  if(gHarpoDebug>1){
    for (UInt_t ndet = 0; ndet < gkNDetectors; ndet++)
      Info("NextEvent","Det %d => %p",ndet,fEvents[ndet]);
    Info("NextEvent","Reco   => %p",fRecoEvent);
  }

  fNevts++;

  return res;
}

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

  fEvent = new HarpoEvent();
  // /* if(fEventHeader==NULL) */ 
  if(!fEventHeader) {
    Warning("GetEvent","Event Header not available");
    fEventHeader = new HarpoEventHeader();
  }//fEventHeader->SetEvtNo(fNevts);

  for(ndet=0;ndet<gkNDetectors;ndet++) {
    if (fDetSet->isExist(ndet)) {
      revt = fEvents[ndet];
      if(gHarpoDebug>1)
        Info("GetEvent","Det %ld exists, %p",ndet,revt);
      //      std::cout << ndet << ": evt " << revt << std::endl;
      fEvent->SetDetEvent(revt,ndet);
    }
  }

  fEvent->SetRecoEvent(fRecoEvent);
  fEvent->SetHeader(fEventHeader);

  for( UInt_t  ndet = 0; ndet < 2; ndet++ ) {
    HarpoFeminosHeader * fh = dynamic_cast<HarpoFeminosHeader *>(fRunHeader->GetHeader(ndet));
    if (fh != NULL) {
      //      std::cout << "Fem Peds " << ndet << std::endl;
      HarpoDccMap *dm = (fEvents[ndet])->GetMap();
      dm->SetPedestals(fh->GetPedestals());
      //fh->GetPedestals()->Print();
      continue;
    }
    HarpoDccHeader * dh = dynamic_cast<HarpoDccHeader *>(fRunHeader->GetHeader(ndet));
    if (dh != NULL) {
      //      std::cout << "Fem Peds " << ndet << std::endl;
      HarpoDccMap *dm = (fEvents[ndet])->GetMap();
      dm->SetPedestals(dh->GetPedestals());
      //fh->GetPedestals()->Print();
      continue;
    }
    Warning("GetEvent","Unknown Header %d",ndet);
  }



  return fEvent;
}

HarpoRunHeader *HarpoReaderRoot::GetRunHeader()
{
  // fRunHeader = fInRootFile->Get("fRunHeader");
  // return fRunHeader;
  if(gHarpoDebug>0) Info("GetRunHeader"," ");
  return fRunHeader;
}


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

  if(gHarpoDebug>0) Info("Init","Open ROOT file %s",gHConfig->GetInFile(0)->Data());
  //fInRootFile = TFile::Open(gHConfig->GetInFile(0)->Data());
  if(gHarpoDebug>0) Info("Init","File open");
  fEventHeader = new HarpoEventHeader();
  if(gHarpoDebug>0) Info("Init","Get TTree HarpoTree");
  fInputTree = new TChain("HarpoTree");
  if(gHConfig->GetInFile(0)->EndsWith(".root")){
    fInRootFile = TFile::Open(gHConfig->GetInFile(0)->Data());
    fInputTree->Add(gHConfig->GetInFile(0)->Data());
    Info("Init","Adding %s",gHConfig->GetInFile(0)->Data());
  }else{

    TString strtmp = *gHConfig->GetInFile(0);
    Ssiz_t pos = strtmp.Last('/');
    TString dirname = strtmp(0,pos);
    TString fileprefix = strtmp(pos+1,strtmp.Length());

    Info("Init","%s  /  %s",dirname.Data(), fileprefix.Data());

    Int_t nfiles = 0;
    TSystemDirectory dir("datadir", dirname.Data());
    TList *files = dir.GetListOfFiles();
    if (files) {
      TSystemFile *file;
      TString fname;
      TIter next(files);
      while ((file=(TSystemFile*)next())) {
         fname = file->GetName();
         if(!fname.BeginsWith(fileprefix.Data()))
            continue;
         if(!fname.EndsWith(".root"))
            continue;
         Info("Init","Include file %s",fname.Data());

         if(nfiles==0)
           fInRootFile = TFile::Open(Form("%s/%s",dirname.Data(),fname.Data()));
         nfiles++;
      }
    }
    if(nfiles==0){
      Warning("Init","No file matching the pattern %s",gHConfig->GetInFile(0)->Data());
      return false;
    }

    Int_t nFiles = fInputTree->Add(Form("%s*.root",gHConfig->GetInFile(0)->Data()));
    Info("Init","Adding %s (%d files)",Form("%s*.root",gHConfig->GetInFile(0)->Data()),nFiles);
  }
  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  Long64_t trdr[gkNDetectors];
  
  fRunHeader = NULL;
  fRunHeader = //new HarpoRunHeader();
    (HarpoRunHeader*)fInRootFile->Get("fRunHeader");

  if (fRunHeader == NULL ) 
      Warning(__FUNCTION__,"HarpoRunHeader not found\n");


  // Check mask
  ULong_t gmask = fDetSet->GetMask();  // global processing mask
  ULong_t hmask = fRunHeader->GetDetMask(); // maks for data in root file
  if ( gmask != hmask ) {
    Info(__FUNCTION__,"Processing mask 0x%lx <> Header mask 0x%lx.\n",gmask,hmask);

    ULong_t diff_mask = ( gmask ^ hmask );
    ULong_t skip_mask = ( hmask & diff_mask ) ; // det to skip
    ULong_t nodet_mask = ( gmask & diff_mask ) ; // no det in root files

    if (skip_mask)
      Info(__FUNCTION__,"Processing skiping det 0x%lx.\n",skip_mask);
    if (nodet_mask)
      Warning(__FUNCTION__,"Defectors  0x%lx not in ROOT file!\n",nodet_mask);

    fDetSet->SetMask(gmask&hmask);
    fRunHeader->SetDetMask(gmask&hmask);

  }
  HarpoDetHeader *rawhdr;
  // HarpoDccHeader *drawhdr;
  // HarpoDeminosHeader *frawhdr;
  ULong_t ndet;
  for(ndet=0;ndet<gkNDetectors;ndet++) {
    if (fDetSet->isExist(ndet)) {
      std::cout << "Get Det Header " << ndet << std::endl;
      rawhdr = (HarpoDetHeader*)fInRootFile->Get(Form("fDetHeader%ld",ndet));
      if(rawhdr !=0 ) {
          if ( ndet < 2 ) {
            if (dynamic_cast<HarpoDccHeader *>(rawhdr) != NULL ) {
              trdr[ndet] = 0;
              gHConfig->SetDetRdrType(static_cast<HarpoDetType>(ndet),0);
              std::cout << "Det " << ndet << " is T2K DCC" << std::endl;
            } else if (dynamic_cast<HarpoFeminosHeader *>(rawhdr) != NULL ) {
              trdr[ndet] = 1;
              gHConfig->SetDetRdrType(static_cast<HarpoDetType>(ndet),1);
              std::cout << "Det " << ndet << " is MINOS DCC" << std::endl;
            } else {
              trdr[ndet] = -1; //  ??Reasonable value??
              gHConfig->SetDetRdrType(static_cast<HarpoDetType>(ndet),-1);
              std::cout << "Det " << ndet << " not a DCC" << std::endl;
            }
          }
        fRunHeader->SetDetHeader(rawhdr,ndet);
        std::cout << "Header Class Name " << rawhdr->ClassName() << std::endl;
      }
    }
  }
  //------------------------------------------------------------------------

  //  fInputTree = (TTree*)fInRootFile->Get("HarpoTree");
  if(gHarpoDebug>0) Info("Init","Initialise event objects");
  //  Long64_t trdr;
  for (UInt_t ndet = 0; ndet < gkNDetectors; ndet++){
    if (gHDetSet->isExist(ndet)){
      //      trdr = gHConfig->GetDetRdrType(static_cast<HarpoDetType>(ndet)); // Det Reader instance type
      HarpoDet dt(ndet);  
      // TBranch * b = fInputTree->GetBranch(dt.GetName());
      // const char * className = b->GetClassName();
      // printf("Branch %s Class %s\n",dt.GetName(),className);
      switch(ndet){
      case XDCC:
        if (trdr[ndet] == 0 ) {
          //std::cout << "Det is T2K XDCC" << std::endl;
          fEvents[ndet] = new HarpoDccEvent();
        }
        else if (trdr[ndet] == 1) {
          //std::cout << "Det is MINOS XDCC" << std::endl;
          fEvents[ndet] = new HarpoFeminosEvent();
        }
        else {
          std::cout << "UNKNOWN XDCC" << std::endl;
          fDetReader[ndet] = NULL;
        }
        break;
      case YDCC:
        if (trdr[ndet] == 0 ) {
          //std::cout << "Det is T2K YDCC" << std::endl;
          fEvents[ndet] = new HarpoDccEvent();
        }
        else if (trdr[ndet] == 1) {
          //std::cout << "Det is MINOS YDCC" << std::endl;
          fEvents[ndet] = new HarpoFeminosEvent();
        }
        else {
          std::cout << "UNKNOWN YDCC" << std::endl;
          fDetReader[ndet] = NULL;
        }
        break;
      case PMM2:
        fEvents[ndet] = new Pmm2Event();
        break;
        // case TEMP_PRESSION,:
        //       break;
      case SIMDET:
        fEvents[ndet] = new HarpoSimEvent();
        break;
      }
      if(gHarpoDebug>0) Info("Init","Set Branch %d (%s)",ndet,dt.GetName());
      fInputTree->SetBranchAddress(dt.GetName(),&fEvents[ndet]);
    }
  }
  if(gHarpoDebug>0) Info("Init","Event objects initialised");
  fInputTree->SetBranchAddress("BRawEventHeader",&fEventHeader);

  fRecoEvent = new HarpoRecoEvent();
  if(gHarpoDebug>0) Info("Init","Set Branch Reco");
  fInputTree->SetBranchAddress("RecoEvent",&fRecoEvent);
  if(gHarpoDebug>0) Info("Init","Reco objects initialised");

  return res;
}

Long_t      HarpoReaderRoot::Loop(Long_t maxevents)
{
  fNevts = gHConfig->GetFromEvent();

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

  Int_t nEntries  = fInputTree->GetEntries();
  if(nEntries<maxevents) maxevents=nEntries;

  while (NextEvent()) {
    //for(Int_t ev = 0; ev<maxevents; ev++){
   if(gHarpoDebug>1) Info("Loop", "########### event = %ld ############",fNevts);
    GetEvent();

    // Pre
    Int_t st = hAnaOK;
    fEvent->GetHeader()->SetEvAnaStatus(st);

    for(Int_t i = 0; i<fNanalyses; i++){
      if(gHarpoDebug>1) Info("Loop", "Start Analysis %d",i);
      fAnalyse[i]->ProcessEvent(fEvent);
      st = fEvent->GetHeader()->GetEvAnaStatus();
      if ( st < hAnaNextEvent ) break;
    }

   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      HarpoReaderRoot::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;

  fNevts = eventNo+1;

  return fInputTree->GetEvent(eventNo);
}

// void      HarpoReaderRoot::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();
//     }
//   }
// }