HarpoAnalyseWriter.cxx

Go to the documentation of this file.

// 
// File HarpoAnalyseWriter.cxx
//
#include "HarpoAnalyseWriter.h"
#include "HarpoConfig.h"
#include "HarpoDetSet.h"
//#include "HarpoDetSet.h"
#include "HarpoEvent.h"
#include "HarpoDccEvent.h"
#include "HarpoFeminosEvent.h"
#include "Pmm2Event.h"
#include "HarpoSimEvent.h"
#include "HarpoRecoEvent.h"

#include "TFile.h"
#include "TROOT.h"

#include <cstdlib>
#include <cstring>
#include <cassert>
#include <iostream>

ClassImp(HarpoAnalyseWriter)

const int bufsize = 16000;
const int splitlevel = 2;
const int compress = 6;   //default ???

HarpoAnalyseWriter::HarpoAnalyseWriter()
 {
  rootFile = NULL;
  fNewFile = NULL;
  fNewTree = NULL;
  fRawBranch = NULL;
  fNewNb = 0;
  fDoAna = true;

  fEvtHdr = NULL ; // event header

  ULong_t ndet;
  for (ndet = 0; ndet < gkNDetectors; ndet++) 
    fDetEvt[ndet] = NULL;
 }

void   HarpoAnalyseWriter::print()
 {
   //   if (! fDoAna ) return;
   ULong_t ndet; // number of detectors
   HarpoEventHeader *hdr = fEvt->GetHeader();
   HarpoDetEvent *devt;
   
   assert(hdr != NULL);
   hdr->print();
   
   for (ndet = 0; ndet < gkNDetectors; ndet++) {
     if (gHDetSet->isExist(ndet)) {
       std::cout << "Print Det " << ndet << std::endl;
       devt = fEvt->GetDetEvent(ndet);
       devt->print();
       std::cout << "** Class ["<< devt->ClassName() << "]" << std::endl;
     }
   }
 }

// Set pointer to new event for analyse, anlyse event , and delete
void   HarpoAnalyseWriter::ProcessEvent(HarpoEvent *event)
{
  //  if ( fDoAna && (event->GetHeader()->GetEvAnaStatus() < hAnaEventWrt) )
  if ( (event->GetHeader()->GetEvAnaStatus() < hAnaEventWrt) )
    return;

  Set(event);
  process();
}

void   HarpoAnalyseWriter::process()
{
  
  if(gHarpoDebug>1)
    std::cout << "HarpoRunHeader = " <<  fRunHeader << std::endl;
   //HarpoDetEvent *devt;
  nEvents++;
  if(gHarpoDebug>1){
    Info("process()","Processing Event %ld",nEvents);
    fEvt->print();
  }
  
  //  if (! fDoAna ) return;
  // ULong_t ndet;
  
  if(nEvents%100 == 0)
    Info("process()"," Processing Event %ld",nEvents);
  
  fEvtHdr = fEvt->GetHeader();
  (void) fEvt->GetDccMap(0); // create map 
  (void) fEvt->GetDccMap(1); //
  for (UInt_t ndet = 0; ndet < gkNDetectors; ndet++){
    if (gHDetSet->isExist(ndet)){
      if(gHarpoDebug>1)
        Info("process()"," Get DetEvent %d",ndet);
      fDetEvt[ndet] = fEvt->GetDetEvent(ndet);
      if(!fDoAna && ndet == SIMDET)
        ((HarpoSimEvent*)fDetEvt[ndet])->DeleteIonisationTracks();
      //       fDetEvt[1] = fEvt->GetDetEvent(1);
    }
  }
  //   Info("process","HasRecoEvent: %d",fEvt->HasRecoEvent());
  
  if(fEvt->HasRecoEvent())
    fRecoEvent = fEvt->GetRecoEvent();
  fNewNb =+ fNewTree->Fill();
  if(gHarpoDebug>0) Info("process","fNewNb = %d",fNewNb);
  if(gHarpoDebug>1)
    fRecoEvent->Summary();
  if(gHarpoDebug>0)
    Info("process()"," Tree Filled ");
  
  
  
  
  if(gHConfig->GetSplit() != 0 && nEvents%gHConfig->GetSplit()==0){
    Info("process","Splitting (evt %li)",nEvents);
    Save();
    //     InitFile(Form("%s_%d.root",rootFile->Data(),fNfiles));
    InitFile(Form("%s_%d.root",rootFile->Data(),fNfiles));
    fNfiles++;
  }
  
  
  //gROOT->ls();
}



 void   HarpoAnalyseWriter::Init()
 {
   fDoAna = gHConfig->isWriterEnabled();

   if(!fDoAna){
     Info("HarpoAnalyseWriter","Do Not Save ionisation data");  
     //     Info("HarpoAnalyseWriter","Do Not run analysis");  
     //return;
   }

   rootFile = gHConfig->GetOutFile();
   if ( rootFile == NULL ) {
     std::cout << gHConfig->GetProgramName() << " " <<
       "No Root Output File name gived, use default" << std::endl;
     rootFile = new TString("hraw");
   }

   fNfiles = 0;
   if(rootFile->EndsWith(".root")){
     Ssiz_t pos = rootFile->Last('.');
     rootFile->Remove(pos,5);
   }
   if(gHConfig->GetSplit()>0)
     InitFile(Form("%s_0.root",rootFile->Data()));
   else
     InitFile(Form("%s.root",rootFile->Data()));
   fNfiles++;
 }

