HarpoSimReader.cxx

Go to the documentation of this file.

// 
// File HarpoSimReader.cxx
//
#include "HarpoSimReader.h"
#include "HarpoConfig.h"
#include "HarpoSimHeader.h"
#include "HarpoDccEvent.h"
#include "HarpoDccHeader.h"
#include "HarpoFeminosEvent.h"
#include "HarpoFeminosHeader.h"
#include "Pmm2Event.h"
#include "Pmm2Header.h"
#include "TpcSimMCEvent.h"
//#include "HarpoPTEvent.h"
#include "HarpoSimEvent.h"
#include "HarpoSimElectronics.h"
#include "HarpoDet.h"
#include <TSystem.h>
#include <TRandom3.h>
//#include <TTree.h>

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

ClassImp(HarpoSimReader)

#define EXPECTED_NB_SIM_FRAGMENTS 18
#define NB_FEM_FRAGMENTS  72
#define MAX_FEM_DATA_BUFFER 2097152


HarpoSimReader::HarpoSimReader() : HarpoDetReader()
{

  std::cout << "Create SIM Reader"  << std::endl;

  //  fCfg = new HarpoConfig();

  //fRunNo = 0;
  fDone = false;
  fError = false;
  //  fPlane = -1;

  fIonType = 1;
  Long64_t iontype;
  if ( ! gHConfig->Lookup("Sim.IonisationType",iontype) )
    Info("HarpoSimReader","Use default simulation for ionisation");
  else
    fIonType = iontype;

  Info("HarpoSimReader","fIontype = %d", fIonType);

  //fInFile = (TString *) NULL;
  fIn = (FILE *) NULL;
  //fNevts = 0;

  //fNewFile = (TFile *) NULL;
  fNewNb = 0;

  fEvents = new HarpoDetEvent *[gkNDetectors];

  //fEvent = (HarpoSimEvent *) NULL;
  //  fRunHeader = (HarpoSimHeader *) NULL;
  fIonisationChamber = (TpcSimIonisationChamber *) NULL;
  fDriftChamber = (TpcSimDriftChamber *) NULL;
  fHarpoChamber = (HarpoSimChamber *) NULL;
  fEventGenerator = (TpcSimEventGenerator *) NULL;
  fHarpoReadout = (HarpoSimReadout *) NULL;
  fDigi = (HarpoSimElectronics *) NULL;

  fSimAnalyse = (HarpoSimAnalyse *) NULL;
  fMaxEvents = LONG_MAX;
  //fPrevTimeStamp = 0;
  isEventSaved = false;

  fTimeStamp = 0;

  //  fDoSave = kTRUE;
}


HarpoSimReader::HarpoSimReader(HarpoConfig *cfg) : HarpoDetReader(cfg)
{
  // Config Function
  fDone = false;
  fError = false;

  fCfg = cfg;

  fIonType = 1;

  Long64_t iontype;
  if ( ! gHConfig->Lookup("Sim.IonisationType",iontype) )
    Info("HarpoSimReader","Use default simulation for ionisation");
  else
    fIonType = iontype;

  //fRunNo = cfg->GetRunNo();
  //fInFile = cfg->GetInFile();
  fIn = (FILE *) NULL;
  //fNevts = 0;
  //fMaxEvents = cfg->GetMaxEvents();
  if(gHarpoDebug>0) Info("HarpoSimReader","fMaxEvents = %li",fMaxEvents);

  TString *FileName = cfg->GetOutFile();
  //  Info("HarpoSimReader","Filename = %s",FileName->Data());
  fNewNb = 0;

  fEvent = (HarpoSimEvent *) NULL;
  fEvents = new HarpoDetEvent *[gkNDetectors];

  //  fRunHeader = (HarpoSimHeader *) NULL;

  fSimAnalyse = (HarpoSimAnalyse *) NULL;

  fIonisationChamber = (TpcSimIonisationChamber *) NULL;
  fDriftChamber = (TpcSimDriftChamber *) NULL;
  fHarpoChamber = (HarpoSimChamber *) NULL;
  fEventGenerator = (TpcSimEventGenerator *) NULL;
  fHarpoReadout = (HarpoSimReadout *) NULL;
  fDigi = (HarpoSimElectronics *) NULL;
 
  if ( FileName != (TString *) NULL )
    SetRootFile(FileName);

  //fPrevTimeStamp = 0;
  isEventSaved = false;

  // fDoSave = kTRUE;
}

