/* $Id: TC2R.cc 5116 2016-01-04 19:09:04Z darko $   */
/* Modified using Claire */
#include <crs2r/TC2R.h>
using namespace crs2r;

#include <crs/TBlock.h>
#include <crs/TSubBlock.h>
#include <crs/MRunHeader.h>
#include <crs/MEventHeader.h>
#include <crs/CParticle.h>
#include <crs/MParticleBlock.h>
#include <crs/MLongitudinalBlock.h>
#include <crs/MEventEnd.h>
#include <crs/MRunEnd.h>
#include <crs/MParticle.h>
#include <crs/MCherenkov.h>
#include <crs/MMuonProductionInfo.h>

#include <crsIO/TRun.h>
#include <crsIO/TShower.h>
#include <crsIO/TParticle.h>
#include <crsIO/TCherenkov.h>
#include <crsIO/TLongitudinal.h>

#include <TFile.h>
#include <TTree.h>
#include <TClonesArray.h>
#include <TObjectTable.h>

#include <iostream>



/*
 */
TC2R::TC2R() :

fFile(0),
  fShower(0),
  fRun(0),
  fCurrentRun(0),
  fCurrentShower(0),
  fParticles(0),
    fCurrentParticle(0),
  fCherenkov(0),
    fCurrentCheren(0),
  fLongitudinal(0),
    fCurrentLong(0)
{}


TC2R::~TC2R() {

  delete fCurrentRun;
  delete fCurrentShower;
    delete fCurrentParticle;
    delete fCurrentCheren;
    delete fCurrentLong;
  delete fFile;
}