void HarpoAnalyseWriter::InitFile(const char* fileName)
{

  //   if (! fDoAna ) return;
   // Dummy Event
   fEvt = new HarpoEvent();
   fEvtHdr = new HarpoEventHeader();
   fEvt->SetHeader(fEvtHdr);
    
   fNewFile = new TFile(fileName,"RECREATE",
                        "Harpo ROOT file");
   fNewFile->SetCompressionLevel(compress);
   fNewTree = new TTree("HarpoTree","Harpo ROOT Tree");
   fNewTree->SetAutoSave();  // autosave when 1 Gbyte written
   fNewTree->SetAutoFlush();  // autosave when 1 Gbyte written
   //TTree::SetBranchStyle(branchStyle);

   // fRawBranch = fNewTree->Branch("BRawEvent","HarpoEvent",
   //                                 &fEvt, bufsize,splitlevel);
   fRawBranch = fNewTree->Branch("BRawEventHeader","HarpoEventHeader",
                                   &fEvtHdr, 32000,99);
   fRecoEvent = new HarpoRecoEvent();
   fNewTree->Branch("RecoEvent","HarpoRecoEvent",
                    &fRecoEvent, 32000,99);
   gHDetSet->print();
   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);
       if(gHarpoDebug>0) Info("Init()","Add Branch %d (%s)",ndet,dt.GetName());
       switch(ndet){
       case XDCC:
        if (trdr == 0 ) {
          std::cout << "Branch for T2K XDCC" << std::endl;
          fDetEvt[ndet] = new HarpoDccEvent();
          fRawBranch = fNewTree->Branch(dt.GetName(),"HarpoDccEvent",
                                        &fDetEvt[ndet], 16000,99);
        }
        else if (trdr == 1) {
          std::cout << "Branch for MINOS XDCC" << std::endl;
          fDetEvt[ndet] = new HarpoFeminosEvent();
          fRawBranch = fNewTree->Branch(dt.GetName(),"HarpoFeminosEvent",
                                        &fDetEvt[ndet], 16000,99);
        }
        else {
          std::cout << "UNKNOWN XDCC" << std::endl;
        }
         break;
       case YDCC:
        if (trdr == 0 ) {
          std::cout << "Branch for T2K YDCC" << std::endl;
          fDetEvt[ndet] = new HarpoDccEvent();
          fRawBranch = fNewTree->Branch(dt.GetName(),"HarpoDccEvent",
                                        &fDetEvt[ndet], 16000,99);
        }
        else if (trdr == 1) {
          std::cout << "Branch for MINOS YDCC" << std::endl;
          fDetEvt[ndet] = new HarpoFeminosEvent();
          fRawBranch = fNewTree->Branch(dt.GetName(),"HarpoFeminosEvent",
                                        &fDetEvt[ndet], 16000,99);
        }
        else {
          std::cout << "UNKNOWN YDCC" << std::endl;
        }
         break;
       case PMM2:
         fDetEvt[ndet] = new Pmm2Event();
         fRawBranch = fNewTree->Branch(dt.GetName(),"Pmm2Event",
                                       &fDetEvt[ndet], 16000,99);
         break;
         // case TEMP_PRESSION,:
         //      fRawBranch = fNewTree->Branch(dt.GetName(),"HarpoDccEvent",
         //                                    &fDetEvt[ndet], 16000,99);
         //      break;
       case SIMDET:
         fDetEvt[ndet] = new HarpoSimEvent();
         fRawBranch = fNewTree->Branch(dt.GetName(),"HarpoSimEvent",
                                       &fDetEvt[ndet], 16000,99);
         break;
       }
     }
     
   }
 }

void   HarpoAnalyseWriter::Save(char * /* mode */)
{
  //  if(!fDoAna) return;

  Info("Save","%s",fNewFile->GetName());

  fRunHeader->print();
  fNewFile->cd();
  
  for (UInt_t ndet = 0; ndet < gkNDetectors; ndet++){
    if (gHDetSet->isExist(ndet)){
      HarpoDetHeader* hdr = fRunHeader->GetHeader(ndet);
      std::cout << "Save Header " << ndet << " " << hdr << std::endl;
      if(hdr) hdr->Write(Form("fDetHeader%d",ndet));
    }
  }
  std::cout << "fRunHeader " << fRunHeader << std::endl;
  fRunHeader->Write("fRunHeader");
  
  
  fNewFile->Write();
  fNewFile->Close();
  printf("fNewFile %s closed %d bytes written\n", fNewFile->GetName(),fNewNb );
}