HarpoSimReader::~HarpoSimReader()
{
  //  std::fclose(fIn);

  if(fIonType)
    fStr.close();

  // if(fSave){
  //   fTree->Write();
  //   fNewFile->Close();
  // }
  //  delete fRunHeader;
}

void HarpoSimReader::SetSimFile(TString *name)
{
  fInFile = name;
}

void HarpoSimReader::SetRootFile(TString *name)
{

  if(gHarpoDebug>0)  Info("SetRootFile","%s",name->Data());
  fNewFileName = name; 
  fSave = kTRUE;
  fNewFile = new TFile(name->Data(),"RECREATE","HARPO SIM Events");
}

void HarpoSimReader::SetAnalyseFunction(HarpoSimAnalyse *func)
{
  fSimAnalyse = func;
}
  
Bool_t HarpoSimReader::NextEvent()
{

  if(gHarpoDebug>1) Info("NextEvent","Generating event %ld",fNevts);
  // HarpoDccEvent* xdccEvent = new HarpoDccEvent();
  // HarpoDccEvent* ydccEvent = new HarpoDccEvent();
  HarpoDetEvent* xdccEvent;// = new HarpoDccEvent();
  HarpoDetEvent* ydccEvent;// = new HarpoDccEvent();



  Pmm2Event* pmm2Event = new Pmm2Event();
  HarpoSimEvent* simEvent = new HarpoSimEvent();
  fHarpoChamber->StartEvent();

  Bool_t isDone = GenerateEvent(simEvent);
  if(!isDone) return isDone;

  fHarpoChamber->FinishEvent();
  HarpoDccMap* mapX = fDigi->GetDccMapX();
  Long64_t trdr = gHConfig->GetDetRdrType(static_cast<HarpoDetType>(0));
  if (trdr == 0 ) {
    //    std::cout << "Det is T2K XDCC" << std::endl;
    xdccEvent = new HarpoDccEvent();
    ((HarpoDccEvent*)xdccEvent)->SetMap(mapX);
    ((HarpoDccEvent*)xdccEvent)->SetHeader(0,fNevts);
    ((HarpoDccEvent*)xdccEvent)->SetZeroSuppressed(kTRUE);
  }
  else if (trdr == 1) {
    //    std::cout << "Det is MINOS XDCC" << std::endl;
    xdccEvent = new HarpoFeminosEvent();
    ((HarpoFeminosEvent*)xdccEvent)->SetMap(mapX);
    ((HarpoFeminosEvent*)xdccEvent)->SetHeader(0,fNevts);
    ((HarpoFeminosEvent*)xdccEvent)->SetZeroSuppressed(kTRUE);
  }
  else {
    std::cout << "UNKNOWN XDCC" << std::endl;
    xdccEvent = NULL;
  }


  HarpoDccMap* mapY = fDigi->GetDccMapY();
  trdr = gHConfig->GetDetRdrType(static_cast<HarpoDetType>(1));
  if (trdr == 0 ) {
    //   std::cout << "Det is T2K YDCC" << std::endl;
    ydccEvent = new HarpoDccEvent();
    ((HarpoDccEvent*)ydccEvent)->SetMap(mapY);
    ((HarpoDccEvent*)ydccEvent)->SetHeader(0,fNevts);
    ((HarpoDccEvent*)ydccEvent)->SetZeroSuppressed(kTRUE);
  }
  else if (trdr == 1) {
    //  std::cout << "Det is MINOS YDCC" << std::endl;
    ydccEvent = new HarpoFeminosEvent();
    ((HarpoFeminosEvent*)ydccEvent)->SetMap(mapY);
    ((HarpoFeminosEvent*)ydccEvent)->SetHeader(0,fNevts);
    ((HarpoFeminosEvent*)ydccEvent)->SetZeroSuppressed(kTRUE);

  }
  else {
    std::cout << "UNKNOWN YDCC" << std::endl;
    ydccEvent = NULL;
  }



  Int_t nTracks = simEvent->GetMCEvent()->GetNtracks();
  if(gHarpoDebug>1) 
    Info("NextEvent","Ntracks = %d",nTracks);
  Int_t pm[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
  for(Int_t tr = 0; tr<nTracks; tr++){
    TpcSimMCTrack* track = simEvent->GetMCEvent()->GetTrack(tr);

    Double_t x[3];
    x[0] = track->GetX0();
    x[1] = track->GetY0();
    x[2] = track->GetZ0();
    Double_t p[3];
    p[0] = track->GetPx();
    p[1] = track->GetPy();
    p[2] = track->GetPz();
    
    if(gHarpoDebug>1) Info("NextEvent","%g %g %g    %g %g %g",x[0],x[1],x[2],p[0],p[1],p[2]);

    Double_t xup[3] = {150,150,300};
    Double_t xdown[3] = {-150,-150,0};
    
    Bool_t scint[6] = {1,1,1,1,1,1};
    
    for(Int_t j = 0; j<3; j++){
      for(Int_t i = 0; i<3; i++){
        if(i==j) continue;
        Double_t xtmp = x[i] + p[i]*(xup[j] + 50 - x[j])/p[j];
        if(gHarpoDebug>1) Info("NextEvent","   %d %d   %g",i,j,xtmp);
        scint[2*j] = scint[2*j] && (xtmp<xup[i] && xtmp>xdown[i]);
        xtmp = x[i] + p[i]*(xdown[j] - 50 - x[j])/p[j];
        //std::cout << "xtmp = " << xtmp << std::endl;
        scint[2*j+1] = scint[2*j+1] && (xtmp<xup[i] && xtmp>xdown[i]);
      }
    }

    // for(Int_t i = 0; i<6; i++){
    //   pm[2*i] = pm[2*i] || scint[i];
    //   pm[2*i+1] = pm[2*i+1] || scint[i];
    // }
    pm[0]  = pm[0]  || scint[3];
    pm[1]  = pm[1]  || scint[3];
    pm[2]  = pm[2]  || scint[4];
    pm[3]  = pm[3]  || scint[4];
    pm[4]  = pm[4]  || scint[5];
    pm[5]  = pm[5]  || scint[5];
    pm[6]  = pm[6]  || scint[1];
    pm[7]  = pm[7]  || scint[1];
    pm[8]  = pm[8]  || scint[2];
    pm[9]  = pm[9]  || scint[2];
    pm[10] = pm[10] || scint[0];
    pm[11] = pm[11] || scint[0];

  }


  Int_t size = 0;
  for(Int_t i = 0; i<12; i++){
    //    if(gRandom->Rndm()<0.5) continue;
    Int_t charge = 800*pm[i];
    //    Int_t charge = Int_t((7*pm[i]+0.)*gRandom->Rndm()*100);
    if(charge < 50) continue;
    Pmm2Mes mes;
    mes._ChannelNum = i+1;
    mes._Charge = charge;
    pmm2Event->AddMes(mes);
    size++;
  }
  pmm2Event->SetTriggerType(Int_t(gRandom->Rndm()*6));
  pmm2Event->SetHeader(size,fNevts);

  // std::cout << "xdccEvent " << xdccEvent << std::endl;
  // std::cout << "ydccEvent " << ydccEvent << std::endl;
  // std::cout << "pmm2Event " << pmm2Event << std::endl;

  // xdccEvent->SetTimeStamp(100*fNevts);
  // ydccEvent->SetTimeStamp(100*fNevts);
  // pmm2Event->SetTimeStamp(100*fNevts);
  // simEvent->SetTimeStamp(100*fNevts);
  fEvents[XDCC] = xdccEvent;
  fEvents[YDCC] = ydccEvent;
  fEvents[PMM2] = pmm2Event;
  fEvents[SIMDET] = simEvent;

  //cout << "Fill Tree" << endl;
  
  fNevts++;

  return isDone;
}


Bool_t HarpoSimReader::GenerateEvent(HarpoSimEvent* simEvent)
{
  if(gHarpoDebug>1) Info("GenerateEvent", "fIonType %d",fIonType);
  switch(fIonType){
  case 1:
    return GenerateEventBasic(simEvent);
  case 2:
    return GenerateEventFromText(simEvent, &fStr);
  default:
    return GenerateEventBasic(simEvent);
  }

  return true;
}


Bool_t HarpoSimReader::GenerateEventBasic(HarpoSimEvent* simEvent)
{
  TpcSimMCEvent* MCEvent = fEventGenerator->GenerateEvent();
  simEvent->SetMCEvent(MCEvent);
  Int_t nTr = MCEvent->GetNtracks();
  if(gHarpoDebug>1) Info("GenerateEventBasic","%d MC tracks generated",nTr);
  for(Int_t iTr = 0; iTr<nTr; iTr++){
    TpcSimMCTrack* mctr = MCEvent->GetTrack(iTr);
    TpcSimIonisationTrack* tr;
    Info("GenerateEvent", "fIonType %d",fIonType);
    tr = fIonisationChamber->GenerateTrack(mctr);
    simEvent->AddIonisationTrack(iTr,tr);
    fHarpoReadout->Reset();
    fHarpoChamber->ProcessTrack(tr);
 }

  return true;
}


Bool_t HarpoSimReader::GenerateEventFromText(HarpoSimEvent* simEvent, std::ifstream *str)
{
  TpcSimMCEvent* MCEvent = new TpcSimMCEvent();
  //if(!str->good()) return false;
  Float_t xg = 0, yg = 0, zg = 0, tg = 0, pgx = 0, pgy = 0, pgz = 0, eg = 0, polax = 0, polay = 0, polaz = 0;
  while(str->good()){
    Char_t line[256];
    str->getline(line,256);
    if(strstr(line,"NEW EVENT") != NULL){
      sscanf(line," NEW EVENT %f %f %f %f %f %f %f %f %f %f %f", &xg, &yg, &zg, &tg, &pgx, &pgy, &pgz, &eg, &polax, &polay, &polaz);
      if(gHarpoDebug>1) Info("GenerateEventFromText:","%s",line);
      break;
    }
    //    Info("GenerateEventFromText","No data left");
    //    return false;
  }
  if(!str->good()) return false;

  if(zg){
    Info("GenerateEventFromText","ADDING PHOTON TRACK (%g %g %g %g) Pola (%g %g %g)",pgx,pgy,pgz,eg,polax,polay,polaz);
    TpcSimMCTrack* mctrack = new TpcSimMCTrack(xg,yg,zg,tg,pgx,pgy,pgz,22);
    mctrack->SetE(eg);
    mctrack->SetPolax(polax);
    mctrack->SetPolay(polay);
    mctrack->SetPolaz(polaz);
    MCEvent->AddTrack(mctrack);
  }

  Int_t iTr = 0;
  simEvent->SetMCEvent(MCEvent);
  while(str->good()){
    TpcSimIonisationTrack* tr;
    tr = fIonisationChamber->GenerateTrackFromText(str);
    if(!tr) break;
    simEvent->AddIonisationTrack(iTr,tr);
    fHarpoReadout->Reset();
    fHarpoChamber->ProcessTrack(tr);
    iTr++;
    // if(tr)
    //tr->Delete();
  }

  return true;
}



HarpoSimEvent *HarpoSimReader::GetEvent()
{
  return NULL;
  //  return fEvent;
}

HarpoDetEvent *HarpoSimReader::GetEvent(ULong_t ndet)
{
  return fEvents[ndet];
}

HarpoDetHeader *HarpoSimReader::GetDetHeader(UInt_t ndet)
{
  return fRunHeader[ndet];
}


Bool_t HarpoSimReader::Init()
{
  fIonisationChamber = new TpcSimIonisationChamber();
  fEventGenerator = new TpcSimEventGenerator();

  fDriftChamber = new TpcSimDriftChamber();

  fHarpoChamber = new HarpoSimChamber();
  fHarpoReadout = new HarpoSimReadout();
  fDigi = new HarpoSimElectronics();
  fHarpoChamber->SetDriftChamber(fDriftChamber);
  fHarpoChamber->SetReadout(fHarpoReadout);
  fHarpoChamber->SetElectronics(fDigi);

  Long64_t trdr = gHConfig->GetDetRdrType(static_cast<HarpoDetType>(0));
  if (trdr == 0 ) {
    //   std::cout << "Det is T2K XDCC" << std::endl;
    fRunHeader[XDCC] = new HarpoDccHeader();
    // HarpoPedestal* pX = new HarpoPedestal();
    // pX->SetMean(0.);
    // pX->SetSigma(fDigi->GetSigmaNoise());
    fRunHeader[YDCC] = new HarpoDccHeader();
    // HarpoPedestal* pY = new HarpoPedestal();
    // pY->SetMean(0.);
    // pY->SetSigma(fDigi->GetSigmaNoise());
    // fRunHeader[1]->SetPedestal(pY);
  }
  else if (trdr == 1) {
    //    std::cout << "Det is MINOS XDCC" << std::endl;
    fRunHeader[XDCC] = new HarpoFeminosHeader();
    ((HarpoFeminosHeader*)fRunHeader[XDCC])->SetPedestals(fDigi->GetPedestals(XDCC));
    fRunHeader[YDCC] = new HarpoFeminosHeader();
    ((HarpoFeminosHeader*)fRunHeader[YDCC])->SetPedestals(fDigi->GetPedestals(YDCC));
  }

  fRunHeader[PMM2] = new Pmm2Header();
  fRunHeader[SIMDET] = new HarpoSimHeader();

  fEvent = 0;
  // if(fSave){
  //   fTree = new TTree("fTree","Simulated HARPO events");
  //   fTree->Bronch("SimEvents", "HarpoSimEvent", &fEvent);
  // }
  if(fIonType == 2){
    TString inputFile = "testrun.dat";
    TString infile;
    if ( ! gHConfig->Lookup("Sim.IonisationFile",infile) ){
      Info("Init","Use default simulation file %s",inputFile.Data());
    } else
      inputFile = infile;
    Info("Init","Reading ionisation points from %s",inputFile.Data());
    fStr.open(inputFile.Data(), std::ifstream::in);
    if(!fStr.good()){
      Warning("Init","Could not read file %s",inputFile.Data());
      return kFALSE;
    }
  }

  return kTRUE;
}
  

Long_t      HarpoSimReader::Loop(Long_t maxevents)
{
  Long_t cnt=0;
  //maxevents = fMaxEvents;
  if(gHarpoDebug>0) Info("Loop","maxevents = %li",maxevents);
  //  HarpoSimEvent *evt;

  //fSimAnalyse->Init();

  // while (NextEvent()) {
  //   // evt = (HarpoSimEvent*)GetEvent();
  //   //fSimAnalyse->ProcessEvent(fEvent);
  //   //delete evt;
  //   if(fDoSave)
  //     fTree->Fill();
  //   fEvent->Delete();
  //   if (++cnt >= maxevents) break;
  // }


  for(Int_t i = 0; i<maxevents; i++){
    Bool_t test = NextEvent();
    if(!test) break;
    if(gHarpoDebug>1) Info("Loop","event = %d",i);
    // if(fSave){
    //   fTree->Fill();
    DeleteEvent();
    cnt++;
  }


  return cnt;

  //fSimAnalyse->Save();
  
}

Long_t      HarpoSimReader::Loop()
{
  if(gHarpoDebug>0) Info("Loop","fMaxEvents = %li",fMaxEvents);
  Long_t maxevents = fMaxEvents;
  return Loop(maxevents);
}

void      HarpoSimReader::DeleteEvent()
{
  if(gHarpoDebug>0) Info("DeleteEvent()"," ");
  std::cout << fEvents[SIMDET] << std::endl;
  fEvents[SIMDET]->Delete();
  std::cout << fEvents[XDCC] << std::endl;
  fEvents[XDCC]->Delete();
  std::cout << fEvents[YDCC] << std::endl;
  fEvents[YDCC]->Delete();
  std::cout << fEvents[PMM2] << std::endl;
  fEvents[PMM2]->Delete();

}