void TC2R::Open(std::string filename,
		 bool thinned) {

  fFile = TFile::Open(filename.c_str(), "RECREATE");
  fFile->cd();
  fThinned = thinned;

  fCurrentShower = new crsIO::TShower();
  fShower = new TTree("sim", "simulated showers");
  fShower->Branch("EventID", &fCurrentShower->EventID);
    fShower->Branch("Energy",&fCurrentShower->Energy);
    fShower->Branch("StartingAltitude",&fCurrentShower->StartingAltitude);
    fShower->Branch("FirstTarget",&fCurrentShower->FirstTarget);
    fShower->Branch("FirstHeight",&fCurrentShower->FirstHeight);
    fShower->Branch("Theta",&fCurrentShower->Theta);
    fShower->Branch("Phi",&fCurrentShower->Phi);
    fShower->Branch("nPhotons",&fCurrentShower->nPhotons);
    fShower->Branch("nElectrons",&fCurrentShower->nElectrons);
    fShower->Branch("nHadrons",&fCurrentShower->nHadrons);
    fShower->Branch("nMuons",&fCurrentShower->nMuons);
    fShower->Branch("nParticlesWritten",&fCurrentShower->nParticlesWritten);
    fShower->Branch("nPhotonsWritten",&fCurrentShower->nPhotonsWritten);
    fShower->Branch("nElectronsWritten",&fCurrentShower->nElectronsWritten);
    fShower->Branch("nHadronsWritten",&fCurrentShower->nHadronsWritten);
    fShower->Branch("nMuonsWritten",&fCurrentShower->nMuonsWritten);
    fShower->Branch("GH_Nmax",&fCurrentShower->GH_Nmax);
    fShower->Branch("GH_t0",&fCurrentShower->GH_t0);
    fShower->Branch("GH_tmax",&fCurrentShower->GH_tmax);
    fShower->Branch("GH_a",&fCurrentShower->GH_a);
    fShower->Branch("GH_b",&fCurrentShower->GH_b);
    fShower->Branch("GH_c",&fCurrentShower->GH_c);
    fShower->Branch("GH_Chi2",&fCurrentShower->GH_Chi2);
    fShower->Branch("nPreshower",&fCurrentShower->nPreshower);
    fShower->Branch("CPUtime",&fCurrentShower->CPUtime);

  fCurrentRun = new crsIO::TRun();
  fRun = new TTree("run", "simulation run");
  fRun->Branch("RunID", &fCurrentRun->RunID);
      fRun->Branch("Date", &fCurrentRun->Date);
      fRun->Branch("Version", &fCurrentRun->Version);
      fRun->Branch("ObservationLevel", &fCurrentRun->ObservationLevel);
      fRun->Branch("ParticleID", &fCurrentRun->ParticleID);
      fRun->Branch("EnergySlope", &fCurrentRun->EnergySlope);
      fRun->Branch("EnergyMin", &fCurrentRun->EnergyMax);
      fRun->Branch("AzimuthMin", &fCurrentRun->AzimuthMin);
      fRun->Branch("AzimuthMax", &fCurrentRun->AzimuthMax);
      fRun->Branch("ZenithMin", &fCurrentRun->ZenithMin);
      fRun->Branch("ZenithMax", &fCurrentRun->ZenithMax);
    fRun->Branch("CutHadrons", &fCurrentRun->CutHadrons);
    fRun->Branch("CutMuons", &fCurrentRun->CutMuons);
    fRun->Branch("CutElectrons", &fCurrentRun->CutElectrons);
    fRun->Branch("CutPhotons", &fCurrentRun->CutPhotons);
    fRun->Branch("AtmosphereA", &fCurrentRun->AtmosphereA);
    fRun->Branch("AtmosphereB", &fCurrentRun->AtmosphereB);
    fRun->Branch("AtmosphereC", &fCurrentRun->AtmosphereC);
    fRun->Branch("BFieldX", &fCurrentRun->BFieldX);
    fRun->Branch("BFieldZ", &fCurrentRun->BFieldZ);
    fRun->Branch("EGS4", &fCurrentRun->EGS4);
    fRun->Branch("NKG", &fCurrentRun->NKG);
    fRun->Branch("Cherenkov", &fCurrentRun->Cherenkov);
    fRun->Branch("Neutrino", &fCurrentRun->Neutrino);
    fRun->Branch("Curved", &fCurrentRun->Curved);
    fRun->Branch("MuonAdditionalInfo", &fCurrentRun->MuonAdditionalInfo);
    fRun->Branch("MuonMultScatteringMoliere", &fCurrentRun->MuonMultScatteringMoliere);
    fRun->Branch("RadialRangeNKG", &fCurrentRun->RadialRangeNKG);
    fRun->Branch("TransitionEnergy", &fCurrentRun->TransitionEnergy);
    fRun->Branch("LowEHadModel", &fCurrentRun->LowEHadModel);
    fRun->Branch("HighEHadModel", &fCurrentRun->HighEHadModel);
        fRun->Branch("VersionSIBYLL_interaction", &fCurrentRun->VersionSIBYLL_interaction);
        fRun->Branch("VersionSIBYLL_crosssection", &fCurrentRun->VersionSIBYLL_crosssection);
        fRun->Branch("VersionQGSJET_interaction", &fCurrentRun->VersionQGSJET_interaction);
        fRun->Branch("VersionQGSJET_crosssection", &fCurrentRun->VersionQGSJET_crosssection);
        fRun->Branch("VersionDPMJET_interaction", &fCurrentRun->VersionDPMJET_interaction);
        fRun->Branch("VersionDPMJET_crosssection", &fCurrentRun->VersionDPMJET_crosssection);
        fRun->Branch("VersionVENUSNEXUS_crosssection", &fCurrentRun->VersionVENUSNEXUS_crosssection);
        fRun->Branch("NFLAIN", &fCurrentRun->NFLAIN);
        fRun->Branch("NFLDIF", &fCurrentRun->NFLDIF);
        fRun->Branch("NFLPI0", &fCurrentRun->NFLPI0);
        fRun->Branch("NFLPIF", &fCurrentRun->NFLPIF);
        fRun->Branch("NFLCHE", &fCurrentRun->NFLCHE);
        fRun->Branch("NFRAGM", &fCurrentRun->NFRAGM);
        fRun->Branch("HadronicThinningFraction", &fCurrentRun->HadronicThinningFraction);
        fRun->Branch("EMThinningFraction", &fCurrentRun->EMThinningFraction);
        fRun->Branch("EMThinningLimit", &fCurrentRun->EMThinningLimit);
        fRun->Branch("RadialThiningRMax", &fCurrentRun->RadialThiningRMax);
    /* note typo */
        fRun->Branch("ViewConeMin", &fCurrentRun->ViewConeMin);
    fRun->Branch("ViewConeMax", &fCurrentRun->ViewConeMax);
    fRun->Branch("nCherenkovDetectorsX", &fCurrentRun->nCherenkovDetectorsX);
    fRun->Branch("nCherenkovDetectorsY", &fCurrentRun->nCherenkovDetectorsY);
    fRun->Branch("GridCherenkovDetectorX", &fCurrentRun->GridCherenkovDetectorX);
    fRun->Branch("GridCherenkovDetectorY", &fCurrentRun->GridCherenkovDetectorY);
    fRun->Branch("LengthCherenkovDetectorX", &fCurrentRun->LengthCherenkovDetectorX);
    fRun->Branch("LengthCherenkovDetectorY", &fCurrentRun->LengthCherenkovDetectorY);
    fRun->Branch("CherenkovOutputSeparate", &fCurrentRun->CherenkovOutputSeparate);
    fRun->Branch("CherenkovBandwidthMin", &fCurrentRun->CherenkovBandwidthMin);
    fRun->Branch("CherenkovBandwidthMax", &fCurrentRun->CherenkovBandwidthMax);
    fRun->Branch("nUseCherenkovEvent", &fCurrentRun->nUseCherenkovEvent);
    fRun->Branch("CherenkovCoreX", &fCurrentRun->CherenkovCoreX);
    fRun->Branch("CherenkovCoreY", &fCurrentRun->CherenkovCoreY);
    fRun->Branch("OrientationArray", &fCurrentRun->OrientationArray);
    fRun->Branch("StepLengthFactorMultiScatter", &fCurrentRun->StepLengthFactorMultiScatter);
    fRun->Branch("Computer", &fCurrentRun->Computer);
    
    fCurrentParticle = new crsIO::TParticle();
    fParticles = new TTree("particle","simulated particles");
    fParticles->Branch("CorsikaID",&fCurrentParticle->CorsikaID);
    fParticles->Branch("ParticleID",&fCurrentParticle->ParticleID);
    fParticles->Branch("ObservationLevel",&fCurrentParticle->ObservationLevel);
    fParticles->Branch("HadronicGeneration",&fCurrentParticle->HadronicGeneration);
    fParticles->Branch("Px",&fCurrentParticle->Px);
    fParticles->Branch("Py",&fCurrentParticle->Py);
    fParticles->Branch("Pz",&fCurrentParticle->Pz);
    fParticles->Branch("x",&fCurrentParticle->x);
    fParticles->Branch("y",&fCurrentParticle->y);
    fParticles->Branch("Time",&fCurrentParticle->Time);
    fParticles->Branch("Weight",&fCurrentParticle->Weight);
    
    fCurrentLong = new crsIO::TLongitudinal();
    fLongitudinal = new TTree("long","longitudinal development");
    fLongitudinal->Branch("Depth",&fCurrentLong->Depth);
    fLongitudinal->Branch("nGammas",&fCurrentLong->nGammas);
    fLongitudinal->Branch("nElectrons",&fCurrentLong->nElectrons);
    fLongitudinal->Branch("nPositrons",&fCurrentLong->nPositrons);
    fLongitudinal->Branch("nMuons",&fCurrentLong->nMuons);
    fLongitudinal->Branch("nAntiMuons",&fCurrentLong->nAntiMuons);
    fLongitudinal->Branch("nHadrons",&fCurrentLong->nHadrons);
    fLongitudinal->Branch("nCharged",&fCurrentLong->nCharged);
    fLongitudinal->Branch("nNuclei",&fCurrentLong->nNuclei);
    fLongitudinal->Branch("nCherenkov",&fCurrentLong->nCherenkov);
    
    fCurrentCheren = new crsIO::TCherenkov();
    fCherenkov = new TTree("cherenkov","cherenkov_particles");
    fCherenkov->Branch("nPhotons",&fCurrentCheren->nPhotons);
        fCherenkov->Branch("x",&fCurrentCheren->x);
    fCherenkov->Branch("y",&fCurrentCheren->y);
    fCherenkov->Branch("u",&fCurrentCheren->u);
    fCherenkov->Branch("v",&fCurrentCheren->v);
    fCherenkov->Branch("Time",&fCurrentCheren->Time);
        fCherenkov->Branch("ProductionHeight",&fCurrentCheren->ProductionHeight);
        fCherenkov->Branch("Weight",&fCurrentCheren->Weight);
    
}


void TC2R::Write(const CREAL *data) {

  crs::TBlock block(data, fThinned);
  Write(block);

}

void TC2R::Write(const crs::TBlock &Block) {

  crs::TBlock::SubBlockConstIterator iSubBlock ;

  for(iSubBlock = Block.GetFirstSubBlock();
       iSubBlock != Block.GetLastSubBlock();
       ++iSubBlock) {

    Write(*iSubBlock);
  }
}

void TC2R::Write(const crs::TSubBlock &SubBlock) {

  if(fFile) {

    switch(SubBlock.GetBlockType()) {

    case crs::TSubBlock::eRUNH:
      fCurrentRun->Clear();
      fCurrentShower->Clear();
      fCurrentRun->AddRunHeader(SubBlock);
      fCurrentShower->AddRunHeader(SubBlock);
      break;

    case crs::TSubBlock::eEVTH:
      fCurrentRun->AddEventHeader(SubBlock);
      fCurrentShower->AddEventHeader(SubBlock);
      break;

    case crs::TSubBlock::ePARTDATA:
      AddParticleBlock(SubBlock);
      break;

    case crs::TSubBlock::eLONG:
      AddLongitudinalBlock(SubBlock);
      break;

    case crs::TSubBlock::eEVTE:
      fCurrentShower->AddEventEnd(SubBlock);

      fShower->Fill();
      break;

    case crs::TSubBlock::eRUNE:
      fRun->Fill();
      //fCurrentRun->Write();
      //fCurrentShower->Wr
      // do nothing
      break;

    case crs::TSubBlock::eNODATA:
      break;
    }

  }

}

/*
  For the inclined shower front sampling
 */


void TC2R::AddParticleBlock(const crs::MParticleBlock &Data) {

  crs::MParticleBlock::ParticleListConstIterator iEntry;
  for(iEntry = Data.ParticlesBegin();
       iEntry != Data.ParticlesEnd();
       ++iEntry) {

    switch(iEntry->GetType()) {

    case crs::eParticle:
    case crs::eNucleus:
    case crs::eMuonProductionInfo:
            *fCurrentParticle = crsIO::TParticle(*iEntry);
            fParticles->Fill();
      break;

    case crs::eCherenkov:
            *fCurrentCheren = crsIO::TCherenkov(*iEntry);
            fCherenkov->Fill();
      break;


    }
  }
}



void TC2R::AddLongitudinalBlock(const crs::MLongitudinalBlock &Data) {

  crs::MLongitudinalBlock::LongitudinalListConstIterator iLong;

  for(iLong = Data.FirstLongitudinal();
       iLong != Data.LastLongitudinal();
       ++iLong) {
      *fCurrentLong = crsIO::TLongitudinal(*iLong);
      fLongitudinal->Fill();

  }
}



void TC2R::Close() {

  // fShower->Write();
  // fRun->Write();

  if(fFile) {
    fFile->Write();
    fFile->Close();
  }
  delete fFile;
  fFile = 0;
}